简介：面向气象、环境科学与太阳能资源评估方向的研究者及数据分析爱好者，提供一套将全国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 代表晴天比例每增加一个单位，辐射增加多少。全阴天 n=0 时 Rs = a×Ra，全晴天 n=N 时 Rs = (a+b)×Ra，两个端点都能对上物理直觉，比拍脑袋乘系数靠谱得多。

2.2 天文辐射与可照时数：三个中间量必须先算对

Ra 的标准算法（FAO-56 里那套）分三步。第一步算日地相对距离倒数 dr，因为地球绕日轨道是椭圆，一年里日地距离变化约 ±3.3%，冬天近、夏天远，对高纬度站点影响不能忽略。第二步算太阳赤纬 δ，也就是太阳直射点纬度，它决定了当天正午太阳高度角。第三步算日落时角 ωs，由纬度与赤纬共同决定，可照时数 N 就是 ωs 的两倍再换算成小时。

三个量的表达式如下，doy 是日序数（1 到 365 或 366）：

dr = 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 灌溉手册里有完整推导，农业气象和蒸散计算领域基本都认它。

注意 ω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 之间，比较常用的缺省组合是 a=0.20、b=0.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(f"Ra={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, sep=r"\s+", header=None)
• 定宽格式：pd.read_fwf(f, widths=[...], header=None)，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.csv

core/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", sep=r"\s+", header=None, 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, dtype=float) phi = np.radians(np.asarray(lat_deg, dtype=float)) 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 返回 NaN，NaN 会沿着乘法一路污染整行结果，最后的 CSV 里就出现莫名其妙的空值，而你还以为只是某个站的数据缺了几天。

返回的两个数组单位分别是 MJ m⁻² d⁻¹ 和小时，后面所有计算都基于这两个量，单位不要在中途改来改去。

4.2 函数二：sunshine_to_radiation 做换算和裁剪

拿到 Ra 和 N，第二步就是套 Ångström–Prescott。这个函数很短，但两个边界处理不能省。

def sunshine_to_radiation(n_hours, ra, n_day, a=0.20, b=0.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, dtype=float) ratio = np.divide(n, n_day, out=np.zeros_like(n, dtype=float), where=n_day > 0) ratio = np.clip(ratio, 0.0, 1.0) # n 超过 N 或为负，归到边界 return ra * (a + b * ratio)

np.divide带where参数的写法是为了处理 n_day 为 0 的情况——极夜时 N=0，直接除会得到 inf 或者警告。out先填 0，where只对 N>0 的位置做除法，剩下的保持 0，逻辑上就是「极夜时晴空比例为 0」，虽然极夜里 Ra 本身也接近 0，结果无论如何都是 0，但至少不会污染下游。

np.clip(ratio, 0, 1)处理的是观测记录问题。日照时数记录仪换型、人工观测误记，都可能出现 n 略大于 N 的情况，比如 N=14.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, a=0.20, b=0.55): """单站装配：输入站点元信息和日值表，返回带辐射列的表。""" df = df.sort_values(["year", "month", "day"]).reset_index(drop=True) 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 out

pd.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(index=False): f = ssd_dir / f"{row.station}.txt" if not f.exists(): skipped.append((row.station, "缺文件")) continue raw = pd.read_csv(f, sep=r"\s+", header=None, names=["station", "year", "month", "day", "ssd", "qc"]) raw["ssd"] = pd.to_numeric(raw["ssd"], errors="coerce") 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_index=True) result.to_csv("output/rs_daily.csv", index=False, float_format="%.3f") print("完成站点:", len(frames), "跳过:", skipped[:5])

itertuples比iterrows快，而且拿到的row.station、row.lat就是列名属性，比row["station"]好读。缺文件不抛异常，记到skipped里继续下一个，这样一趟跑完你能拿到完整的跳过清单，而不是跑到第 500 个站崩掉重来。

pd.to_numeric(..., errors="coerce")是第二道保险，万一某个文件里混进了非数字字符，这一句会把它变成 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 用的是「a=0.25、b=0.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δ 的绝对值会大于 1，numpy 返回 NaN，并且给出 RuntimeWarning。如果日志被刷屏，这个警告很容易被忽略。

解决就是在 arccos 之前加np.clip(x, -1.0, 1.0)。截断后极昼时 ωs=π、N=24 小时，极夜时 ωs=0、N=0，物理上是对的。如果站点确实在极圈内，还要考虑极昼期间日照时数记录本身是否有效，这时候更稳妥的做法是对 N=24 的日子单独打标记，而不是直接输出。

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. 二次开发：把系数、单位、时间分辨率换成你自己项目要的

跑通只是起点，真正用起来总要改点什么。这个项目里最常改的三个地方是经验系数、输出单位和时间尺度，每个都有对应的改法，不用动核心公式。

分月系数是我改得最多的一项。季风气候区全年一组系数，雨季低估、旱季高估，误差能到 15%。把 a、b 换成 12 个值的数组，函数里用月份做索引，改动不到十行：

def sunshine_to_radiation_monthly(n_hours, ra, n_day, month, a_arr, b_arr): """分月系数的 Ångström-Prescott，a_arr/b_arr 为长度 12 的数组。""" idx = np.asarray(month, dtype=int) - 1 # 月份 1~12 映射到 0~11 a = np.asarray(a_arr, dtype=float)[idx] b = np.asarray(b_arr, dtype=float)[idx] ratio = np.clip(np.divide(n_hours, n_day, out=np.zeros_like(n_hours, dtype=float), where=n_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 交换，改系数只需要重跑计算那一步，几分钟出结果。从那以后我每次接这类换算任务，都强制先把「纯函数 + 独立批量脚本 + 体检脚本」三件套搭好再动手写公式，希望帮到你。

本文还有配套的精品资源，点击获取