全国840站日照时数批量转日总太阳辐射Python实践

简介面向气象、气候与太阳能资源评估场景的Python数据处理脚本可将全国840个气象站点的日照时数观测值换算为日总太阳辐射适合科研人员、数据分析爱好者及具备基础Python能力的学习者使用。压缩包内共3个文件以2个CSV数据文件和1个Python脚本为主整体约275KBCSV分别承载输入观测数据与转换后的输出结果脚本承担读取、清洗、计算和结果保存等流程。已有6351人学习说明该工具在相关数据处理练习中具有一定参考价值。脚本围绕日期、经纬度与日照时数等字段展开包含缺失值处理、经验公式转换及结果导出思路读者可据此理解日照时数与日总太阳辐射之间的能量转换关系并在此基础上扩展站点分布图、批量处理和误差检查等分析用于气候变化、能源利用与气候模型构建的初级实践。1. 全国840个气象站点日照时数转日总太阳辐射一份拿来就能跑的Python代码手头拿到的是全国840个气象站点的逐日日照时数单位是小时可模型要的是日总太阳辐射单位是 MJ/m²。中间隔着天文辐射、日照百分率、经验系数三层换算单站手算还能忍840 个站点乘上多年逐日记录Excel 直接躺平。这套 Python 代码干的就是这件事把日照时数批量转成日总太阳辐射输出成可以直接进模型的规整表格。它适合三类人做农业气象和蒸散发的做光伏资源初筛的以及手上只有气象站观测、拿不到辐射产品的研究者。脚本只依赖 numpy 和 pandas核心逻辑不到两百行参数放在配置里换站点网、换系数都不用动主流程。2. 技术选型与数据链路从一份站点表推到可复现的辐射序列日照时数是个很“台站”的变量观测历史长、覆盖站点多、质量控制相对成熟。但它不是辐射必须经过一次经验换算。这一章先把“为什么这么算”和“数据怎么组织”讲清楚后面再看代码就不容易迷路。2.1 为什么走 Ångström-Prescott而不是直接下载再分析辐射产品常见的替代方案是拿再分析辐射产品直接插值到站点位置。这条路能走通但有几个现实麻烦一是格点分辨率粗山区和地形复杂地带偏差明显插值到站点后误差不好交代二是数据量大几百个站点、几十年逐日下载和存储是额外成本三是版本更新后序列可能变化复现性差。Ångström-Prescott 形式简单Rs (a b·n/N)·Ra。其中Ra是日地天文辐射只跟纬度和日序有关纯几何量算得准n/N是日照百分率直接来自台站观测a、b是经验系数FAO-56 给了一组通用值 0.25 和 0.50。这套方法的优点不是“最准”而是可复现、可解释、可标定——你若手上有辐射观测可以自己拟合 a、b 把精度拉回来没有的话用通用值误差量级也是业内公认可接受的。我一般会这么定有辐射观测站的区域优先标定没条件的用通用值但一定要在输出里带上“系数来源”这一列。不然半年后回头看结果完全想不起当时用的是哪套系数这就是典型的后悔药场景。另一个容易忽略的点是Ra只跟纬度和日序有关跟经度、海拔基本无关海拔对辐射的影响主要体现在大气透过率上经验公式这一层通常不单独修正。所以 840 个站点里纬度相同或相近的站点Ra的序列是同一条只是日期对齐方式不同。这个性质是后面做向量化的关键。2.2 脚本的目录结构与每个文件的职责工程目录按“计算与 IO 分离”来切这样标定系数、换数据源都不会污染核心公式sunshine2rs/ ├── config.yaml # 路径、系数、单位开关 ├── data/ │ ├── stations.csv # 站点元数据: id, lat, lon, elev │ └── sunshine/ # 日照时数, 按年份分文件 ├── src/ │ ├── radiation.py # 纯计算函数, 不碰文件 │ ├── io_utils.py # 读取、清洗、写出 │ └── calibrate.py # a、b 系数最小二乘标定 ├── run.py # 命令行入口 └── tests/ └── test_radiation.py # 已知输入的回归测试radiation.py里只放函数输入是 numpy 数组输出也是数组不读文件、不写日志。这么做的好处是测试极其轻量给一组固定的纬度、日序、日照时数断言输出等于手算值就行。io_utils.py承担所有脏活——编码、缺失值标记、单位换算、列名对齐。真正出问题时你能立刻判断是“数据脏”还是“公式错”不用在一个三百行的脚本里靠 print 大法定位。config.yaml只放四类东西输入输出路径、a/b系数、日照时数的原始单位小时还是 0.1 小时、输出单位MJ/m²/d 还是 W/m²。把这些从代码里抽出来是为了换区域跑的时候不用改代码、不用重新测试。2.3 三类输入数据的字段约定数据链路一共三张表字段约定必须先定死表关键字段类型说明stations.csvstation_idstr站号一律按字符串存避免前导零丢失lat / lonfloat十进制度北纬为正elevfloat米当前版本只做记录不参与计算sunshine_YYYY.csvstation_idstr必须和站点表能对上datedate格式统一为 YYYY-MM-DDssdfloat日照时数单位由配置指定输出表station_id, date与输入对齐Ra / N / ratio / Rsfloat中间量全部保留便于排查站号用字符串是血泪经验。有些站号带字母后缀有些数据源存成数值后在导出环节被截断等到 merge 的时候发现对不上几百个站点静默丢失还不报错。用字符串存、merge 后用indicatorTrue检查_merge列是最省事的防线。日照时数的原始单位也必须显式声明。台站数据里既有以小时为单位的浮点值也有以 0.1 小时为单位的整数值混在一起算结果会差十倍而且不会报错——数值看起来仍然“像辐射”只是整体偏大。这种错误比程序崩溃难查得多。3. 核心算法拆解日照时数如何一步步变成 MJ/m²这一章把公式拆成四步每一步都给可执行代码。看完你应该能拿一张纸自己推一遍也能判断输出结果是否离谱。3.1 日地天文辐射 Ra所有换算的基准面Ra表示在没有大气的假设下水平面上一天能接收到的太阳辐射总量单位 MJ/m²/d。计算过程是三个中间量加一个合成式import numpy as np GSC 0.0820 # 太阳常数, MJ m-2 min-1 def extraterrestrial_radiation(doy, lat_deg): 计算日地天文辐射 Ra (MJ m-2 d-1) doy : 年日序, 1-366, 支持 numpy 数组 lat_deg : 纬度, 十进制度, 北纬为正 lat np.deg2rad(lat_deg) # 日地相对距离倒数平方, 无量纲 dr 1 0.033 * np.cos(2 * np.pi * doy / 365) # 太阳赤纬, 弧度 delta 0.409 * np.sin(2 * np.pi * doy / 365 - 1.39) # 日落时角, 弧度; clip 处理极昼极夜导致的越界 x np.clip(-np.tan(lat) * np.tan(delta), -1.0, 1.0) ws np.arccos(x) # 合成 Ra Ra (24 * 60 / np.pi) * GSC * dr * ( ws * np.sin(lat) * np.sin(delta) np.cos(lat) * np.cos(delta) * np.sin(ws) ) return Ra, ws, dr, deltadr是日地距离修正项年中在 0.97 到 1.03 之间摆动看起来不起眼但在冬季高纬度站点能带来百分之几的差异。delta是太阳赤纬用 FAO-56 的经验式即可精度足够。ws是日落时角必须做clip否则在极昼极夜时段arccos的入参超出[-1, 1]会返回nan然后整条序列被污染。系数24*60/π的来历是把“分钟”和“弧度”两个单位统一到“天”和“MJ/m²”。这个系数不要记直接照抄但要知道它为什么在那儿——单位不统一是这类公式最常见的翻车点。3.2 理论日照时数 N 与日照百分率 n/N有了日落时角理论日照时数N就是白昼长度def daylight_hours(ws): 由日落时角(弧度)换算理论日照时数 N (小时) return 24 / np.pi * wsws取 π 时N等于 24 小时取 0 时N等于 0物理上是对的。N是分母所以必须先处理N 0的情况。极夜站点虽然少见但全国 840 个站点里高纬度站点在冬季确实可能接近 0直接相除会得到inf。日照百分率的计算要同时处理两个边界def sunshine_ratio(ssd, N): 日照百分率 n/N, 结果裁剪到 [0, 1] safe_N np.where(N 1e-6, N, np.nan) ratio ssd / safe_N # 观测的日照时数理论上不可能超过白昼长度, 超出视为记录异常 ratio np.clip(ratio, 0.0, 1.0) return ratio为什么上界要裁到 1因为台站观测的日照时数偶尔会因为仪器灵敏度、记录舍入等原因略微超过N比如N算出来 13.2 小时、观测记录 13.5 小时。如果不裁ratio大于 1Rs会超过(ab)·Ra物理上讲不通裁到 1 相当于认为“全天有日照”这是可接受的保守处理。负值和缺测标记必须在更早的清洗环节处理掉。常见的缺测标记是-999、9999、空字符串如果漏进计算ratio会被裁到 0Rs变成a·Ra看起来“有值但偏小”反而比报错更难发现。3.3 a、b 系数最容易拍脑袋的一个参数FAO-56 推荐的通用值是a 0.25、b 0.50含义是全阴天n/N 0时地面接收的短波辐射约为天文辐射的 25%晴天n/N 1时约为 75%。这两个数字在很多气候区都能用但偏差不稳定。国内常见的做法是按气候区取标定值比如干旱区b偏大、湿润区a偏小。手上有辐射观测站的话最稳妥的办法是自己拟合def fit_angstrom(ratio, rs_obs, ra): 用辐射观测标定 a、b ratio : 日照百分率 n/N rs_obs : 观测日总辐射 MJ m-2 d-1 ra : 天文辐射 MJ m-2 d-1 返回 (a, b)并保证 a b 1 y np.asarray(rs_obs) / np.asarray(ra) # 晴空指数 x np.asarray(ratio) mask np.isfinite(x) np.isfinite(y) slope, intercept np.polyfit(x[mask], y[mask], 1) a, b float(intercept), float(slope) if a b 1.0: # 物理约束兜底 b 1.0 - a return a, b这段代码没有做交叉验证实际用的时候至少要留一年数据做独立检验算一下 MBE 和 RMSE。系数拟合出来的数值如果a为负或b大于 1基本可以判定样本里有异常值或者Ra单位错了不要急着用。另一个常见误区是拿“日照时数”直接当n/N用。日照时数的单位是小时N的单位也是小时两者相除才是无量纲的百分率。有些数据源直接提供“日照百分率”单位 %这时要除以 100别直接塞进公式。3.4 向量化落地840 个站点一次算完840 个站点乘上三十年逐日大概是九百万行量级。用for循环一站一天地算Python 层面至少几分钟用 pandas 整表向量化通常十几秒能跑完。关键在于“广播”把站点纬度 merge 到日照长表上让每一行都带上自己的纬度然后一次性对整列做 numpy 运算。def compute_rs(df, a0.25, b0.50, ssd_colssd, lat_collat): 整表计算日总太阳辐射, 要求 df 已含 date 与 lat 列 doy df[date].dt.dayofyear.to_numpy() lat df[lat_col].to_numpy(dtypefloat) Ra, ws, _, _ extraterrestrial_radiation(doy, lat) N daylight_hours(ws) ratio sunshine_ratio(df[ssd_col].to_numpy(dtypefloat), N) df df.copy() df[Ra] Ra df[N] N df[ratio] ratio df[Rs] (a b * ratio) * Ra return df四个中间量全部保留在输出里是为了后面自检。Rs出问题时你能一眼看出是Ra偏了、N偏了还是ratio被裁了。把中间量丢掉只留Rs等于把整个计算过程变成黑匣子。向量化的代价是内存九百万行乘上几个 float64 列大约几百 MB普通工作站扛得住。如果数据量再大一倍就按年份分批跑每年一个 parquet 文件最后再拼——这个改动只需要在入口脚本里加一层循环不影响radiation.py。4. 代码逐段走读一条记录从读入到写出的完整链路这一章按实际运行顺序走一遍读入、清洗、合并计算、写出校验。每一步都给代码你可以照着改。4.1 读入与清洗把宽表和长表统一成日期索引台站数据的形态不统一有的是“一行一天、多列站点”的宽表有的是“一行一个站点一天”的长表。脚本统一往长表靠因为长表更适合 merge 和分组。import pandas as pd def load_stations(path): st pd.read_csv(path, dtype{station_id: str}) st[station_id] st[station_id].str.strip() st[lat] pd.to_numeric(st[lat], errorscoerce) st st.dropna(subset[lat]) # 纬度合法性检查, 超出范围的站点直接剔除并记录 bad st[(st[lat] -90) | (st[lat] 90)] if len(bad): print(f[warn] 剔除纬度异常站点 {len(bad)} 个: {bad[station_id].tolist()[:5]}) return st[(st[lat] -90) (st[lat] 90)].copy()dtype{station_id: str}这一句必须写不能用默认推断。默认推断会把纯数字站号读成 int64合并时一旦另一张表是字符串就会得到全空的 merge 结果而且不报错。日照时数的清洗重点是缺测标记和单位def load_sunshine(path, unit_scale1.0): df pd.read_csv(path, dtype{station_id: str}) df[date] pd.to_datetime(df[date], errorscoerce) ssd pd.to_numeric(df[ssd], errorscoerce) # 统一缺测标记为 NaN ssd ssd.mask(ssd.isin([-999, -9999, 9999, 99999])) ssd ssd.mask(ssd 0) df[ssd] ssd * unit_scale # 0.1 小时制数据在此处乘 0.1 return df.dropna(subset[date, station_id])unit_scale从配置读进来。这一步看着平淡但它决定了整个结果的量级是否正确属于“改一个数字、全盘皆错”的地方务必在配置里写清楚注释。4.2 合并与计算站点纬度广播与逐日循环的取舍合并是整条链路最容易出事的一步。我的习惯是显式检查_merge列def build_daily_table(sunshine, stations): df sunshine.merge( stations[[station_id, lat, lon]], onstation_id, howleft, indicatorTrue ) missing df[df[_merge] left_only][station_id].unique() if len(missing): print(f[warn] {len(missing)} 个站号在站点表中不存在: {missing[:5]}) df df[df[_merge] both].drop(columns_merge) return df.sort_values([station_id, date]).reset_index(dropTrue)left_only的站号要么是站号写错了要么是站点表不完整。直接丢掉会让输出“看起来正常但少了站点”而 840 个站点丢十几个肉眼根本发现不了。打印出来是成本最低的兜底。合并完成后再调compute_rs整表一次算完。这里不需要按站点分组循环因为公式里每个站点只用到自己的纬度而纬度已经广播到每一行上了。如果确实需要分组比如不同气候区用不同系数用groupby加apply会比for循环快但要注意apply返回的索引结构。更稳的写法是给每行先打上“系数组”标签然后用map把a、b映射成列再做向量化计算——这样完全不需要分组。4.3 写出与自检三条勾稽关系卡住大部分错误输出前加一段自检能抓出八成以上的错误def sanity_check(df, tol1e-6): issues {} # 1. 日照百分率必须落在 [0,1] issues[ratio_out_of_range] int( ((df[ratio] -tol) | (df[ratio] 1 tol)).sum() ) # 2. Rs 不应超过 (ab)*Ra, 用 0.8 做宽松上界 issues[rs_gt_ra] int((df[Rs] 0.81 * df[Ra] tol).sum()) # 3. 观测日照时数不应超过理论白昼 issues[ssd_gt_N] int((df[ssd] df[N] 0.5).sum()) return issues三条检查分别对应三类问题ratio越界说明清洗没做好Rs过大说明系数或单位有问题ssd N说明观测记录本身异常或者N算错了。阈值给 0.81 是因为ab最大 0.75留一点余量给系数标定后的情况。输出格式建议用 parquet字段类型不会在反复读写中漂移。CSV 在中文路径、编码、float精度上都容易出岔子尤其是站点名带中文的时候。5. 避坑与排查跑出来的辐射值不对时先按这个顺序查下面几条都是实际跑数据时踩过的按“输入、计算、输出”三段排。每条都按现象、原因、解决的顺序写照着对号入座就行。5.1 输入侧缺测、单位、站号对不齐现象一输出结果整体偏大十倍左右但空间分布看起来很正常。原因日照时数原始单位是 0.1 小时制配置里的unit_scale没改直接按小时算了。ratio被裁到 1 之后Rs几乎恒等于0.75·Ra看起来“处处晴天”。 解决先抽样看原始数据的取值范围正常日照时数在 0 到 14 之间出现 0 到 140 的量级就是 0.1 小时制。把unit_scale改成 0.1重跑并对比ratio的分布。现象二输出少了十几个站点没有任何报错。原因站号在两张表里类型不一致一张是字符串、一张是整数merge 后left_only被静默丢弃。 解决所有读入都强制dtype{station_id: str}merge 时加indicatorTrue把left_only的站号打印出来。现象三某些站点整年Rs全是NaN。原因缺测标记-999没有过滤被当成真实值参与运算后又触发了NaN传播或者日期列里有非法格式to_datetime转换失败变成NaT。 解决在清洗函数里统一把-999 / -9999 / 9999 / 99999映射为NaN并对date列检查isna().sum()大于 0 就打印前几行原始记录。5.2 计算侧角度弧度、闰年、极昼极夜现象四所有站点的Ra都在冬季出现nan。原因高纬度站点在冬季-tan(lat)·tan(delta)超出[-1, 1]arccos返回nan。 解决对入参做np.clip(x, -1, 1)。裁剪后ws变成 0 或 π对应极夜或极昼物理上是合理的。注意裁剪之后N可能为 0ratio计算时要用np.where保护分母。现象五2 月 29 日附近的Ra比前后几天有轻微跳变。原因公式里用的是 365 天周期闰年时doy在 3 月 1 日之后整体偏移一天导致delta和dr有微小错位。 解决这是方法本身的近似量级在千分之几通常可以接受。如果要求严格可以把doy换算成“年内比例”而不是整数日序或者干脆用 366 做分母并接受另一侧的偏差。别在这里过度纠结先确认误差是否影响下游结论。现象六Rs出现负值或极端大值。原因纬度传成了弧度又做了一次deg2rad或者站点表里的纬度符号搞反南半球站点写成正值。 解决在extraterrestrial_radiation入口加断言lat_deg必须在[-90, 90]输出前检查Rs的分位数正常日总辐射在 0 到 40 MJ/m²/d 之间超过这个范围基本可以判定为输入错误。5.3 输出侧单位换算与物理边界现象七结果交给下游模型后蒸散发算出来大得离谱。原因Rs单位是 MJ/m²/d但下游模型期望的是 W/m²量纲没换。 解决在配置里加output_unit开关MJ 转 W/m² 乘1e6/86400 ≈ 11.574如果下游要的是等效蒸发量mm/d除以潜热系数 2.45。换算集中在一个函数里做别散落在各处。现象八同一站点不同年份的Rs序列拼起来后在年初有明显断层。原因按年分文件跑doy是每年重新从 1 开始算的本身没问题断层通常来自date解析时年份被推断错或者文件排序时按字符串排了date列。 解决所有date读入后立刻to_datetime排序统一用sort_values([station_id, date])。写文件时不要用字符串日期用真正的日期类型。6. 二次开发换成自己的站点网、自己的系数再补一套校验习惯这份代码最容易改的地方恰恰也是最该改的地方。换站点网只需要替换stations.csv保证station_id和lat两列存在且纬度是十进制度。经度和海拔在计算里不参与但建议保留——以后接大气透过率修正或者做空间插值时用得上。如果站点量从 840 涨到几千只需要在入口按年份或按站点分批读compute_rs本身不用动。换系数分两种情形。有辐射观测的用fit_angstrom拟合然后按气候区分组落表分组依据典型做法注意点按气候区每个区一套 a、b分区边界要固定别年年改按站点每站一套 a、b观测年限短的站点容易过拟合全局一套直接用 0.25 / 0.50误差可接受但要在文档里写清楚没有观测的老老实实用通用值但要在输出表里加一列coef_source取值fao56或local。这列以后会救你很多次。加时间聚合很简单groupby之后求和或求平均即可但要注意日总辐射的月值应该是“月内日值之和”而不是“日值平均”单位是 MJ/m²/月。这个区别在给作物模型供数时经常被搞混一旦搞混整个生育期的辐射累积会差三十倍。真正值得投入精力的是校验环节。我现在固定做三件事一是随机抽三个站点用手算核对一整月的Ra确认公式没抄错二是把Rs/Ra的分布画出来正常应该集中在 0.4 到 0.7 之间出现双峰或者长尾就说明ratio的清洗有问题三是把输出和邻近的辐射观测站做一次年尺度对比看偏差是否稳定。前两条是自检第三条才是真正意义上的验证。说个教训。有次批量跑完结果看起来一切正常Rs的均值、分位数都在合理区间我就直接交出去了。两周后下游反馈说某几个站点的季节曲线反了——查下来是站点表里那几个站的纬度符号写错了lat存成了负值Ra在南半球模式下算出来刚好“数值合理但季节相反”。从那以后我每次读站点表都强制走一遍纬度范围检查和符号分布统计把“有多少站在南半球”打印出来跟预期比对。这套习惯没什么技术含量但它挡住的问题比任何优化都多。希望帮到你。本文还有配套的精品资源点击获取