简介面向气象、环境科学与太阳能资源评估场景的Python处理工具用于将全国840个气象站点观测的日照时数批量换算为日总太阳辐射。脚本覆盖数据读取、缺失与异常值预处理、经验公式转换及逐站点结果输出适合具备基础Python与pandas能力的科研人员和数据分析爱好者参考复用。压缩包为rar格式共3个文件含2个CSV数据文件与1个Python脚本整体约275KB便于解压后对照输入输出理解计算流程。已有6348人学习下载其在气象数据批处理与太阳能资源评估方面具有一定参考价值。读者可获得一套可运行、可修改的日照时数转日总太阳辐射实现范例并借以梳理站点数据组织、公式落地和结果校验思路也可扩展可视化、批量处理或接入自有站点数据。1. 只有日照时数没有辐射840 站点批量补出日总辐射的实战起点手里只有日照时数、没有辐射是气象数据二次加工里最常见的缺口。全国 840 个气象站点逐日观测里日照时数是常规项目太阳辐射却只有寥寥几十个辐射观测站才配得上仪器。做作物生长模型、参考蒸散量估算、光伏资源初筛的时候缺了这一列日总太阳辐射整条链路就卡在中间。这个项目解决的就是这段缺口用一份全国气象站点数据配一套 python 脚本把逐日日照时数按 Ångström-Prescott 关系换算成日总太阳辐射统一输出成 MJ/m²/d顺带给出 W/m² 和天文辐射上限供质量校验。适合已经拿到站点日照时数表、要批量补出辐射列的人也适合想把换算过程吃透、自己改系数的人。2. 技术选型与数据约定Ångström-Prescott 为什么比辐射传输模型更适合这张表2.1 从日照时数到日总辐射的物理链路先把链路讲清楚不然代码写出来也是照抄。日照时数 n 的定义是“太阳直射辐照度超过 120 W/m² 的累计时长”单位小时。它本身不是能量但相对日照 n/NN 是当天可照时长是一个相当靠得住的云量代理量晴空时 n 逼近 N比值趋近 1全阴天 n 接近 0比值也趋近 0。Ångström 最早给出的就是这条线性关系后来 Prescott 把上限从“晴天辐射”换成便于计算的天文辐射 Ra就有了现在通用的形式H / Ra a b * (n / N)这里的 a 可以理解成全天阴时的透射率b 是相对日照每增加一个单位所带来的额外透射ab 大致对应晴空透射率上限。默认值取 0.25 和 0.50 时ab0.75意味着晴空条件下地面能收到天文辐射的 75%这个量级在中国大部分地区是合理的。整条链路需要三样输入站点纬度算 Ra 和 N、日期算日序 J、赤纬 δ、日照时数 n。经度在纯 Ångström 公式里其实用不上但它在判定日界、做地方时修正的时候是必需的所以数据表里最好留着。注意日照时数必须和 Ra 落在同一天。跨时区站点如果按北京时统计日照而日序按世界时算西部站点会整体错位一天误差不大但很隐蔽往往要等到和实测辐射对比时才暴露。2.2 技术栈pandas 加 numpy算这张表够用了840 个站点乘上一年 365 天也就三十万行量级一张普通宽表而已。这个体量下选型的核心不是性能是可读性和可复现性。不用 xarray 和 NetCDF。xarray 的优势在多维格点、按维度切片、懒加载站点数据是典型的长表结构用 xarray 属于拿高射炮打蚊子反而让后来接手的人多学一套索引语法。除非你后面要把结果插值成格点再和再分析资料对比那时候再转。也不用 pvlib。pvlib 的强项是晴空辐照度模型、太阳位置算法、组件建模Ångström 这种经验回归反而不在它的主路径上。为了一个线性公式拉进来一整套依赖在服务器上装环境的时候容易翻车。更不需要 numba 或 Cython。整个计算是纯向量化的三角函数加四则运算numpy 广播一开三十万行几百毫秒就跑完了加编译层属于自我感动。依赖清单就三样pandas 负责读写和日期处理numpy 负责向量化计算pathlib 负责文件遍历。如果要做并行再加标准库的 concurrent.futures仍然是零外部依赖。2.3 输入输出字段约定数据格式不统一是这类批量脚本最常见的返工原因。我一般会在动手前先把字段钉死脚本里用 dtype 显式声明别指望 pandas 自动推断——它把站点号读成整数、把缺测码读成正常值的事我见过不止一次。输入表长这样字段类型单位说明station_idstring-站点号保留前导零必须按字符串读latfloat度纬度北纬为正中国境内取 15~55lonfloat度经度东经为正主要用于日界核对elevfloat米海拔Ångström 公式本身不用但做气压订正时会用到datedatetime-观测日期必须能解析成 YYYY-MM-DDsshfloat小时日照时数缺测用 NaN输出表在原始列基础上追加四列ra_mj天文辐射MJ/m²/d、daylen_h可照时长 N小时、rs_mj换算后的日总辐射MJ/m²/d、rs_wm2同一物理量的功率表达W/m²。两个单位都输出是因为下游模型的口味不一样——作物模型多用 MJ/m²/d光伏和建筑能耗模拟习惯看 W/m²。换算关系是 1 MJ/m²/d 10⁶/86400 W/m² ≈ 11.574 W/m²这个常数写进代码里比每次现算靠谱。3. 核心算法拆解从日序 J 一路推到日总辐射 Rs3.1 目录结构脚本按“天文量 / 经验公式 / 批处理”三层拆开每个文件职责单一改系数的时候不用翻遍全文solar_rad/ ├── config.py # 系数、路径、缺测码等全部可调参数 ├── astro.py # 日序、赤纬、日落时角、Ra、日长 N ├── angstrom.py # Ångström-Prescott 主公式与单位换算 ├── batch_run.py # 遍历原始目录、写出结果、生成质检报告 ├── data/ │ ├── raw/ # 按站点或按年切的原始 csv │ └── out/ # 结果 csv └── tests/ └── test_astro.py # 用已知日期对 Ra 做回归测试有人喜欢把所有东西塞进一个 200 行的 main.py前期确实快但等到要按分区换系数、要加并行、要写单元测试的时候那个文件会迅速变成黑匣子。三层拆分多花的十分钟后面能省回来。3.2 日序、赤纬与日落时角三个天文量是后续所有计算的底座。日序 J 从日期直接取赤纬 δ 和日落时角 ωs 用 FAO-56 推荐的近似式形式简单且在中国范围内精度足够# astro.py import numpy as np SOLAR_CONST 0.0820 # MJ m-2 min-1FAO-56 采用的太阳常数 def declination(doy): 赤纬 δ输入日序 J(1~366)返回弧度 return 0.409 * np.sin(2.0 * np.pi * doy / 365.0 - 1.39) def sunset_angle(doy, lat_deg): 日落时角 ωs单位弧度对极昼极夜做安全裁剪 phi np.deg2rad(lat_deg) delta declination(doy) # -tan(phi)*tan(delta) 超出 [-1,1] 时物理上就是极昼/极夜 x np.clip(-np.tan(phi) * np.tan(delta), -1.0, 1.0) return np.arccos(x)declination里的 0.409 是赤纬振幅弧度约 23.45°1.39 是相位修正项这两个常数来自对日地几何的拟合不需要每年改。sunset_angle里的np.clip是关键纬度高的站点在冬至前后-tan(φ)tan(δ)可能大于 1arccos 会直接返回 NaN一传十十传百整列结果全废。裁剪到 [-1,1] 之后arccos 返回 0 或 π对应日长 0 或 24 小时物理意义正好对上极夜和极昼。3.3 天文辐射 Ra 与可照时长 N有了 ωsRa 和 N 都是闭式解。Ra 的公式看着长其实就是“太阳常数 × 日地距离修正 × 日照半球上的投影积分”# astro.py续 def extraterrestrial_radiation(doy, lat_deg): 天文辐射 RaMJ m-2 d-1 phi np.deg2rad(lat_deg) delta declination(doy) ws sunset_angle(doy, lat_deg) # 日地距离修正因子 dr1 月初最大 dr 1.0 0.033 * np.cos(2.0 * np.pi * doy / 365.0) # 24*60/pi ≈ 458.37把 min-1 的太阳常数换成 d-1 ra (24.0 * 60.0 / np.pi) * SOLAR_CONST * dr * ( ws * np.sin(phi) * np.sin(delta) np.cos(phi) * np.cos(delta) * np.sin(ws) ) return ra def day_length(doy, lat_deg): 可照时长 N单位小时 ws sunset_angle(doy, lat_deg) return 24.0 / np.pi * wsdr这一项不能省。地球在 1 月初处于近日点日地距离比 7 月初近约 3.3%对应到 Ra 上就是北半球冬季晴天反而能收到更多天文辐射。省掉dr冬季结果会系统性偏低 3% 左右做长序列统计的时候这个偏差会一直存在。day_length的系数 24/π 来自 ωs 是半日角从日出到正午的角距乘 2 得全天角距再按 2π 对应 24 小时换算化简就是 24/π。注意这两个函数都是按“向量”设计的doy和lat_deg传 numpy 数组或 pandas Series 都行别在循环里逐行调用那样等于把向量化的优势全扔了。3.4 Ångström-Prescott 主公式与单位换算最后一步最薄但也最容易被人改错# angstrom.py import numpy as np def angstrom_prescott(ra, n, n_max, a0.25, b0.50): ra: 天文辐射 MJ m-2 d-1n: 日照时数 hn_max: 可照时长 N h ratio np.clip(n / n_max, 0.0, 1.0) # 相对日照强制落在 [0,1] return (a b * ratio) * ra def to_wm2(rs_mj): MJ m-2 d-1 转 W m-2 return rs_mj * 1e6 / 86400.0ratio的裁剪不是可选项。日照计有阈值和时间响应误差实测日照时数偶尔会略超 N比值变成 1.02 这种数字套进公式后辐射反而超过 ab 的晴空上限。裁剪到 1 之后最坏情况就是取到晴空值不会出现物理上讲不通的结果。同时我会在批处理里统计被裁剪的行数这个比例本身就是一个很好的数据质量指标——如果某个站一年里有 20% 的天数被裁剪那多半不是天气问题是仪器或统计口径问题。默认系数 0.25/0.50 是 FAO-56 给全球推荐的缺省值用在没有本地标定系数的场景属于合理妥协。但中国幅员跨度大西北干旱区和四川盆地的透射特性差得很远全国一套系数跑下来区域系统偏差能到 10% 以上这一块留到第 6 章展开。4. 批量跑通840 个站点从原始表到结果表的完整流程4.1 数据读取与缺测清洗原始数据最常见的形态是一个大 csv或者按站点切成的几百个小文件。不管哪种读取这一步要把三件事做掉字段类型钉死、缺测码归一、非法值剔除。# batch_run.py片段一 import numpy as np import pandas as pd DTYPES { station_id: string, # 必须字符串否则 000123 会变成 123 lat: float32, lon: float32, elev: float32, ssh: float32, } def load_station_csv(path): df pd.read_csv(path, dtypeDTYPES, parse_dates[date]) # 常见缺测编码32766、999999、-9999统一归一到 NaN df[ssh] pd.to_numeric(df[ssh], errorscoerce) df.loc[df[ssh].isin([32766, 999999, -9999]), ssh] np.nan # 物理上不可能的值直接判缺测 df.loc[(df[ssh] 0) | (df[ssh] 24.0), ssh] np.nan return df.sort_values([station_id, date]).reset_index(dropTrue)station_id用 string 而不是默认的 int64这一点属于典型的血泪经验。站点号普遍带前导零pandas 一推断就变成整数后面按站点分组、和台账表 join 的时候全都对不上而且报错信息不会提示你哪里错了只会静默地少匹配一堆行。缺测码那一行按实际数据源调整不同来源的编码习惯不一样。ssh 24这个判断是兜底一天最多 24 小时超过这个数只可能是编码错误或者单位错有些表用分钟。4.2 向量化计算与结果拼装清洗完之后整表一次性送进计算不要写for index, row in df.iterrows()# batch_run.py片段二 from astro import extraterrestrial_radiation, day_length from angstrom import angstrom_prescott, to_wm2 def compute(df, a0.25, b0.50): doy df[date].dt.dayofyear.to_numpy() lat df[lat].to_numpy() ra extraterrestrial_radiation(doy, lat) # MJ m-2 d-1 n_max day_length(doy, lat) # h rs angstrom_prescott(ra, df[ssh].to_numpy(), n_max, a, b) out df.copy() out[ra_mj] ra.astype(float32) out[daylen_h] n_max.astype(float32) out[rs_mj] rs.astype(float32) out[rs_wm2] to_wm2(rs).astype(float32) return outdayofyear直接从 datetime 列取不用手工拼月份天数闰年自动处理。astype(float32)是体积优化三十万行四列float64 大约 10 MBfloat32 砍一半精度对于辐射量有效数字三到四位完全够。如果后面要做多年序列拼接这个习惯能显著降低内存压力。compute里没有做任何缺测判断——np.nan参与三角函数运算的结果还是nan会自然传播到输出列这正是我们想要的。不要在计算前把 NaN 填 0填 0 等于告诉公式“这天全阴”会污染统计。4.3 目录遍历、结果写出与质检报告批量跑的关键不是算得快是出错时能定位到哪个文件、哪个站# batch_run.py片段三 from pathlib import Path RAW_DIR Path(data/raw) OUT_DIR Path(data/out) def run_all(a0.25, b0.50): OUT_DIR.mkdir(parentsTrue, exist_okTrue) report [] for fp in sorted(RAW_DIR.glob(*.csv)): try: df load_station_csv(fp) res compute(df, a, b) res.to_csv(OUT_DIR / f{fp.stem}_rs.csv, indexFalse) report.append({ file: fp.name, rows: len(res), nan_ssh: int(res[ssh].isna().sum()), rs_max: float(res[rs_mj].max()), clip_ratio: float((res[ssh] / res[daylen_h]).gt(1).mean()), }) except Exception as exc: report.append({file: fp.name, rows: -1, error: repr(exc)}) return pd.DataFrame(report)这里的try/except包住单文件而不是整个循环是为了让一个坏文件不至于中断整批任务。跑完之后把report存成 csv按clip_ratio降序看一眼异常站点会自己浮上来。rs_max那列也有用正常站点全年的日总辐射峰值大致落在 25~32 MJ/m²/d 之间西北高原夏季最高如果出现 45 这种数字多半是纬度字段错了或者日照时数超界没裁干净。5. 避坑与排查五条把结果带偏的真实翻车记录5.1 输入侧的三个数据坑第一坑日照时数大于可照时长。现象是clip_ratio莫名其妙地高某些站高达 15%换算出来的辐射贴着晴空上限。原因通常是日照计阈值设置偏高比如按 200 W/m² 判定、仪器时间响应慢导致日落后还在累计或者统计日界不统一。解决办法分两步代码里对n/N做clip(0,1)兜底业务上把裁剪比例作为质检指标输出超过 5% 的站点回查原始分钟数据处理方式。第二坑纬度字段是字符串形式的度分秒。现象是脚本不报错但 Ra 明显不对或者np.deg2rad收到 object 类型直接抛异常。原因是一些台账表把经纬度存成116°28′这种格式或者拆成了lat_deg、lat_min两列。解决办法是在读取阶段强制float32转换转换失败的行打印出来人工处理如果源数据确实是度分秒写一个小函数d m/60 s/3600先归一。别指望 pandas 猜。第三坑日期解析成 DD/MM 还是 MM/DD。现象是大部分数据看起来正常个别月份 Ra 曲线整体偏移。原因是源表的日期格式是01/02/2020这种含糊写法parse_dates按美式习惯解析成了 1 月 2 日实际是 2 月 1 日。解决办法是pd.read_csv时显式传dayfirstTrue或直接用pd.to_datetime(..., format%Y-%m-%d)指定格式。这个坑最阴的地方在于它不报错结果看上去也对只是四季的相位悄悄挪了。5.2 计算侧的两个公式边界问题第四坑高纬站点冬季整段 NaN。现象是东北、新疆北部站点在 11 月到 1 月出现大量 NaN。原因就是第 3.2 节说的arccos越界-tan(φ)tan(δ)在极夜条件下大于 1。解决办法是np.clip(x, -1, 1)裁剪后 ωs 返回 0Ra 返回 0日长返回 0物理上都对。补充一句中国境内真正意义上的极夜很少但在纬度 50° 上下、冬至前后数值上恰好会顶到边界不加 clip 迟早会遇到。第五坑全国一套系数导致区域系统偏差。现象是整体相关性看着不错但分区统计时西北偏高、西南偏低MBE 能到 2 MJ/m²/d。原因是默认 a0.25、b0.50 是全球平均意义下的取值中国西北干旱区大气透明度高、a 值偏小四川盆地云量大、b 值偏小。解决办法不是简单调一个全国常数而是分区标定具体做法放在下一章。临时凑合的话至少把干旱区和湿润区分别给一套系数误差能砍掉一半。6. 二次开发系数本地化、并行加速与交叉验证6.1 用实测辐射反标定 a、bÅngström-Prescott 是线性形式标定就退化成一次最小二乘回归把实测日总辐射除以同日 Ra 得到纵轴相对日照 n/N 作横轴斜率是 b截距是 a。前提是你手上有辐射观测站的数据哪怕是几十个站也足够按气候区分组拟合。import numpy as np def fit_coefficients(rs_obs, ra, n, n_max, min_days200): rs_obs: 实测日总辐射 MJ m-2 d-1返回 (a, b, r2, n_used) x np.clip(n / n_max, 0.0, 1.0) y rs_obs / ra mask np.isfinite(x) np.isfinite(y) (ra 0) if mask.sum() min_days: return None # 样本太少别硬拟合 x, y x[mask], y[mask] b, a np.polyfit(x, y, 1) # polyfit 返回的是 [斜率, 截距] pred a b * x ss_res ((y - pred) ** 2).sum() ss_tot ((y - y.mean()) ** 2).sum() r2 1.0 - ss_res / ss_tot return float(a), float(b), float(r2), int(mask.sum())min_days200这个门槛是必要的。样本少于两三百天时拟合出来的 b 可能大到 0.7 以上a 甚至是负数虽然数学上成立但物理上讲不通。真遇到这种情况宁可回退到区域默认值也不要硬用。拟合完记得看一眼 ab正常应该落在 0.68~0.80 之间超出这个范围就要怀疑实测辐射数据本身的质量控制。实践中更稳的做法是分气候区拟合比如干旱区、半干旱区、湿润区各一套每套再用邻近辐射站的样本池化。单个辐射站的样本量往往不够池化之后稳定性明显改善。系数表和站点-分区映射表都放进config.py别硬编码在函数里。6.2 按文件并行与内存控制840 个文件如果是按站点切的天然适合进程级并行。pandas 的 GIL 释放不彻底多线程对纯计算帮助有限用ProcessPoolExecutor更直接from concurrent.futures import ProcessPoolExecutor from pathlib import Path def run_parallel(raw_dir, max_workers8): files sorted(Path(raw_dir).glob(*.csv)) with ProcessPoolExecutor(max_workersmax_workers) as pool: results list(pool.map(process_one_file, files, chunksize4)) return pd.DataFrame(results)chunksize4是在任务数和进程数之间找平衡。每个文件只有一两千行任务粒度太细的话进程间通信和序列化的开销会盖过计算本身给每个进程一次发四五个文件的活调度开销摊薄下来就很可观了。max_workers别超过物理核数IO 密集型的场景才需要往上加。跑完之后用pd.concat合并结果时如果总行数上了千万级改成边算边写 parquet 分片全量塞进内存迟早会爆。6.3 交叉验证与结果落地脚本跑通只是开始还得知道误差有多大。我一般会留出有辐射观测的站点把换算结果和实测值对一遍看三个指标MBE平均偏差反映系统误差、RMSE均方根误差反映离散度、R²解释方差。用默认系数时多数区域 R² 能到 0.85 以上RMSE 在 2.5~3.5 MJ/m²/d 区间分区标定后 RMSE 通常能压到 2 MJ/m²/d 以内。这些数字不是拿来交差的是后续改系数时的对照基准——没有基准任何“优化”都无从判断好坏。落地层面的一个具体建议结果表里务必保留ra_mj和daylen_h。很多人只输出rs_mj等到半年后有人来问“这个值靠不靠谱”你已经没法回溯了。留着中间量任何一个可疑数值都能当场复算这对二次开发的价值远大于表大那么几 MB。从那以后我每次跑批量换算都会先抽三个站、跨四季各取十天把 Ra 的日变化曲线画出来看一遍确认相位和量级都对再放开跑全量——这个习惯帮我拦下过两次时区和日界错位。希望帮到你。本文还有配套的精品资源点击获取