840站点日照时数转日总太阳辐射Python批处理与参数标定

简介面向气象、气候变化与太阳能资源评估方向的研究者和数据分析爱好者这份资源提供了一套把全国840个气象站点观测日照时数换算为日总太阳辐射的Python实现方案。压缩包共3个文件、约275KB其中1个py脚本承担读取站点日照时数、清洗缺失与异常记录、按经验公式完成日照时数到日总太阳辐射的转换并输出结果的完整流程2个csv分别对应输入站点数据与计算后的辐射结果便于直接对照验证。已有6351人学习。读者可借此掌握批量处理多站点气象表格的写法理解日照时数与地表能量平衡、太阳能资源评估之间的换算逻辑并把它当作可复用的自动化模板替换站点或时间范围后快速产出日辐射序列为气候模型构建与资源普查提供基础数据。1. 840 个站点、十年日照时数先把日总太阳辐射算准再谈批量拿到一份 840 个气象站点的日照时数长表第一反应往往是乘一个系数了事结果同一套系数在低纬和高纬站点上能差出 20% 以上。这套 Python 脚本处理的就是日照时数转日总太阳辐射按站点纬度、日期逐日算日长 N 和大气顶层辐射 H0再用 Angstrom-Prescott 公式得到日总太阳辐射 H输出仍然是站号-日期-辐射的长表。它适合农业气象蒸散量计算、光伏资源初筛、生态模型的前处理也适合手上只有日照时数、没有辐射观测的团队。下面从公式、代码、参数标定到踩坑按能复现的顺序拆开。2. 公式选型Angstrom-Prescott、日长 N 与大气顶层辐射 H02.1 为什么不能拿日照时数直接乘固定系数日照时数 n 只说明云遮了多久不说明太阳高度角和日地距离。Angstrom-Prescott 把它写成H/H0 a b(n/N)H0 是大气顶层日辐射N 是日长。固定系数等于把 H0 和 N 都当常数纬度 20° 和 45° 的站点在夏至的 H0 能差 10 MJ/m2/day 以上。把 n 直接乘 0.5 这种写法低纬站点可能勉强接近高纬站点在冬季会明显偏高夏季又偏低。血泪经验是先把 H0 和 N 算出来再谈系数。这个公式的物理含义很直白n/N 是日照百分率a 代表全云条件下的相对辐射b 代表日照增加时辐射的斜率。ab 通常在 0.7 到 0.8 之间超过 0.9 就要怀疑系数来源或者数据单位。公式本身不复杂难的是逐日逐站把 N 和 H0 算对并且把缺测、单位、merge 这些脏活处理干净。变量含义单位备注n日照时数小时日累计值不是分钟N日长小时由纬度和日序算H0大气顶层日辐射MJ/m2/day由日地距离和太阳高度角算H日总太阳辐射MJ/m2/day目标输出a、bAngstrom 系数无量纲默认 0.25、0.50φ站点纬度度南半球为负δ太阳赤纬弧度随日序变化ωs日落时角弧度高纬要处理极昼极夜2.2 日长 N 与大气顶层辐射 H0 的计算链路逐日计算按六步走顺序不要乱。第一步取日序 doy闰年总天数 D 用 366平年用 365。第二步算太阳赤纬δ 0.409 * sin(2π * doy / D - 1.39)。第三步把纬度转弧度算x -tan(φ) * tan(δ)把 x 截断到 [-1, 1]再用ωs arccos(x)得到日落时角。第四步算日长N 24 / π * ωs。第五步算日地距离修正dr 1 0.033 * cos(2π * doy / D)。第六步算大气顶层辐射H0 (24 * 60 / π) * 0.0820 * dr * (ωs * sinφ * sinδ cosφ * cosδ * sinωs)其中 0.0820 是太阳常数单位 MJ m-2 min-1。手算一个例子纬度 35°日序 166。δ 约 0.406 radωs 约 1.878 radN 约 14.35 小时dr 约 0.968H0 约 41.6 MJ/m2/day。如果日照时数 n8.2 小时n/N0.571取 a0.25、b0.50H41.6*(0.250.50*0.571)22.3 MJ/m2/day。这个手算结果可以当作脚本的第一条单元测试。注意x不做截断时高纬站点在冬季会出现arccos定义域越界结果直接 NaN。截断不是可选项是必选项。2.3 a、b 系数怎么选默认值与本地标定没有本地辐射观测时先用默认 0.25 和 0.50 跑通全流程再谈优化。这个组合在 FAO-56 里被广泛引用适合作为基准。偏干旱站点云少、辐射强a 可以略高到 0.28 至 0.32b 保持在 0.45 至 0.50偏湿润站点云多a 可以降到 0.18 至 0.22b 略高到 0.52 至 0.58。这些范围只是筛查起点不是定论。有辐射观测的站点直接把H_obs/H0对n/N做线性回归截距是 a斜率是 b。样本至少覆盖一年四季否则冬季样本少会把系数带偏。回归后检查 a 是否在 0.1 至 0.4b 是否在 0.3 至 0.7ab 是否在 0.6 至 0.9。超出范围的先回头查单位不要急着调系数。场景ab说明无本地观测、默认0.250.50先跑通后标定偏干旱站点0.28-0.320.45-0.50云少a 略高偏湿润站点0.18-0.220.52-0.58云多b 略高有辐射观测最小二乘最小二乘用 H/H0 对 n/N 回归3. Python 脚本落地长表合并、向量化计算与 840 站点批处理3.1 输入表结构与字段约定推荐把站点表和日照时数表分开。站点表一行一个站日照表一行一个站点一天。字段名可以按实际数据改但脚本里要固定映射不要一边跑一边改列名。站号统一按字符串读避免001被读成1。日期统一转datetime64缺测码统一替换成NaN不要让-999进入公式。列名类型示例说明site_idstringA001站点唯一编号datedate2023-06-15日值日期ssdfloat8.2日照时数小时latfloat35.0纬度度lonfloat110.0经度度可选elevfloat1200海拔米可选长表的好处是 840 个站点可以一次 merge 完不需要循环读 840 个文件。数据量按 840 站 × 365 天 × 30 年估算大约 900 万行pandas 分块读完全能处理。关键是 merge 前把站点表去重merge 时加validatemany_to_one跑完断言行数没有膨胀。3.2 核心计算函数从日期和纬度生成 N、H0、H把公式封装成一个纯函数输入 DataFrame输出带N、H0、n_N、H的 DataFrame。下面这版是向量化写法840 个站点一次算不跑apply。import numpy as np import pandas as pd def add_radiation( df: pd.DataFrame, lat_col: str lat, ssh_col: str ssd, date_col: str date, a: float 0.25, b: float 0.50, clip_ratio: bool True, ) - pd.DataFrame: 按 Angstrom-Prescott 公式计算日总太阳辐射 H。 out df.copy() # 1. 日期与日序闰年用 366 d pd.to_datetime(out[date_col], errorscoerce) out[doy] d.dt.dayofyear out[days_in_year] np.where(d.dt.is_leap_year, 366, 365) # 2. 纬度、赤纬、日落时角 phi np.deg2rad(out[lat_col].astype(float)) delta 0.409 * np.sin( 2 * np.pi * out[doy] / out[days_in_year] - 1.39 ) x -np.tan(phi) * np.tan(delta) x np.clip(x, -1.0, 1.0) # 防止高纬 arccos 越界 omega_s np.arccos(x) # 3. 日长 N小时 out[N] 24.0 / np.pi * omega_s # 4. 日地距离修正与大气顶层辐射 H0MJ/m2/day dr 1 0.033 * np.cos(2 * np.pi * out[doy] / out[days_in_year]) gsc 0.0820 # MJ m-2 min-1 out[H0] ( (24 * 60 / np.pi) * gsc * dr * ( omega_s * np.sin(phi) * np.sin(delta) np.cos(phi) * np.cos(delta) * np.sin(omega_s) ) ) # 5. 日照百分率与日总辐射 n out[ssh_col].astype(float) ratio n / out[N] if clip_ratio: ratio ratio.clip(lower0, upper1) out[n_N] ratio out[H] out[H0] * (a b * ratio) return out函数里d.dt.is_leap_year决定 D 取 365 还是 366这一步不做2 月 29 日附近的 H0 会有小跳变。np.clip(x, -1, 1)处理的是高纬冬季tanφ * tanδ越界不截断就会整行 NaN。gsc0.0820的单位是 MJ m-2 min-1和24*60配套最后得到 MJ/m2/day。clip_ratio默认开启把 n/N 限制在 0 到 1如果数据里 n 大于 N先怀疑单位不要靠截断掩盖问题。参数a、b控制辐射水平默认 0.25 和 0.50。lat_col、ssh_col、date_col用来适配不同列名不要在函数里写死中文列名。clip_ratioFalse只在排查异常时用输出保留原始 n/N方便看看到底是哪些站点越界。提示如果日照时数列单位是分钟先在外面除以 60再进函数。函数内部不猜单位猜错一次后面全错。3.3 840 站点批处理入口与输出入口脚本负责读站点表、分块读日照表、merge、调函数、落盘。下面这版按chunksize分块适合几百万行以上的长表。输出只保留关键列避免中间列把文件撑大。import argparse from pathlib import Path import pandas as pd from radiation import add_radiation # 上面的函数放在 radiation.py def main() - None: parser argparse.ArgumentParser() parser.add_argument(--sites, requiredTrue) parser.add_argument(--ssh, requiredTrue) parser.add_argument(--out, requiredTrue) parser.add_argument(--a, typefloat, default0.25) parser.add_argument(--b, typefloat, default0.50) parser.add_argument(--chunksize, typeint, default500_000) args parser.parse_args() sites pd.read_csv(args.sites, dtype{site_id: str}) sites sites[[site_id, lat]].drop_duplicates(site_id) assert sites[site_id].is_unique, 站点表 site_id 不唯一 pieces [] for chunk in pd.read_csv(args.ssh, dtype{site_id: str}, chunksizeargs.chunksize): chunk[date] pd.to_datetime(chunk[date], errorscoerce) chunk chunk.dropna(subset[date]) chunk[ssd] pd.to_numeric(chunk[ssd], errorscoerce) chunk.loc[chunk[ssd] -900, ssd] pd.NA # 常见缺测码 n_before len(chunk) chunk chunk.merge( sites, onsite_id, howleft, validatemany_to_one ) assert len(chunk) n_before, merge 后行数膨胀检查站点表重复 chunk add_radiation(chunk, aargs.a, bargs.b) pieces.append( chunk[[site_id, date, ssd, N, H0, n_N, H]] ) result pd.concat(pieces, ignore_indexTrue) Path(args.out).parent.mkdir(parentsTrue, exist_okTrue) result.to_csv(args.out, indexFalse, float_format%.4f) print(frows{len(result)}, out{args.out}) if __name__ __main__: main()chunksize控制单次读入行数50 万行对普通笔记本比较稳。validatemany_to_one保证站点表一侧唯一merge 后行数不变断言失败直接停不要带着膨胀结果往下算。缺测码 -900统一置空这一步放在 merge 前避免缺测值被当成真实纬度或日照。输出列保留ssd原始值、N、H0、n_N、H方便后面抽查。命令行这样跑python run_radiation.py \ --sites data/sites.csv \ --ssh data/ssd_long.csv \ --out output/H_daily.csv \ --a 0.25 --b 0.50 \ --chunksize 500000参数--sites和--ssh分别指向站点表和日照长表--out是输出路径。--a、--b先给默认值后面标定完再换。--chunksize按内存调数据量大就降到 20 万数据量小可以加到 100 万。跑完先看rows是否等于输入行数不等就先查 merge 和缺测过滤。4. 参数标定与校验a/b 系数、单位换算和异常值拦截4.1 有辐射观测时用最小二乘标定 a/b有辐射观测的站点把观测辐射统一到 MJ/m2/day再用H_obs/H0对n/N做线性回归。回归前过滤n/N在 0 到 1、H_obs/H0在 0.05 到 1.2 的样本。np.polyfit返回的是斜率和截距斜率对应 b截距对应 a顺序不要写反。import numpy as np import pandas as pd def fit_angstrom( df: pd.DataFrame, h_obs_col: str H_obs, h0_col: str H0, n_col: str ssd, N_col: str N, ) - tuple[float, float]: 用 H_obs/H0 对 n/N 线性回归返回 a, b。 x (df[n_col] / df[N_col]).clip(0, 1) y df[h_obs_col] / df[h0_col] mask ( np.isfinite(x) np.isfinite(y) (x 0) (x 1) (y 0.05) (y 1.2) ) if mask.sum() 30: raise ValueError(有效样本不足 30先查单位或缺测) slope, intercept np.polyfit(x[mask], y[mask], 1) return float(intercept), float(slope)mask里同时过滤 x 和 y缺测、负值、异常高值都会被挡掉。样本少于 30 直接报错不要用几个点硬拟合。返回顺序是a, b即截距在前、斜率在后。标定完把ab打出来正常在 0.6 到 0.9如果a是负数或者b大于 1先查观测辐射单位是不是 W/m2 没转成 MJ/m2/day。4.2 无观测时的系数迁移与合理性检查没有辐射观测的站点用相近气候区的标定系数或者先用默认值。迁移系数时不要跨气候带硬套至少按干旱/湿润、低纬/高纬分两组。跑完做四类检查n_N是否在 0 到 1H0是否在 0 到 45 MJ/m2/dayH/H0是否在 0.05 到 0.9H是否为非负。任何一条不满足先标记不要直接删。def flag_radiation(df: pd.DataFrame) - pd.DataFrame: out df.copy() out[flag_nN] ~out[n_N].between(0, 1) out[flag_H0] ~out[H0].between(0, 45) out[flag_ratio] ~(out[H] / out[H0]).between(0.05, 0.9) out[flag_H] out[H] 0 out[flag_any] out[ [flag_nN, flag_H0, flag_ratio, flag_H] ].any(axis1) return out检查项合理范围异常处理n/N0-1查日照时数单位再截断H00-45 MJ/m2/day查纬度符号、日序、闰年H/H00.05-0.9查 a/b 系数和观测单位H非负查缺测码是否参与计算N0-24 小时高纬极昼极夜单独标记4.3 输出单位与常见误用Angstrom-Prescott 算出来的是 MJ/m2/day很多模型要 W/m2。换算系数是1e6 / 86400约 11.574。也就是H_Wm2 H_MJ * 11.574。不要把 MJ 和 cal/cm2 混着用1 cal/cm2 约 0.04184 MJ/m2差两个数量级。# MJ/m2/day 转 W/m2 df[H_Wm2] df[H] * 1e6 / 86400 # 如果原始观测是 cal/cm2/day先转 MJ/m2/day df[H_obs_MJ] df[H_obs_cal] * 0.04184单位转换放在输出前不要在中途来回转。常见误用有三种把日照时数分钟当小时n/N 直接翻倍把 H0 当 H 输出结果偏高把 W/m2 当 MJ/m2/day 喂给蒸散模型结果偏小两个数量级。每次换数据源先把单位列出来比事后调系数省事。5. 避坑与排查日照时数转辐射最常翻车的五处5.1 n/N 大于 1先查日照时数单位再查日长现象输出里n_N出现 1.3、2.0H甚至超过H0。原因通常不是公式错而是日照时数单位混了有的表按分钟存有的按小时存还有的用 0.1 小时为单位。另一个原因是日长 N 算错比如纬度用了弧度又转了一次。解决先在数据字典里确认ssd单位分钟就除以 60进公式前把负值和缺测码置空n_N截断到 0 到 1 只作为兜底不要用它掩盖单位问题。跑完统计n_N 1的站点数超过 1% 就停下来查源表。5.2 闰年日序用 3652 月 29 日附近 H0 会跳现象闰年 2 月 29 日前后H0出现不连续或者 12 月 31 日的赤纬和往年对不上。原因赤纬公式里的 D 固定写成了 365闰年应该用 366。解决用d.dt.is_leap_year生成days_in_year再代入赤纬和日地距离修正。不要用doy/365一把梭。这个错误很小但在时间序列里会变成一个尖刺后面做月均值时容易被当成异常值删掉。5.3 高纬站点 arccos 出现 NaNtanφ tanδ 越界现象高纬站点冬季N和H0整行 NaN低纬站点正常。原因x -tanφ * tanδ超出 [-1, 1]arccos定义域越界。解决计算前x np.clip(x, -1.0, 1.0)。截断后ωs接近 0 或 π对应极夜或极昼日长接近 0 或 24 小时。极昼极夜不是脏数据应该在输出里加标记不要直接删。南半球纬度给负值tanφ和sinφ会自动变号不要手动改公式。5.4 merge 后行数膨胀站点表站号重复现象输出行数比输入多或者同一个站号同一天出现两条记录。原因站点表site_id有重复merge 变成多对多。解决merge 前drop_duplicates(site_id)merge 时加validatemany_to_onemerge 后断言行数等于输入行数。站号按字符串读避免A001和a001被当成两个站。如果站点表里一个站有多行经纬度先确认哪一行是有效坐标不要直接全量 merge。5.5 缺测码 -999 参与计算H 出现负值或巨大值现象H出现 -200 或者 1e6 量级的值n_N严重偏离。原因缺测码-999、-9999没有替换直接进了公式。解决读入后统一replace([-999, -9999], np.nan)或者用pd.to_numeric(errorscoerce)加阈值过滤。缺测行不要删保留ssd为空、H为空加flag_missing列后面统计缺测率。缺测率超过 20% 的站点标定系数时不要用。6. 进阶用法分块并行与日常任务化别让 840 站跑成黑匣子840 个站点如果每个站点 30 年日值大约 900 万行单进程 pandas 也能跑但每次调参都全量重算会很拖。我一般先用--sample-sites 10抽 10 个站点跑通再用分块并行。分块不要按行随机切按站点切保证一个站点的日期连续。参数建议如下数据规模分块方式建议小于 100 万行单进程chunksize 50 万100 万到 1000 万行按站点分块8 进程左右超过 1000 万行按站点加年份分区写 parquet别全量 CSV按站点并行的骨架很短每个进程拿到一个站点的子表调同一个add_radiation最后 concat。注意进程数不要超过物理核数840 个站点分 8 组比开 840 个进程稳。from multiprocessing import Pool def run_one(item): site_id, sub item return add_radiation(sub, a0.25, b0.50) if __name__ __main__: groups list(df.groupby(site_id)) with Pool(8) as pool: parts pool.map(run_one, groups) result pd.concat(parts, ignore_indexTrue)日常任务化时把a、b、chunksize、sample_sites写进配置文件不要每次改脚本。输出加一列run_id记录这次用的系数和日期后面对比不同版本结果时能追。最省事的验证方法是抽 10 个站点用手算的 N、H0 和脚本对一遍误差在 0.1 以内再放全量。从那以后我每次全量前都强制走一遍抽样对账确认H0在 0 到 45、n_N在 0 到 1、H非负再让 840 个站点跑到底。希望帮到你。本文还有配套的精品资源点击获取