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 与大气顶层辐射 H0

2.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 代表日照增加时辐射的斜率。a+b 通常在 0.7 到 0.8 之间，超过 0.9 就要怀疑系数来源或者数据单位。公式本身不复杂，难的是逐日逐站把 N 和 H0 算对，并且把缺测、单位、merge 这些脏活处理干净。

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 rad，N 约 14.35 小时，dr 约 0.968，H0 约 41.6 MJ/m2/day。如果日照时数 n=8.2 小时，n/N=0.571，取 a=0.25、b=0.50，H=41.6*(0.25+0.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.32，b 保持在 0.45 至 0.50；偏湿润站点云多，a 可以降到 0.18 至 0.22，b 略高到 0.52 至 0.58。这些范围只是筛查起点，不是定论。

有辐射观测的站点，直接把H_obs/H0对n/N做线性回归，截距是 a，斜率是 b。样本至少覆盖一年四季，否则冬季样本少会把系数带偏。回归后检查 a 是否在 0.1 至 0.4，b 是否在 0.3 至 0.7，a+b 是否在 0.6 至 0.9。超出范围的，先回头查单位，不要急着调系数。

3. Python 脚本落地：长表合并、向量化计算与 840 站点批处理

3.1 输入表结构与字段约定

推荐把站点表和日照时数表分开。站点表一行一个站，日照表一行一个站点一天。字段名可以按实际数据改，但脚本里要固定映射，不要一边跑一边改列名。站号统一按字符串读，避免001被读成1。日期统一转datetime64，缺测码统一替换成NaN，不要让-999进入公式。

长表的好处是 840 个站点可以一次 merge 完，不需要循环读 840 个文件。数据量按 840 站 × 365 天 × 30 年估算，大约 900 万行，pandas 分块读完全能处理。关键是 merge 前把站点表去重，merge 时加validate="many_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], errors="coerce") 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. 日地距离修正与大气顶层辐射 H0，MJ/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(lower=0, upper=1) 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。gsc=0.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_ratio=False只在排查异常时用，输出保留原始 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", required=True) parser.add_argument("--ssh", required=True) parser.add_argument("--out", required=True) parser.add_argument("--a", type=float, default=0.25) parser.add_argument("--b", type=float, default=0.50) parser.add_argument("--chunksize", type=int, default=500_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}, chunksize=args.chunksize): chunk["date"] = pd.to_datetime(chunk["date"], errors="coerce") chunk = chunk.dropna(subset=["date"]) chunk["ssd"] = pd.to_numeric(chunk["ssd"], errors="coerce") chunk.loc[chunk["ssd"] <= -900, "ssd"] = pd.NA # 常见缺测码 n_before = len(chunk) chunk = chunk.merge( sites, on="site_id", how="left", validate="many_to_one" ) assert len(chunk) == n_before, "merge 后行数膨胀，检查站点表重复" chunk = add_radiation(chunk, a=args.a, b=args.b) pieces.append( chunk[["site_id", "date", "ssd", "N", "H0", "n_N", "H"]] ) result = pd.concat(pieces, ignore_index=True) Path(args.out).parent.mkdir(parents=True, exist_ok=True) result.to_csv(args.out, index=False, float_format="%.4f") print(f"rows={len(result)}, out={args.out}") if __name__ == "__main__": main()

chunksize控制单次读入行数，50 万行对普通笔记本比较稳。validate="many_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，即截距在前、斜率在后。标定完把a+b打出来，正常在 0.6 到 0.9；如果a是负数或者b大于 1，先查观测辐射单位是不是 W/m2 没转成 MJ/m2/day。

4.2 无观测时的系数迁移与合理性检查

没有辐射观测的站点，用相近气候区的标定系数，或者先用默认值。迁移系数时不要跨气候带硬套，至少按干旱/湿润、低纬/高纬分两组。跑完做四类检查：n_N是否在 0 到 1，H0是否在 0 到 45 MJ/m2/day，H/H0是否在 0.05 到 0.9，H是否为非负。任何一条不满足，先标记，不要直接删。

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(axis=1) return out

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.0，H甚至超过H0。原因通常不是公式错，而是日照时数单位混了：有的表按分钟存，有的按小时存，还有的用 0.1 小时为单位。另一个原因是日长 N 算错，比如纬度用了弧度又转了一次。解决：先在数据字典里确认ssd单位，分钟就除以 60；进公式前把负值和缺测码置空；n_N截断到 0 到 1 只作为兜底，不要用它掩盖单位问题。跑完统计n_N > 1的站点数，超过 1% 就停下来查源表。

5.2 闰年日序用 365：2 月 29 日附近 H0 会跳

现象：闰年 2 月 29 日前后H0出现不连续，或者 12 月 31 日的赤纬和往年对不上。原因：赤纬公式里的 D 固定写成了 365，闰年应该用 366。解决：用d.dt.is_leap_year生成days_in_year，再代入赤纬和日地距离修正。不要用doy/365一把梭。这个错误很小，但在时间序列里会变成一个尖刺，后面做月均值时容易被当成异常值删掉。

5.3 高纬站点 arccos 出现 NaN：tanφ 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 时加validate="many_to_one"，merge 后断言行数等于输入行数。站号按字符串读，避免A001和a001被当成两个站。如果站点表里一个站有多行经纬度，先确认哪一行是有效坐标，不要直接全量 merge。

5.5 缺测码 -999 参与计算：H 出现负值或巨大值

现象：H出现 -200 或者 1e6 量级的值，n_N严重偏离。原因：缺测码-999、-9999没有替换，直接进了公式。解决：读入后统一replace([-999, -9999], np.nan)，或者用pd.to_numeric(errors="coerce")加阈值过滤。缺测行不要删，保留ssd为空、H为空，加flag_missing列，后面统计缺测率。缺测率超过 20% 的站点，标定系数时不要用。

6. 进阶用法：分块并行与日常任务化，别让 840 站跑成黑匣子

840 个站点如果每个站点 30 年日值，大约 900 万行，单进程 pandas 也能跑，但每次调参都全量重算会很拖。我一般先用--sample-sites 10抽 10 个站点跑通，再用分块并行。分块不要按行随机切，按站点切，保证一个站点的日期连续。参数建议如下：

按站点并行的骨架很短，每个进程拿到一个站点的子表，调同一个add_radiation，最后 concat。注意进程数不要超过物理核数，840 个站点分 8 组比开 840 个进程稳。

from multiprocessing import Pool def run_one(item): site_id, sub = item return add_radiation(sub, a=0.25, b=0.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_index=True)

日常任务化时，把a、b、chunksize、sample_sites写进配置文件，不要每次改脚本。输出加一列run_id，记录这次用的系数和日期，后面对比不同版本结果时能追。最省事的验证方法是抽 10 个站点，用手算的 N、H0 和脚本对一遍，误差在 0.1 以内再放全量。从那以后我每次全量前都强制走一遍抽样对账，确认H0在 0 到 45、n_N在 0 到 1、H非负，再让 840 个站点跑到底。希望帮到你。

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