简介面向气象、环境科学与太阳能资源评估方向的研究者及数据分析爱好者提供一套将全国840个气象站点观测日照时数SSD批量换算为日总太阳辐射SR的Python处理方案解决逐站点手工计算效率低、口径不统一的问题属于可直接套用的实用脚本类资源。压缩包共3个文件含2个CSV数据表分别对应输入日照时数与示例输出结果便于对照验证和1个Python主脚本整体约275KB轻量易部署。该资源已有6348人学习下载说明其在同类气象数据处理场景中具备较高的参考认可度。脚本围绕读取、预处理、转换计算与结果输出组织流程可据此理解经验公式中地理位置、季节与大气透明度等因素的接入方式并借助pandas批量处理缺失值、异常值与日期格式问题还可扩展matplotlib绘图与并行加速。读者既能获得可复用的换算模板也能通过输入输出样例快速核对计算结果适合作为气象数据批处理的入门与练手案例。1. 840 个站点的日照时数要变成辐射值先得算清分母某开发者手上有一份日照时数的日值数据八年的记录准备喂给一个蒸散模型结果卡在最后一步模型要的是日总太阳辐射单位 MJ/m²数据里只有每天日照了几个小时。他先试了乘一个固定系数跟辐射观测比对夏天差了将近三成冬天又反过来偏低。这不是参数没调好是方法从根上就不对日照时数只反映晴天占了多少同样的 8 小时高纬度夏天只占可照时数的五成低纬度却接近七成落到地面的日总太阳辐射压根不是一回事。这份 Python 脚本处理的就是这一步读入 840 个气象站点的站点表与日照时数日值文件逐站逐日先算天文辐射和可照时数再用 Ångström–Prescott 公式换算成日总太阳辐射输出带单位标注的 CSV。适合手上只有日照时数、又急需辐射驱动场的场景比如蒸散计算、作物模型、光伏潜力粗估。难点不在代码量在于几个中间量的单位和边界条件这篇就按「公式 → 数据 → 代码 → 排错 → 改造」的顺序拆。2. 公式选型Ångström–Prescott、天文辐射与 a、b 系数的取值2.1 固定系数法错在哪分母随纬度和季节变把一个站点的日照时数直接乘 0.5 或者查一张「月份-系数」表是最省事也最容易出错的做法。问题在于日照时数本身不是能量它只是「太阳实际露脸多久」的时长统计要变成能量必须知道这段时间里大气层顶本来能接收多少辐射。这个「本来能接收多少」就是天文辐射 Ra它只跟日地距离、太阳赤纬和站点纬度有关跟云、气溶胶、海拔统统无关。所以换算链路是两步先把 Ra 和可照时数 N 算出来再用实测日照时数 n 除以 N 得到一个 0 到 1 之间的晴空比例最后乘上经验系数。Ångström–Prescott 形式最简单Rs Ra × (a b × n/N)其中 Rs 是日总太阳辐射Ra 是天文辐射n 是实测日照时数N 是天文可照时数a、b 是站点或区域的经验系数。这个式子的物理含义很直白a 代表全阴天时仍能到达地面的那部分散射辐射b 代表晴天比例每增加一个单位辐射增加多少。全阴天 n0 时 Rs a×Ra全晴天 nN 时 Rs (ab)×Ra两个端点都能对上物理直觉比拍脑袋乘系数靠谱得多。2.2 天文辐射与可照时数三个中间量必须先算对Ra 的标准算法FAO-56 里那套分三步。第一步算日地相对距离倒数 dr因为地球绕日轨道是椭圆一年里日地距离变化约 ±3.3%冬天近、夏天远对高纬度站点影响不能忽略。第二步算太阳赤纬 δ也就是太阳直射点纬度它决定了当天正午太阳高度角。第三步算日落时角 ωs由纬度与赤纬共同决定可照时数 N 就是 ωs 的两倍再换算成小时。三个量的表达式如下doy 是日序数1 到 365 或 366dr 1 0.033 × cos(2π × doy / 365) δ 0.409 × sin(2π × doy / 365 − 1.39) ωs arccos(−tan φ × tan δ) Ra (24×60/π) × Gsc × dr × [ωs × sinφ × sinδ cosφ × cosδ × sinωs] N 24/π × ωs这里 Gsc 是太阳常数取 0.0820 MJ m⁻² min⁻¹前面那个 24×60/π 是把分钟辐射累积成全天并把角度换成弧度的组合系数。φ 是站点纬度北纬为正、南纬为负必须换成弧度再进三角函数。这套公式在 FAO-56 灌溉手册里有完整推导农业气象和蒸散计算领域基本都认它。符号含义单位典型量级Gsc太阳常数MJ m⁻² min⁻¹0.0820dr日地相对距离倒数无量纲0.97 ~ 1.03δ太阳赤纬rad−0.41 ~ 0.41φ站点纬度弧度rad取决于站点ωs日落时角rad0 ~ πRa天文辐射MJ m⁻² d⁻¹5 ~ 42N可照时数h8 ~ 16中低纬注意 ωs 那一行arccos 的入参是 −tanφ×tanδ在中高纬度夏季可能小于 −1冬季可能大于 1这就是极昼极夜的数学表现。直接扔给 arccos 会返回 NaN批量跑 840 个站点时只要有那么几个高纬站整张输出表就会莫名其妙出现空值这个坑在第 5 章还会展开。2.3 a、b 系数怎么定先用默认值再谈标定a、b 不是普适常数它们随气候区、海拔、气溶胶状况变化。国内多数台站拟合出来的 a 落在 0.18 到 0.25 之间b 落在 0.50 到 0.60 之间比较常用的缺省组合是 a0.20、b0.55。如果你的研究区没有现成的标定值用这组缺省值先跑通、看量级对不对比纠结系数更有意义。真要标定思路是拿同一站点、同一时段的实测辐射 Rs_obs 和算出来的 Ra、n/N 做线性回归以 n/N 为自变量、Rs_obs/Ra 为因变量截距就是 a斜率就是 b。至少要有三年逐日数据回归的 R² 一般能到 0.85 以上达不到就说明该站的日照时数记录有问题或者辐射观测本身有系统性偏差。分月标定是进阶做法因为雨季和旱季的云型不同全年一组系数在季风区会明显低估雨季、高估旱季。2.4 先拿两个特殊日验算确认公式没抄错写完 Ra 函数不要急着跑 840 个站。挑一个纬度 40°N 的站点、日序数 172夏至前后手算或用独立脚本验一遍量级对不对一眼就能看出来。import numpy as np # 只验证单点纬度 40 度日序数 172夏至前后 lat, doy 40.0, 172 phi np.radians(lat) dr 1 0.033 * np.cos(2 * np.pi * doy / 365) delta 0.409 * np.sin(2 * np.pi * doy / 365 - 1.39) ws np.arccos(-np.tan(phi) * np.tan(delta)) ra (24 * 60 / np.pi) * 0.0820 * dr * ( ws * np.sin(phi) * np.sin(delta) np.cos(phi) * np.cos(delta) * np.sin(ws) ) print(fRa{ra:.2f} MJ/m2/d, N{24/np.pi*ws:.2f} h) # 期望输出Ra≈41.8 MJ/m2/d, N≈14.8 h参数说明lat是纬度度数doy是日序数dr用 365 做分母是 FAO 的标准写法闰年换 366 只影响小数第三位可以忽略。这段代码的价值不在计算在于给你一个可对照的锚点40°N 夏至 Ra 应该在 41 到 42 MJ m⁻² d⁻¹ 之间可照时数 14.8 小时左右。如果跑出来是 4.18八成是 Gsc 少乘了一个 10如果是 418那是把 24×60/π 重复乘了。再从反面验一个同一纬度、日序数 355冬至前后Ra 应该在 12 到 13 MJ m⁻² d⁻¹N 约 9.2 小时。夏冬两个端点的比值大约 3.4 倍这个比例关系比绝对数值更能说明公式写对了没有。3. 输入数据与目录结构站点表、日照日值文件和输出约定3.1 站点表三个字段必须能对上840 个站点的元信息至少要包含站号、纬度、经度海拔可选但建议留着后面做气压订正或者质量筛查用得上。站号一定要按字符串读不要让它自动变成整数——有些站号带前导零pd.read_csv默认会把它转成 int一读进来就少了一位跟日照文件里的站号对不上最后的结果是全部站点都被跳过而日志里只打印一行「缺文件」排查起来非常费劲。站点表的列名建议统一成station, lat, lon, elev。如果原始表里是「经度」「纬度」「台站名」这种中文列名在读入后立刻改掉别让中文列名流进后面的函数itertuples取属性时中文名会变成位置索引代码可读性直接崩掉。海岸线附近的站点要特别留意符号西经、南纬必须转成负数国内站点虽然只涉及东经北纬但如果以后要把代码迁到其他区域符号约定最好一开始就定死。还有一类常见情况是站点表里经纬度用度分秒表示比如116°28′这种字符串。这种列 pandas 会读成 object直接进np.radians就报错。写一个转换函数用正则把度分秒拆出来再合并成十进制度转换完打印前五行人工核对一遍比后面查错省事。3.2 日照时数日值文件格式和缺测码先摸清楚日值数据集常见的组织方式是每站一个文本文件或者全国一个文件按站号排序。文本里每一行是一条日记录字段大致是站号、年、月、日、日照时数、质量控制码。分隔符可能是空格、制表符也可能是定宽格式。这三种情况在 pandas 里对应三种读法选错了不会报错只会读出一列全是 NaN 或者一堆错位的字段。空格或制表符分隔pd.read_csv(f, sepr\s, headerNone)定宽格式pd.read_fwf(f, widths[...], headerNone)widths 要按实际字符位置数出来有表头且字段名规整直接read_csv加dtype{station: str}日照时数的缺测码是最容易出问题的地方。日值数据集里常用 32766 表示缺测也有用 999999 或者 −9999 的不同批次、不同要素的约定可能不一样。处理原则很简单先df[ssd].describe()看一眼最大值如果出现 32766 这种明显不是小时数的值就是缺测码确认后统一替换成pd.NA再参与计算。千万别用「大于 24 就截断成 24」这种想当然的处理那会把缺测硬生生变成「全天满日照」辐射值直接顶到天文辐射上限。质量控制码也别急着丢。质控码为 0 或 1 通常表示数据正常或可疑但可用2 以上往往意味着该条记录有问题。初跑阶段可以保留这一列后面统计一下各质控码的占比如果某站的异常码比例超过 10%这个站的输出单独标记出来别跟正常站混在一张表里。3.3 目录布局与命名约定840 个站点的数据量不算大单站八年逐日也就三千行左右合计两百多万行一台普通笔记本十几秒能跑完。但目录结构如果不提前定好跑到一半缺一个文件、多一个重复站号重跑的成本比想象中高。我习惯按下面这样摆project/ ├── data/ │ ├── stations.csv # 站号,lat,lon,elev │ └── ssd/ # 每站一个 txt文件名即站号 │ ├── 5xxxx.txt │ └── ... ├── core/ │ └── radiation.py # 纯计算函数不碰 IO ├── run_all.py # 批量调度 └── output/ └── rs_daily.csvcore/radiation.py里只放数学函数输入是数组、输出是数组不读文件也不写文件。这样做的好处是单测好写改系数的时候不会牵一发动全身。run_all.py负责遍历站点、读文件、调函数、拼结果。文件名直接用站号省掉在站点表里再维护一个文件路径列也避免了路径里带空格或中文导致的编码问题。3.4 先跑一遍数据体检再进计算正式计算之前写一个十几行的体检脚本把三类问题提前暴露站点表有没有重复站号、日照文件有没有缺失、日照时数的分布有没有异常。这一步花五分钟能省掉后面几个小时的翻车排查。from pathlib import Path import pandas as pd sites pd.read_csv(data/stations.csv, dtype{station: str}) print(站点数:, len(sites), 重复站号:, sites[station].duplicated().sum()) print(纬度范围:, sites[lat].min(), sites[lat].max()) ssd_dir Path(data/ssd) missing [s for s in sites[station] if not (ssd_dir / f{s}.txt).exists()] print(缺文件站点数:, len(missing), missing[:10]) # 抽查一个站看日照时数的分布 sample pd.read_csv(ssd_dir / f{sites[station].iloc[0]}.txt, sepr\s, headerNone, names[station, year, month, day, ssd, qc]) print(sample[ssd].describe()) print(疑似缺测码出现次数:, (sample[ssd] 32766).sum())逻辑上分三段第一段核对站点表自身的一致性重复站号和纬度越界是最常见的两类脏数据第二段用集合差找出缺文件的站注意这里没做异常捕获文件不存在就是不存在早暴露早处理第三段抽一个站看ssd的统计量max如果是 32766 就说明缺测码没处理min如果是负数说明有异常记录。参数上dtype{station: str}是必须的不写这一句站号读成整数后跟文件名比对会全部失败。names里的字段顺序要跟原始文件严格对齐字段数不匹配时 pandas 会报ParserError报错信息里会给出实际列数照着改就行。这一步跑完如果缺文件站点数是 0、ssd的 max 在 16 以内、min 不小于 0就可以放心进核心计算了。4. 核心模块拆解三个函数走完从小时数到 MJ/m² 的链路4.1 函数一extraterrestrial_radiation 算 Ra 和 N这个函数是整个链路的地基输入日序数和纬度输出天文辐射与可照时数两个数组。写成向量化形式一次能算完一个站八年的所有日比逐行循环快两个数量级。import numpy as np GSC 0.0820 # 太阳常数, MJ m-2 min-1 MJ_TO_CAL 23.89 # 1 MJ m-2 23.89 cal cm-2 def extraterrestrial_radiation(doy, lat_deg): 逐日天文辐射 Ra(MJ m-2 d-1) 与可照时数 N(h)。 doy : 日序数标量或数组1~366 lat_deg : 纬度北正南负单位度 doy np.asarray(doy, dtypefloat) phi np.radians(np.asarray(lat_deg, dtypefloat)) dr 1.0 0.033 * np.cos(2.0 * np.pi * doy / 365.0) delta 0.409 * np.sin(2.0 * np.pi * doy / 365.0 - 1.39) x -np.tan(phi) * np.tan(delta) x np.clip(x, -1.0, 1.0) # 极昼极夜保护避免 arccos 出 NaN ws np.arccos(x) ra (24.0 * 60.0 / np.pi) * GSC * dr * ( ws * np.sin(phi) * np.sin(delta) np.cos(phi) * np.cos(delta) * np.sin(ws) ) n_day 24.0 / np.pi * ws return ra, n_day逐段说明。np.asarray把标量和列表统一成数组这样后面全是逐元素运算调用方传一个数或者传一整列都不会出错。dr和delta两行是纯周期函数分母用 365 是 FAO 的标准写法闰年用 366 会让结果差千分之几对辐射量级没有实质影响不必为此加分支。np.clip(x, -1.0, 1.0)这一句是整段代码里最容易被忽略、也最容易出事的地方。纬度超过 66.5° 的站点在夏季会出现 x −1冬季出现 x 1此时数学上对应极昼和极夜物理上 ωs 应该取 π 和 0。不截断的话 arccos 返回 NaNNaN 会沿着乘法一路污染整行结果最后的 CSV 里就出现莫名其妙的空值而你还以为只是某个站的数据缺了几天。返回的两个数组单位分别是 MJ m⁻² d⁻¹ 和小时后面所有计算都基于这两个量单位不要在中途改来改去。4.2 函数二sunshine_to_radiation 做换算和裁剪拿到 Ra 和 N第二步就是套 Ångström–Prescott。这个函数很短但两个边界处理不能省。def sunshine_to_radiation(n_hours, ra, n_day, a0.20, b0.55): Ångström-Prescott 换算Rs Ra * (a b * n/N)。 n_hours : 实测日照时数(h)缺测应已在外部剔除 ra : 天文辐射(MJ m-2 d-1) n_day : 可照时数(h) a, b : 经验系数缺省 0.20 / 0.55 n np.asarray(n_hours, dtypefloat) ratio np.divide(n, n_day, outnp.zeros_like(n, dtypefloat), wheren_day 0) ratio np.clip(ratio, 0.0, 1.0) # n 超过 N 或为负归到边界 return ra * (a b * ratio)np.divide带where参数的写法是为了处理 n_day 为 0 的情况——极夜时 N0直接除会得到 inf 或者警告。out先填 0where只对 N0 的位置做除法剩下的保持 0逻辑上就是「极夜时晴空比例为 0」虽然极夜里 Ra 本身也接近 0结果无论如何都是 0但至少不会污染下游。np.clip(ratio, 0, 1)处理的是观测记录问题。日照时数记录仪换型、人工观测误记都可能出现 n 略大于 N 的情况比如 N14.8 而记录写了 15.1。这种记录不裁掉的话a b×1.02会算出比天文辐射还大的地面辐射物理上不可能。裁到 1 相当于把它当成全天满日照是误差最小的处理方式。如果这种记录占比很高比如超过 3%就要回头查数据源了光在代码里裁是不够的。系数 a、b 作为默认参数放在签名里调用时可以逐站传入标定值不用改函数体。这也是把计算逻辑单独抽成一个模块的好处。4.3 函数三process_station 负责单站装配和落盘字段前两个函数是纯数学第三个函数负责把一年的数据装配起来生成最终要写进 CSV 的字段。这里要一次性把单位都算好别等到后面再补。import pandas as pd def process_station(meta, df, a0.20, b0.55): 单站装配输入站点元信息和日值表返回带辐射列的表。 df df.sort_values([year, month, day]).reset_index(dropTrue) dates pd.to_datetime(df[[year, month, day]]) doy dates.dt.dayofyear.to_numpy() ra, n_day extraterrestrial_radiation(doy, meta[lat]) rs sunshine_to_radiation(df[ssd].to_numpy(), ra, n_day, a, b) out df[[station, year, month, day]].copy() out[doy] doy out[ra_mj] np.round(ra, 3) out[n_day_h] np.round(n_day, 2) out[ssd_h] df[ssd].to_numpy() out[rs_mj] np.round(rs, 3) out[rs_cal] np.round(rs * MJ_TO_CAL, 1) # cal cm-2 d-1 out[rs_wm2] np.round(rs / 0.0864, 1) # 日均辐射通量 return outpd.to_datetime那两行是为了拿到日序数。有的数据集只给年月日不给日序数用 datetime 转是最稳的办法闰年 2 月 29 日会自动算成 60年底 12 月 31 日是 366不用自己写月份天数表。自己手写[0,31,59,...]累加数组碰上闰年就整体错一天这种错误在小样本上看不出来跑完八年才发现整条序列偏移非常难查。输出字段里刻意保留了三个单位的辐射rs_mj是主流模型要的 MJ m⁻² d⁻¹rs_cal是老一辈气象资料和部分农业手册用的 cal cm⁻² d⁻¹rs_wm2是日均辐射通量跟再分析产品的单位对得上。多算两列的成本几乎为零但能省掉后面反复换算、反复怀疑单位的麻烦。np.round只是为了让 CSV 好看如果下游还要做精确计算建议去掉 round 直接用原值。meta[lat]是单站的纬度df是已经去掉缺测、按日期排序的日值表。函数内部没有做任何 IO也不改传入的 df返回的是新表这样批量循环里出问题时容易定位到具体哪个站。4.4 批量调度840 个站点怎么稳定跑完批量脚本要做的事很机械读站点表、遍历、读文件、剔缺测、调函数、拼结果。真正需要注意的是异常处理和日志跑两百多万行的时候中间某一行报错不能把整个任务带崩。from pathlib import Path import pandas as pd from core.radiation import process_station ssd_dir Path(data/ssd) sites pd.read_csv(data/stations.csv, dtype{station: str}) frames, skipped [], [] for row in sites.itertuples(indexFalse): f ssd_dir / f{row.station}.txt if not f.exists(): skipped.append((row.station, 缺文件)) continue raw pd.read_csv(f, sepr\s, headerNone, names[station, year, month, day, ssd, qc]) raw[ssd] pd.to_numeric(raw[ssd], errorscoerce) raw.loc[raw[ssd] 32766, ssd] pd.NA # 缺测码剔除 raw raw.dropna(subset[ssd]) if raw.empty: skipped.append((row.station, 全为缺测)) continue frames.append(process_station({station: row.station, lat: row.lat}, raw)) result pd.concat(frames, ignore_indexTrue) result.to_csv(output/rs_daily.csv, indexFalse, float_format%.3f) print(完成站点:, len(frames), 跳过:, skipped[:5])itertuples比iterrows快而且拿到的row.station、row.lat就是列名属性比row[station]好读。缺文件不抛异常记到skipped里继续下一个这样一趟跑完你能拿到完整的跳过清单而不是跑到第 500 个站崩掉重来。pd.to_numeric(..., errorscoerce)是第二道保险万一某个文件里混进了非数字字符这一句会把它变成 NaN 而不是让整列变成 object。raw.loc[raw[ssd] 32766, ssd] pd.NA处理缺测码如果你的数据源用的是 999999改这个阈值就行。最后float_format%.3f控制输出精度避免 CSV 里出现一长串浮点尾数。840 个站、八年逐日总行数大约两百多万pd.concat一次性拼完内存压力不大普通笔记本 8G 内存够用。如果站数翻到几千或者年份拉长到三十年就得改成分批写文件再合并或者直接上polars。不过那是另一个话题了当前规模没必要提前优化。5. 常见问题排查批量跑 840 站时最容易翻车的五个点5.1 结果整体偏大或偏小两成现象所有站点的辐射值跟观测比对系统性偏高或偏低 15% 到 25%各站方向一致不是随机误差。原因大概率是单位混用。常见的有三种一是 Ra 用了 MJ 但输出标成 cal/cm² 没换算差了 23.89 倍二是把日均辐射通量 W/m² 直接当成了日累积量 MJ/m²/d两者差 0.0864 倍三是系数 a、b 用的是「a0.25、b0.50」这类老文献值但公式写成了a b×n/N之外的形式比如漏了 b 只乘了 n。解决办法是先算一个全晴天样本对照取某站夏季日照时数接近 N 的几天Rs 应该略小于 Ra比值在 0.70 到 0.75 之间。如果算出来比值是 0.03 或者 17就是单位问题如果比值在 0.5 左右那是系数问题。这个对照法比逐项查代码快得多我一般先做这个再动代码。5.2 高纬度站点算出 NaN现象输出表里部分站点的rs_mj是空值集中在纬度偏高的站而且往往集中在夏季或冬季的某几天。原因就是 2.2 节提到的arccos越界。纬度接近或超过极圈时−tanφ×tanδ 的绝对值会大于 1numpy 返回 NaN并且给出 RuntimeWarning。如果日志被刷屏这个警告很容易被忽略。解决就是在 arccos 之前加np.clip(x, -1.0, 1.0)。截断后极昼时 ωsπ、N24 小时极夜时 ωs0、N0物理上是对的。如果站点确实在极圈内还要考虑极昼期间日照时数记录本身是否有效这时候更稳妥的做法是对 N24 的日子单独打标记而不是直接输出。5.3 日照时数大于可照时数现象算出来的 Rs 比 Ra 还大或者n/N大于 1某些日子辐射值明显不合理。原因有两类。一类是观测记录问题仪器换型、人工读数误记日照时数记成了 15 小时而当天 N 只有 14.2。另一类是公式里的 N 算错了比如纬度符号搞反北纬传成了负值算出来的可照时数跟实际对不上。解决办法分两步先在sunshine_to_radiation里用np.clip(ratio, 0, 1)兜底保证输出不会超过 Ra然后统计n N的记录占比超过 3% 就说明数据源有问题要回去查那个站的原始记录。单纯靠 clip 掩盖过去等于把数据质量问题藏起来了后面做趋势分析会出偏差。5.4 缺测码被当成真实值参与计算现象某几天的辐射值巨大可能是正常值的十几倍集中在特定站点、特定时段。原因就是 32766 这类缺测码没剔除直接进了乘法。32766 小时除以 N 得到几千的比值再乘 Ra结果能到几万 MJ/m²一眼就能看出来不对但如果只做汇总统计没逐行看很容易被平均值掩盖。解决是在读取阶段就把缺测码替换成 NA 并 drop 掉不要留到计算阶段再处理。同时建议在输出里加一列有效标记记录该行是实测还是缺测插补后面做时间序列分析时能区分开。缺测率高的站点比如超过 20%建议整站不进主表单独输出到另一个文件。5.5 站点表经纬度是度分秒字符串现象读站点表时lat列是 object 类型进np.radians直接抛 TypeError或者被 pandas 静默转成 NaN 后算出全站一样的辐射值。原因是原始站点表用「116°28′」这类格式记录经纬度read_csv无法识别成数字。这种情况在整理过的数据集里很常见尤其是从网页或者 PDF 转录过来的表。解决办法是写一个解析函数用正则(\d)[°度](\d)[′分]?拆出度和分再合并成度 分/60。转换完打印前五行跟原始表对照确认没有把纬度和经度列搞反。经纬度搞反是最隐蔽的错代码不报错结果看起来也「像那么回事」但整个中国的站点辐射分布会完全错位。我一般会在体检脚本里加一句检查国内站点的经度应该在 70 到 140 之间纬度在 15 到 55 之间超出范围就报警。6. 二次开发把系数、单位、时间分辨率换成你自己项目要的跑通只是起点真正用起来总要改点什么。这个项目里最常改的三个地方是经验系数、输出单位和时间尺度每个都有对应的改法不用动核心公式。改造点改哪里怎么改验证方式分区标定 a、bprocess_station的入参站点表里加a_coef、b_coef两列循环时传进去与实测辐射回归看 R² 和偏差分月系数sunshine_to_radiation内部系数改成 12 元素数组按month索引逐月偏差是否收敛输出 cal/cm²输出字段保留rs_cal或改MJ_TO_CAL精度与历史气象资料对照月尺度汇总run_all.py之后按站号年月 groupby 求均值注意是日均值不是月总量月总量除以天数应接近日均换数据源读取段改分隔符、列名、缺测码阈值体检脚本先跑一遍分月系数是我改得最多的一项。季风气候区全年一组系数雨季低估、旱季高估误差能到 15%。把 a、b 换成 12 个值的数组函数里用月份做索引改动不到十行def sunshine_to_radiation_monthly(n_hours, ra, n_day, month, a_arr, b_arr): 分月系数的 Ångström-Prescotta_arr/b_arr 为长度 12 的数组。 idx np.asarray(month, dtypeint) - 1 # 月份 1~12 映射到 0~11 a np.asarray(a_arr, dtypefloat)[idx] b np.asarray(b_arr, dtypefloat)[idx] ratio np.clip(np.divide(n_hours, n_day, outnp.zeros_like(n_hours, dtypefloat), wheren_day 0), 0.0, 1.0) return ra * (a b * ratio)这里的idx一定要减 1月份从 1 开始而数组从 0 开始忘了减就会取到下一个月的系数冬天用成春天的值误差不大不小正好让你查不出来。系数数组的长度也要校验传进来 11 个值的时候索引 11 会越界加一句assert len(a_arr) 12能省不少事。另一个值得做的改造是加一层结果校验。算完不要直接拿去用挑三到五个有辐射观测的站点把算出来的 Rs 和实测值放一起看散点和偏差。我习惯算三个指标平均偏差MBE、均方根误差RMSE、相关系数。日照时数换算的辐射日尺度 RMSE 一般在 2 到 3 MJ m⁻² d⁻¹相关系数 0.9 以上算合格。如果 RMSE 超过 4先查缺测处理再查系数最后才怀疑公式本身。血泪经验是别在同一个脚本里既读数据又算公式又画图改一个系数要重跑整条链路等得人心焦。计算模块和可视化模块彻底分开中间用 CSV 交换改系数只需要重跑计算那一步几分钟出结果。从那以后我每次接这类换算任务都强制先把「纯函数 独立批量脚本 体检脚本」三件套搭好再动手写公式希望帮到你。本文还有配套的精品资源点击获取
简介面向气象、环境科学与太阳能资源评估方向的研究者及数据分析爱好者提供一套将全国840个气象站点观测日照时数SSD批量换算为日总太阳辐射SR的Python处理方案解决逐站点手工计算效率低、口径不统一的问题属于可直接套用的实用脚本类资源。压缩包共3个文件含2个CSV数据表分别对应输入日照时数与示例输出结果便于对照验证和1个Python主脚本整体约275KB轻量易部署。该资源已有6348人学习下载说明其在同类气象数据处理场景中具备较高的参考认可度。脚本围绕读取、预处理、转换计算与结果输出组织流程可据此理解经验公式中地理位置、季节与大气透明度等因素的接入方式并借助pandas批量处理缺失值、异常值与日期格式问题还可扩展matplotlib绘图与并行加速。读者既能获得可复用的换算模板也能通过输入输出样例快速核对计算结果适合作为气象数据批处理的入门与练手案例。1. 840 个站点的日照时数要变成辐射值先得算清分母某开发者手上有一份日照时数的日值数据八年的记录准备喂给一个蒸散模型结果卡在最后一步模型要的是日总太阳辐射单位 MJ/m²数据里只有每天日照了几个小时。他先试了乘一个固定系数跟辐射观测比对夏天差了将近三成冬天又反过来偏低。这不是参数没调好是方法从根上就不对日照时数只反映晴天占了多少同样的 8 小时高纬度夏天只占可照时数的五成低纬度却接近七成落到地面的日总太阳辐射压根不是一回事。这份 Python 脚本处理的就是这一步读入 840 个气象站点的站点表与日照时数日值文件逐站逐日先算天文辐射和可照时数再用 Ångström–Prescott 公式换算成日总太阳辐射输出带单位标注的 CSV。适合手上只有日照时数、又急需辐射驱动场的场景比如蒸散计算、作物模型、光伏潜力粗估。难点不在代码量在于几个中间量的单位和边界条件这篇就按「公式 → 数据 → 代码 → 排错 → 改造」的顺序拆。2. 公式选型Ångström–Prescott、天文辐射与 a、b 系数的取值2.1 固定系数法错在哪分母随纬度和季节变把一个站点的日照时数直接乘 0.5 或者查一张「月份-系数」表是最省事也最容易出错的做法。问题在于日照时数本身不是能量它只是「太阳实际露脸多久」的时长统计要变成能量必须知道这段时间里大气层顶本来能接收多少辐射。这个「本来能接收多少」就是天文辐射 Ra它只跟日地距离、太阳赤纬和站点纬度有关跟云、气溶胶、海拔统统无关。所以换算链路是两步先把 Ra 和可照时数 N 算出来再用实测日照时数 n 除以 N 得到一个 0 到 1 之间的晴空比例最后乘上经验系数。Ångström–Prescott 形式最简单Rs Ra × (a b × n/N)其中 Rs 是日总太阳辐射Ra 是天文辐射n 是实测日照时数N 是天文可照时数a、b 是站点或区域的经验系数。这个式子的物理含义很直白a 代表全阴天时仍能到达地面的那部分散射辐射b 代表晴天比例每增加一个单位辐射增加多少。全阴天 n0 时 Rs a×Ra全晴天 nN 时 Rs (ab)×Ra两个端点都能对上物理直觉比拍脑袋乘系数靠谱得多。2.2 天文辐射与可照时数三个中间量必须先算对Ra 的标准算法FAO-56 里那套分三步。第一步算日地相对距离倒数 dr因为地球绕日轨道是椭圆一年里日地距离变化约 ±3.3%冬天近、夏天远对高纬度站点影响不能忽略。第二步算太阳赤纬 δ也就是太阳直射点纬度它决定了当天正午太阳高度角。第三步算日落时角 ωs由纬度与赤纬共同决定可照时数 N 就是 ωs 的两倍再换算成小时。三个量的表达式如下doy 是日序数1 到 365 或 366dr 1 0.033 × cos(2π × doy / 365) δ 0.409 × sin(2π × doy / 365 − 1.39) ωs arccos(−tan φ × tan δ) Ra (24×60/π) × Gsc × dr × [ωs × sinφ × sinδ cosφ × cosδ × sinωs] N 24/π × ωs这里 Gsc 是太阳常数取 0.0820 MJ m⁻² min⁻¹前面那个 24×60/π 是把分钟辐射累积成全天并把角度换成弧度的组合系数。φ 是站点纬度北纬为正、南纬为负必须换成弧度再进三角函数。这套公式在 FAO-56 灌溉手册里有完整推导农业气象和蒸散计算领域基本都认它。符号含义单位典型量级Gsc太阳常数MJ m⁻² min⁻¹0.0820dr日地相对距离倒数无量纲0.97 ~ 1.03δ太阳赤纬rad−0.41 ~ 0.41φ站点纬度弧度rad取决于站点ωs日落时角rad0 ~ πRa天文辐射MJ m⁻² d⁻¹5 ~ 42N可照时数h8 ~ 16中低纬注意 ωs 那一行arccos 的入参是 −tanφ×tanδ在中高纬度夏季可能小于 −1冬季可能大于 1这就是极昼极夜的数学表现。直接扔给 arccos 会返回 NaN批量跑 840 个站点时只要有那么几个高纬站整张输出表就会莫名其妙出现空值这个坑在第 5 章还会展开。2.3 a、b 系数怎么定先用默认值再谈标定a、b 不是普适常数它们随气候区、海拔、气溶胶状况变化。国内多数台站拟合出来的 a 落在 0.18 到 0.25 之间b 落在 0.50 到 0.60 之间比较常用的缺省组合是 a0.20、b0.55。如果你的研究区没有现成的标定值用这组缺省值先跑通、看量级对不对比纠结系数更有意义。真要标定思路是拿同一站点、同一时段的实测辐射 Rs_obs 和算出来的 Ra、n/N 做线性回归以 n/N 为自变量、Rs_obs/Ra 为因变量截距就是 a斜率就是 b。至少要有三年逐日数据回归的 R² 一般能到 0.85 以上达不到就说明该站的日照时数记录有问题或者辐射观测本身有系统性偏差。分月标定是进阶做法因为雨季和旱季的云型不同全年一组系数在季风区会明显低估雨季、高估旱季。2.4 先拿两个特殊日验算确认公式没抄错写完 Ra 函数不要急着跑 840 个站。挑一个纬度 40°N 的站点、日序数 172夏至前后手算或用独立脚本验一遍量级对不对一眼就能看出来。import numpy as np # 只验证单点纬度 40 度日序数 172夏至前后 lat, doy 40.0, 172 phi np.radians(lat) dr 1 0.033 * np.cos(2 * np.pi * doy / 365) delta 0.409 * np.sin(2 * np.pi * doy / 365 - 1.39) ws np.arccos(-np.tan(phi) * np.tan(delta)) ra (24 * 60 / np.pi) * 0.0820 * dr * ( ws * np.sin(phi) * np.sin(delta) np.cos(phi) * np.cos(delta) * np.sin(ws) ) print(fRa{ra:.2f} MJ/m2/d, N{24/np.pi*ws:.2f} h) # 期望输出Ra≈41.8 MJ/m2/d, N≈14.8 h参数说明lat是纬度度数doy是日序数dr用 365 做分母是 FAO 的标准写法闰年换 366 只影响小数第三位可以忽略。这段代码的价值不在计算在于给你一个可对照的锚点40°N 夏至 Ra 应该在 41 到 42 MJ m⁻² d⁻¹ 之间可照时数 14.8 小时左右。如果跑出来是 4.18八成是 Gsc 少乘了一个 10如果是 418那是把 24×60/π 重复乘了。再从反面验一个同一纬度、日序数 355冬至前后Ra 应该在 12 到 13 MJ m⁻² d⁻¹N 约 9.2 小时。夏冬两个端点的比值大约 3.4 倍这个比例关系比绝对数值更能说明公式写对了没有。3. 输入数据与目录结构站点表、日照日值文件和输出约定3.1 站点表三个字段必须能对上840 个站点的元信息至少要包含站号、纬度、经度海拔可选但建议留着后面做气压订正或者质量筛查用得上。站号一定要按字符串读不要让它自动变成整数——有些站号带前导零pd.read_csv默认会把它转成 int一读进来就少了一位跟日照文件里的站号对不上最后的结果是全部站点都被跳过而日志里只打印一行「缺文件」排查起来非常费劲。站点表的列名建议统一成station, lat, lon, elev。如果原始表里是「经度」「纬度」「台站名」这种中文列名在读入后立刻改掉别让中文列名流进后面的函数itertuples取属性时中文名会变成位置索引代码可读性直接崩掉。海岸线附近的站点要特别留意符号西经、南纬必须转成负数国内站点虽然只涉及东经北纬但如果以后要把代码迁到其他区域符号约定最好一开始就定死。还有一类常见情况是站点表里经纬度用度分秒表示比如116°28′这种字符串。这种列 pandas 会读成 object直接进np.radians就报错。写一个转换函数用正则把度分秒拆出来再合并成十进制度转换完打印前五行人工核对一遍比后面查错省事。3.2 日照时数日值文件格式和缺测码先摸清楚日值数据集常见的组织方式是每站一个文本文件或者全国一个文件按站号排序。文本里每一行是一条日记录字段大致是站号、年、月、日、日照时数、质量控制码。分隔符可能是空格、制表符也可能是定宽格式。这三种情况在 pandas 里对应三种读法选错了不会报错只会读出一列全是 NaN 或者一堆错位的字段。空格或制表符分隔pd.read_csv(f, sepr\s, headerNone)定宽格式pd.read_fwf(f, widths[...], headerNone)widths 要按实际字符位置数出来有表头且字段名规整直接read_csv加dtype{station: str}日照时数的缺测码是最容易出问题的地方。日值数据集里常用 32766 表示缺测也有用 999999 或者 −9999 的不同批次、不同要素的约定可能不一样。处理原则很简单先df[ssd].describe()看一眼最大值如果出现 32766 这种明显不是小时数的值就是缺测码确认后统一替换成pd.NA再参与计算。千万别用「大于 24 就截断成 24」这种想当然的处理那会把缺测硬生生变成「全天满日照」辐射值直接顶到天文辐射上限。质量控制码也别急着丢。质控码为 0 或 1 通常表示数据正常或可疑但可用2 以上往往意味着该条记录有问题。初跑阶段可以保留这一列后面统计一下各质控码的占比如果某站的异常码比例超过 10%这个站的输出单独标记出来别跟正常站混在一张表里。3.3 目录布局与命名约定840 个站点的数据量不算大单站八年逐日也就三千行左右合计两百多万行一台普通笔记本十几秒能跑完。但目录结构如果不提前定好跑到一半缺一个文件、多一个重复站号重跑的成本比想象中高。我习惯按下面这样摆project/ ├── data/ │ ├── stations.csv # 站号,lat,lon,elev │ └── ssd/ # 每站一个 txt文件名即站号 │ ├── 5xxxx.txt │ └── ... ├── core/ │ └── radiation.py # 纯计算函数不碰 IO ├── run_all.py # 批量调度 └── output/ └── rs_daily.csvcore/radiation.py里只放数学函数输入是数组、输出是数组不读文件也不写文件。这样做的好处是单测好写改系数的时候不会牵一发动全身。run_all.py负责遍历站点、读文件、调函数、拼结果。文件名直接用站号省掉在站点表里再维护一个文件路径列也避免了路径里带空格或中文导致的编码问题。3.4 先跑一遍数据体检再进计算正式计算之前写一个十几行的体检脚本把三类问题提前暴露站点表有没有重复站号、日照文件有没有缺失、日照时数的分布有没有异常。这一步花五分钟能省掉后面几个小时的翻车排查。from pathlib import Path import pandas as pd sites pd.read_csv(data/stations.csv, dtype{station: str}) print(站点数:, len(sites), 重复站号:, sites[station].duplicated().sum()) print(纬度范围:, sites[lat].min(), sites[lat].max()) ssd_dir Path(data/ssd) missing [s for s in sites[station] if not (ssd_dir / f{s}.txt).exists()] print(缺文件站点数:, len(missing), missing[:10]) # 抽查一个站看日照时数的分布 sample pd.read_csv(ssd_dir / f{sites[station].iloc[0]}.txt, sepr\s, headerNone, names[station, year, month, day, ssd, qc]) print(sample[ssd].describe()) print(疑似缺测码出现次数:, (sample[ssd] 32766).sum())逻辑上分三段第一段核对站点表自身的一致性重复站号和纬度越界是最常见的两类脏数据第二段用集合差找出缺文件的站注意这里没做异常捕获文件不存在就是不存在早暴露早处理第三段抽一个站看ssd的统计量max如果是 32766 就说明缺测码没处理min如果是负数说明有异常记录。参数上dtype{station: str}是必须的不写这一句站号读成整数后跟文件名比对会全部失败。names里的字段顺序要跟原始文件严格对齐字段数不匹配时 pandas 会报ParserError报错信息里会给出实际列数照着改就行。这一步跑完如果缺文件站点数是 0、ssd的 max 在 16 以内、min 不小于 0就可以放心进核心计算了。4. 核心模块拆解三个函数走完从小时数到 MJ/m² 的链路4.1 函数一extraterrestrial_radiation 算 Ra 和 N这个函数是整个链路的地基输入日序数和纬度输出天文辐射与可照时数两个数组。写成向量化形式一次能算完一个站八年的所有日比逐行循环快两个数量级。import numpy as np GSC 0.0820 # 太阳常数, MJ m-2 min-1 MJ_TO_CAL 23.89 # 1 MJ m-2 23.89 cal cm-2 def extraterrestrial_radiation(doy, lat_deg): 逐日天文辐射 Ra(MJ m-2 d-1) 与可照时数 N(h)。 doy : 日序数标量或数组1~366 lat_deg : 纬度北正南负单位度 doy np.asarray(doy, dtypefloat) phi np.radians(np.asarray(lat_deg, dtypefloat)) dr 1.0 0.033 * np.cos(2.0 * np.pi * doy / 365.0) delta 0.409 * np.sin(2.0 * np.pi * doy / 365.0 - 1.39) x -np.tan(phi) * np.tan(delta) x np.clip(x, -1.0, 1.0) # 极昼极夜保护避免 arccos 出 NaN ws np.arccos(x) ra (24.0 * 60.0 / np.pi) * GSC * dr * ( ws * np.sin(phi) * np.sin(delta) np.cos(phi) * np.cos(delta) * np.sin(ws) ) n_day 24.0 / np.pi * ws return ra, n_day逐段说明。np.asarray把标量和列表统一成数组这样后面全是逐元素运算调用方传一个数或者传一整列都不会出错。dr和delta两行是纯周期函数分母用 365 是 FAO 的标准写法闰年用 366 会让结果差千分之几对辐射量级没有实质影响不必为此加分支。np.clip(x, -1.0, 1.0)这一句是整段代码里最容易被忽略、也最容易出事的地方。纬度超过 66.5° 的站点在夏季会出现 x −1冬季出现 x 1此时数学上对应极昼和极夜物理上 ωs 应该取 π 和 0。不截断的话 arccos 返回 NaNNaN 会沿着乘法一路污染整行结果最后的 CSV 里就出现莫名其妙的空值而你还以为只是某个站的数据缺了几天。返回的两个数组单位分别是 MJ m⁻² d⁻¹ 和小时后面所有计算都基于这两个量单位不要在中途改来改去。4.2 函数二sunshine_to_radiation 做换算和裁剪拿到 Ra 和 N第二步就是套 Ångström–Prescott。这个函数很短但两个边界处理不能省。def sunshine_to_radiation(n_hours, ra, n_day, a0.20, b0.55): Ångström-Prescott 换算Rs Ra * (a b * n/N)。 n_hours : 实测日照时数(h)缺测应已在外部剔除 ra : 天文辐射(MJ m-2 d-1) n_day : 可照时数(h) a, b : 经验系数缺省 0.20 / 0.55 n np.asarray(n_hours, dtypefloat) ratio np.divide(n, n_day, outnp.zeros_like(n, dtypefloat), wheren_day 0) ratio np.clip(ratio, 0.0, 1.0) # n 超过 N 或为负归到边界 return ra * (a b * ratio)np.divide带where参数的写法是为了处理 n_day 为 0 的情况——极夜时 N0直接除会得到 inf 或者警告。out先填 0where只对 N0 的位置做除法剩下的保持 0逻辑上就是「极夜时晴空比例为 0」虽然极夜里 Ra 本身也接近 0结果无论如何都是 0但至少不会污染下游。np.clip(ratio, 0, 1)处理的是观测记录问题。日照时数记录仪换型、人工观测误记都可能出现 n 略大于 N 的情况比如 N14.8 而记录写了 15.1。这种记录不裁掉的话a b×1.02会算出比天文辐射还大的地面辐射物理上不可能。裁到 1 相当于把它当成全天满日照是误差最小的处理方式。如果这种记录占比很高比如超过 3%就要回头查数据源了光在代码里裁是不够的。系数 a、b 作为默认参数放在签名里调用时可以逐站传入标定值不用改函数体。这也是把计算逻辑单独抽成一个模块的好处。4.3 函数三process_station 负责单站装配和落盘字段前两个函数是纯数学第三个函数负责把一年的数据装配起来生成最终要写进 CSV 的字段。这里要一次性把单位都算好别等到后面再补。import pandas as pd def process_station(meta, df, a0.20, b0.55): 单站装配输入站点元信息和日值表返回带辐射列的表。 df df.sort_values([year, month, day]).reset_index(dropTrue) dates pd.to_datetime(df[[year, month, day]]) doy dates.dt.dayofyear.to_numpy() ra, n_day extraterrestrial_radiation(doy, meta[lat]) rs sunshine_to_radiation(df[ssd].to_numpy(), ra, n_day, a, b) out df[[station, year, month, day]].copy() out[doy] doy out[ra_mj] np.round(ra, 3) out[n_day_h] np.round(n_day, 2) out[ssd_h] df[ssd].to_numpy() out[rs_mj] np.round(rs, 3) out[rs_cal] np.round(rs * MJ_TO_CAL, 1) # cal cm-2 d-1 out[rs_wm2] np.round(rs / 0.0864, 1) # 日均辐射通量 return outpd.to_datetime那两行是为了拿到日序数。有的数据集只给年月日不给日序数用 datetime 转是最稳的办法闰年 2 月 29 日会自动算成 60年底 12 月 31 日是 366不用自己写月份天数表。自己手写[0,31,59,...]累加数组碰上闰年就整体错一天这种错误在小样本上看不出来跑完八年才发现整条序列偏移非常难查。输出字段里刻意保留了三个单位的辐射rs_mj是主流模型要的 MJ m⁻² d⁻¹rs_cal是老一辈气象资料和部分农业手册用的 cal cm⁻² d⁻¹rs_wm2是日均辐射通量跟再分析产品的单位对得上。多算两列的成本几乎为零但能省掉后面反复换算、反复怀疑单位的麻烦。np.round只是为了让 CSV 好看如果下游还要做精确计算建议去掉 round 直接用原值。meta[lat]是单站的纬度df是已经去掉缺测、按日期排序的日值表。函数内部没有做任何 IO也不改传入的 df返回的是新表这样批量循环里出问题时容易定位到具体哪个站。4.4 批量调度840 个站点怎么稳定跑完批量脚本要做的事很机械读站点表、遍历、读文件、剔缺测、调函数、拼结果。真正需要注意的是异常处理和日志跑两百多万行的时候中间某一行报错不能把整个任务带崩。from pathlib import Path import pandas as pd from core.radiation import process_station ssd_dir Path(data/ssd) sites pd.read_csv(data/stations.csv, dtype{station: str}) frames, skipped [], [] for row in sites.itertuples(indexFalse): f ssd_dir / f{row.station}.txt if not f.exists(): skipped.append((row.station, 缺文件)) continue raw pd.read_csv(f, sepr\s, headerNone, names[station, year, month, day, ssd, qc]) raw[ssd] pd.to_numeric(raw[ssd], errorscoerce) raw.loc[raw[ssd] 32766, ssd] pd.NA # 缺测码剔除 raw raw.dropna(subset[ssd]) if raw.empty: skipped.append((row.station, 全为缺测)) continue frames.append(process_station({station: row.station, lat: row.lat}, raw)) result pd.concat(frames, ignore_indexTrue) result.to_csv(output/rs_daily.csv, indexFalse, float_format%.3f) print(完成站点:, len(frames), 跳过:, skipped[:5])itertuples比iterrows快而且拿到的row.station、row.lat就是列名属性比row[station]好读。缺文件不抛异常记到skipped里继续下一个这样一趟跑完你能拿到完整的跳过清单而不是跑到第 500 个站崩掉重来。pd.to_numeric(..., errorscoerce)是第二道保险万一某个文件里混进了非数字字符这一句会把它变成 NaN 而不是让整列变成 object。raw.loc[raw[ssd] 32766, ssd] pd.NA处理缺测码如果你的数据源用的是 999999改这个阈值就行。最后float_format%.3f控制输出精度避免 CSV 里出现一长串浮点尾数。840 个站、八年逐日总行数大约两百多万pd.concat一次性拼完内存压力不大普通笔记本 8G 内存够用。如果站数翻到几千或者年份拉长到三十年就得改成分批写文件再合并或者直接上polars。不过那是另一个话题了当前规模没必要提前优化。5. 常见问题排查批量跑 840 站时最容易翻车的五个点5.1 结果整体偏大或偏小两成现象所有站点的辐射值跟观测比对系统性偏高或偏低 15% 到 25%各站方向一致不是随机误差。原因大概率是单位混用。常见的有三种一是 Ra 用了 MJ 但输出标成 cal/cm² 没换算差了 23.89 倍二是把日均辐射通量 W/m² 直接当成了日累积量 MJ/m²/d两者差 0.0864 倍三是系数 a、b 用的是「a0.25、b0.50」这类老文献值但公式写成了a b×n/N之外的形式比如漏了 b 只乘了 n。解决办法是先算一个全晴天样本对照取某站夏季日照时数接近 N 的几天Rs 应该略小于 Ra比值在 0.70 到 0.75 之间。如果算出来比值是 0.03 或者 17就是单位问题如果比值在 0.5 左右那是系数问题。这个对照法比逐项查代码快得多我一般先做这个再动代码。5.2 高纬度站点算出 NaN现象输出表里部分站点的rs_mj是空值集中在纬度偏高的站而且往往集中在夏季或冬季的某几天。原因就是 2.2 节提到的arccos越界。纬度接近或超过极圈时−tanφ×tanδ 的绝对值会大于 1numpy 返回 NaN并且给出 RuntimeWarning。如果日志被刷屏这个警告很容易被忽略。解决就是在 arccos 之前加np.clip(x, -1.0, 1.0)。截断后极昼时 ωsπ、N24 小时极夜时 ωs0、N0物理上是对的。如果站点确实在极圈内还要考虑极昼期间日照时数记录本身是否有效这时候更稳妥的做法是对 N24 的日子单独打标记而不是直接输出。5.3 日照时数大于可照时数现象算出来的 Rs 比 Ra 还大或者n/N大于 1某些日子辐射值明显不合理。原因有两类。一类是观测记录问题仪器换型、人工读数误记日照时数记成了 15 小时而当天 N 只有 14.2。另一类是公式里的 N 算错了比如纬度符号搞反北纬传成了负值算出来的可照时数跟实际对不上。解决办法分两步先在sunshine_to_radiation里用np.clip(ratio, 0, 1)兜底保证输出不会超过 Ra然后统计n N的记录占比超过 3% 就说明数据源有问题要回去查那个站的原始记录。单纯靠 clip 掩盖过去等于把数据质量问题藏起来了后面做趋势分析会出偏差。5.4 缺测码被当成真实值参与计算现象某几天的辐射值巨大可能是正常值的十几倍集中在特定站点、特定时段。原因就是 32766 这类缺测码没剔除直接进了乘法。32766 小时除以 N 得到几千的比值再乘 Ra结果能到几万 MJ/m²一眼就能看出来不对但如果只做汇总统计没逐行看很容易被平均值掩盖。解决是在读取阶段就把缺测码替换成 NA 并 drop 掉不要留到计算阶段再处理。同时建议在输出里加一列有效标记记录该行是实测还是缺测插补后面做时间序列分析时能区分开。缺测率高的站点比如超过 20%建议整站不进主表单独输出到另一个文件。5.5 站点表经纬度是度分秒字符串现象读站点表时lat列是 object 类型进np.radians直接抛 TypeError或者被 pandas 静默转成 NaN 后算出全站一样的辐射值。原因是原始站点表用「116°28′」这类格式记录经纬度read_csv无法识别成数字。这种情况在整理过的数据集里很常见尤其是从网页或者 PDF 转录过来的表。解决办法是写一个解析函数用正则(\d)[°度](\d)[′分]?拆出度和分再合并成度 分/60。转换完打印前五行跟原始表对照确认没有把纬度和经度列搞反。经纬度搞反是最隐蔽的错代码不报错结果看起来也「像那么回事」但整个中国的站点辐射分布会完全错位。我一般会在体检脚本里加一句检查国内站点的经度应该在 70 到 140 之间纬度在 15 到 55 之间超出范围就报警。6. 二次开发把系数、单位、时间分辨率换成你自己项目要的跑通只是起点真正用起来总要改点什么。这个项目里最常改的三个地方是经验系数、输出单位和时间尺度每个都有对应的改法不用动核心公式。改造点改哪里怎么改验证方式分区标定 a、bprocess_station的入参站点表里加a_coef、b_coef两列循环时传进去与实测辐射回归看 R² 和偏差分月系数sunshine_to_radiation内部系数改成 12 元素数组按month索引逐月偏差是否收敛输出 cal/cm²输出字段保留rs_cal或改MJ_TO_CAL精度与历史气象资料对照月尺度汇总run_all.py之后按站号年月 groupby 求均值注意是日均值不是月总量月总量除以天数应接近日均换数据源读取段改分隔符、列名、缺测码阈值体检脚本先跑一遍分月系数是我改得最多的一项。季风气候区全年一组系数雨季低估、旱季高估误差能到 15%。把 a、b 换成 12 个值的数组函数里用月份做索引改动不到十行def sunshine_to_radiation_monthly(n_hours, ra, n_day, month, a_arr, b_arr): 分月系数的 Ångström-Prescotta_arr/b_arr 为长度 12 的数组。 idx np.asarray(month, dtypeint) - 1 # 月份 1~12 映射到 0~11 a np.asarray(a_arr, dtypefloat)[idx] b np.asarray(b_arr, dtypefloat)[idx] ratio np.clip(np.divide(n_hours, n_day, outnp.zeros_like(n_hours, dtypefloat), wheren_day 0), 0.0, 1.0) return ra * (a b * ratio)这里的idx一定要减 1月份从 1 开始而数组从 0 开始忘了减就会取到下一个月的系数冬天用成春天的值误差不大不小正好让你查不出来。系数数组的长度也要校验传进来 11 个值的时候索引 11 会越界加一句assert len(a_arr) 12能省不少事。另一个值得做的改造是加一层结果校验。算完不要直接拿去用挑三到五个有辐射观测的站点把算出来的 Rs 和实测值放一起看散点和偏差。我习惯算三个指标平均偏差MBE、均方根误差RMSE、相关系数。日照时数换算的辐射日尺度 RMSE 一般在 2 到 3 MJ m⁻² d⁻¹相关系数 0.9 以上算合格。如果 RMSE 超过 4先查缺测处理再查系数最后才怀疑公式本身。血泪经验是别在同一个脚本里既读数据又算公式又画图改一个系数要重跑整条链路等得人心焦。计算模块和可视化模块彻底分开中间用 CSV 交换改系数只需要重跑计算那一步几分钟出结果。从那以后我每次接这类换算任务都强制先把「纯函数 独立批量脚本 体检脚本」三件套搭好再动手写公式希望帮到你。本文还有配套的精品资源点击获取
