简介面向全国840个气象站点的日照时数转日总太阳辐射Python代码专为气象科研人员、环境科学研究者及有一定Python基础的数据分析学习者设计。脚本覆盖数据读取、缺失值与异常值清洗、基于林格曼方程等经验公式的辐射换算、结果输出以及matplotlib可视化等完整流程并包含错误处理与性能优化可应用于气候变化研究、太阳能资源评估和气候模型构建等实际场景。压缩包共3个文件包括1个Python脚本和2个csv样例数据分别存放输入日照时数与转换后的输出结果便于对照验证整体仅275KB。已有6336人学习下载。该代码既可充当自动化处理840个站点日辐射数据的实用工具也可作为学习pandas数据操作、科学计算与可视化的实践范例帮助读者深入理解日照时数到日总太阳辐射的能量转换机制。 840个气象站点、几十年的日照时数记录要一口气转换成日总太阳辐射这个需求听起来就是把一个公式套进循环里跑一遍但真正动手时会发现数据不是攥在一个文件里的系数不是一个固定值通吃全国的算出来的数字也不是靠眼睛看就能判断靠不靠谱的。这篇文章我就把这套流程完整摊开讲包括数据预处理、天文辐射计算、Angström-Prescott模型的落地写法以及840个站点批量处理时我踩过的各种坑。适合做光伏资源评估、农业气候区划、建筑能耗模拟的朋友参考。1. 为什么要用日照时数来推太阳辐射1.1 直测辐射数据的缺口远比想象中大太阳总辐射的直接观测依赖日射强度计这类仪器本身不便宜而且对维护要求非常高石英罩要定期清洁、灵敏度会漂移、内部的干燥剂失效就得换一不留神数据质量就崩了。所以全国真正长期稳定发布总辐射观测的站点数量并不多要做省级网格、全国尺度的资源图谱光靠这些站点根本撑不起空间分辨率。相比之下日照时数的观测几乎每个国家级气象站都在做。日照计的记录相对皮实数据序列长、站点密度高全国840个站点的日照时数资料基本上可以把各个气候亚区都覆盖到。所以很多工程场景下日照时数成了太阳辐射的“代理变量”用经验模型把它换算成总辐射属于一种非常成熟的工程妥协。1.2 日照时数转辐射的经典模型日照时数和太阳辐射之间最经典的桥梁就是Angström-Prescott回归方程H H0 × (a b × n/N)其中H是目标日总太阳辐射H0是大气层顶的日天文辐射量也叫地外辐射n是当天实测日照时数N是当天理论可照时数也就是白昼时长a和b是经验系数。模型的核心逻辑很直观日照占比n/N越高云越少实际到达地面的辐射占天文辐射的比例就越高。这个模型从二十世纪二十年代提出至今经过无数区域修正依然在光伏、气象、农业领域大量使用。它最大的优点是只依赖日照时数一个常规观测要素就能计算而且只要标定好本地系数精度相当可观。业界也尝试过引入气温、湿度、降水做更复杂的机器学习模型但数据获取成本和可解释性都不如这个经典公式。840个站点要做统一批处理用这个模型是最合适的起点。2. 840个站点的原始数据先整理成能跑模型的样子2.1 数据文件长什么样读进来要注意什么国内数据集常见的形式是“一站一个文件”文件名就是区站号文件内是类似这样按空格或者固定宽度排列的文本57083 2015 01 01 0.0 57083 2015 01 02 8.2 57083 2015 01 03 -9999列的含义通常是区站号、年、月、日、日照时数单位小时。不同来源的数据缺失值编码不一样常见的有-9999、9999、32766这一类的占位符读进来之后第一步就是把它们替换成NaN。还有一个极其容易踩的细节有些数据集里日照时数会乘以0.1存储也就是说表里存的是8.2实际要除以10才是0.82小时不看元数据直接算结果会离谱到天际。先准备两份输入一份是站点元信息区站号、纬度、经度、海拔一份是长表观测记录。下面是我惯用的读法import pandas as pd import numpy as np # 站点信息 stations pd.read_csv( stations_840.csv, encodinggbk, dtype{station_id: str}, ) # 日照时数长表所有站点合并后 # 假设列名station_id, year, month, day, sunshine_hour obs pd.read_csv( sunshine_hour_all.csv, sepr\s, names[station_id, year, month, day, sunshine_hour], dtype{station_id: str}, ) obs[date] pd.to_datetime( obs[[year, month, day]] ) obs[sunshine_hour] obs[sunshine_hour].replace( [-9999, 9999, 32766], np.nan ) obs obs.dropna(subset[sunshine_hour]) obs obs[(obs[sunshine_hour] 0) (obs[sunshine_hour] 24)]最后这个范围过滤很关键日照时数不可能小于0也不应该超过24小时一些脏数据比如28.5、-3.2这时候就暴露了。2.2 合并站点信息和观测记录模型需要用到纬度和经度所以要先把站点元信息合并到观测长表里。合并之后还有一个动作不能省检查每个站点经纬度是否有缺失、是否有重复的日期记录。df obs.merge(stations, onstation_id, howleft) df df.dropna(subset[lat, lon]) # 检查同一个站点是否有重复日期 dupe df.duplicated(subset[station_id, date], keepFalse) if dupe.sum(): print(f发现 {dupe.sum()} 条重复日期记录请先人工核对) df df.drop_duplicates(subset[station_id, date])840个站点合并完行数通常是几百万级但这一步在pandas里也就是几秒的事。真正花时间的是检查那些“看起来成功但实际失败”的站点比如某些站点某一年文件缺失、某些站点经纬度为零值这些不查清楚后面算出来的辐射值就是一堆心里没底的数字。3. 辐射估算模型的核心计算一步一步拆给你看3.1 日地外辐射H0的计算Angström-Prescott里的H0不是一个常数它随着日期、纬度变化是天文位置决定的。计算H0需要三个量日地距离修正系数dr、太阳赤纬角δ、日落时角ωs。为了在840个站点上做向量化计算我建议全部用numpy批量处理不要写for循环逐步调用math库。先计算一年中的第几天然后算日角doy df[date].dt.dayofyear.to_numpy() # 日角弧度 gamma 2 * np.pi * (doy - 1) / 365.0 # 日地距离修正系数 dr 1 0.033 * np.cos(gamma) # 太阳赤纬角Spencer公式弧度 dec ( 0.006918 - 0.399912 * np.cos(gamma) 0.070257 * np.sin(gamma) - 0.006758 * np.cos(2 * gamma) 0.000907 * np.sin(2 * gamma) - 0.002697 * np.cos(3 * gamma) 0.00148 * np.sin(3 * gamma) )太阳赤纬角还可以用更简单的Cooper公式近似但Spencer公式的精度对日尺度计算来说更有保障尤其在春秋分附近两者差异能达到0.5度以上这对辐射值的影响不是可以忽略的。日地外辐射总量的标准公式是H0 (24×3600×Gsc/π) × dr × (ωs×sinφ×sinδ cosφ×cosδ×sinωs)其中Gsc是太阳常数工程计算里取1367 W/m²φ是纬度弧度。最终结果除以10^6转成MJ/m²/day。用numpy直接这样写lat_rad np.deg2rad(df[lat].to_numpy()) Gsc 1367.0 omega_0 np.arccos( np.clip(-np.tan(lat_rad) * np.tan(dec), -1.0, 1.0) ) H0 ( (24.0 * 3600.0 * Gsc / np.pi) * dr * ( omega_0 * np.sin(lat_rad) * np.sin(dec) np.cos(lat_rad) * np.cos(dec) * np.sin(omega_0) ) / 1e6 )这里有一个非常重要的细节-tan(lat)*tan(dec)在数学上可能落在[-1,1]区间之外高纬度地区尤其常见如果直接丢进arccos会返回NaN。用np.clip把输入夹到[-1,1]再计算能避免一整个站点因为冬季极端日期而全部报废。3.2 理论白昼时长N的计算白昼时长N本质上由日落时角直接换算卫星轨道上太阳中心在地平线上的持续时间单位是小时N (24/π) × ωs对应到代码就是daylight_hours 24.0 / np.pi * omega_0注意N的物理含义极昼情况下ωs等于π弧度N就是24小时极夜情况下ωs是0N就是0。由于前面已经做了clip这两个极端值在计算上天然成立不需要额外判断。但后续算n/N的时候N等于0的那一天会变成除零需要单独处理。3.3 a和b系数怎么取才不至于瞎算这是整个流程里最容易出争议的地方。Angström-Prescott方程的a、b系数并不是普适常数不同气候区、不同季节差异很大。全国840个站点横跨湿润、半干旱、干旱、高原等多个气候区用一个固定系数套全国误差在部分地区会非常难看。文献里的经典取值月尺度上很多研究用a0.18、b0.55这个组合在不少区域适用但并不全球通用。日尺度的估算通常推荐另一个量级。我做实际项目时习惯按气候区给几套参考系数然后再用靠近站点的辐射观测站做校验调整。适用区域ab备注华北平原、东北中部0.200.52春季风沙天误差偏大长江中下游0.180.50梅雨季容易低估四川盆地、云贵东部0.200.44多云地区日照-辐射关系弱青藏高原0.250.48建议再做海拔修正西北干旱区新疆、甘肃0.220.55晴天比例高相关性好如果你有站点附近辐射站的实测资料强烈建议做一次本地回归把自己地区的a、b直接标定出来。方法很简单把实测辐射除以H0得到y把日照占比n/N作为x做一元线性回归截距就是a斜率就是b。这套自标定流程在代码里也就是np.polyfit两行的事但效果远好于随便套文献值。4. 全国840个站点的批处理实现4.1 向量化计算告别840个for循环数据量在百万级时最直观的想法是for循环遍历每个站点然后内层再遍历每一天840个站点跑下来可能要等一阵子。实际上完全可以一步到位把所有站点所有日期的数据一次性丢进numpy数组做向量化计算几百万行也只需几秒到十几秒。把经纬度、日期、日照时数都合并到同一个DataFrame后按上面的方式算出H0和daylight_hours然后直接算估算辐射sunshine_hour df[sunshine_hour].to_numpy() # 日照占比注意极夜N0的情况 ratio np.divide( sunshine_hour, daylight_hours, outnp.zeros_like(sunshine_hour, dtypenp.float64), wheredaylight_hours 0, ) # 给每个站点挂上对应区域的a、b示例全国统一初值 a_val 0.20 b_val 0.52 ghi_estimate H0 * (a_val b_val * ratio) # 物理量约束辐射值不为负 ghi_estimate np.maximum(ghi_estimate, 0.0) df[ghi_estimate] ghi_estimatenp.divide里的where参数是个容易被忽略但是极其关键的操作它只允许白天时长大于0的日期参与除法剩下的地方直接用0填充省去一次mask处理。这样写出来的代码没有循环也没有额外的内存爆炸风险。4.2 输出结果的组织方式算完的结果建议不要只憋在一个超大CSV里后续按站点、按日期切片取数会很痛苦。我习惯按站点拆分输出一个站点一个文件或者整体存成parquet格式。下面这个写法实用且不容易出错import os os.makedirs(output/ghi_daily, exist_okTrue) for station_id, g in df.groupby(station_id): out g[[date, ghi_estimate]].sort_values(date) out.to_csv( foutput/ghi_daily/{station_id}.csv, indexFalse, float_format%.2f, )如果你后续想要直接做全国面板分析推荐用df.to_parquet(ghi_daily.parquet)一次性落盘Polars、Dask读取都比CSV快一个量级。我个人一般是两种都存CSV用来给业务同事看parquet留给自己做二次计算。5. 结果验算不要让看起来合理的数字骗了你5.1 用邻近辐射观测站做基准对比算完辐射值之后最忌直接拿去画图甚至写报告。840个站点里可能就有那么几十个站在数据质量上有隐藏问题必须做验证。最有效的办法是找到那些同时具备日照时数和实测总辐射的站点把模型估算值和实测值放在一起看偏差。抽样选十几个覆盖不同气候区的站点计算平均绝对误差MAE、均方根误差RMSE和相对偏差from sklearn.metrics import mean_absolute_error, mean_squared_error # est: 模型估算值, obs: 实测辐射值 mae mean_absolute_error(obs, est) rmse mean_squared_error(obs, est, squaredFalse) mape np.mean(np.abs((obs - est) / obs)) * 100 print(fMAE: {mae:.2f} MJ/m2/day) print(fRMSE: {rmse:.2f} MJ/m2/day) print(fMAPE: {mape:.1f}%)日尺度的Angström-Prescott模型误差到15%~25%都不算意外这是模型的正常波动除非你想做融资级别的精算否则不必太焦虑。但如果你发现某个站点MAPE超过了40%大概率不是模型问题而是这个站点的日照时数或辐射观测本身存在系统误差这时候要去查原始数据而不是硬调系数。5.2 通过物理自检快速淘汰异常在没有实测辐射站对照的区域至少要过一遍物理合理性检查。几个最实用的判断点日辐射估计值不应该超过同一天H0的值更不可能为负值。晴天日照时数接近N时日辐射估算值理论上应在H0的60%~75%区间如果出现接近H0甚至超过H0的数值说明系数设置或数据单位出了问题。同一站点相邻日期的辐射值不应出现剧烈震荡日照时数从0跳到12辐射值从2跳到25是合理的但日照时数都是10的连续两天辐射值却从20掉到5就要警惕。还可以按站点查看年总量均值。全国除四川盆地、贵州一带外年总辐射量多数在4500~6500 MJ/m²/year这个区间如果某站点算出来年总量只有2000或者冲到8000不要怀疑气候特殊先去检查数据。6. 我踩过的几个坑以及对应的解决办法6.1 arccos函数越界导致整年数据报废这是所有踩坑里影响面最大的一个。没加clip之前我把全国数据一次性丢进去算结果青藏高原和东北北部一大批冬季日期的H0变成了NaN一侧是坦坦荡荡的报错另一侧是几万个NaN在DataFrame里安静潜伏。如果后面直接groupby聚合这些NaN会自动被跳过你根本察觉不到直到画图时才发现一片空白。用np.clip把三角函数的输入限制在[-1,1]是唯一可靠的解法。数学上-tanφ·tanδ本来就该落在这个区间超出部分完全是由浮点误差和极端纬度造成的伪异常clip掉不会丢信息。6.2 极昼极夜地区除以零高纬度站点在冬至前后可能出现N0的日子此时不管实际日照时数是多少日照占比n/N都会变成inf或者NaN。用普通除法pandas会给出inf和警告后面再算辐射值就全部被污染。我习惯用np.divide的where参数处理原因前文提过它可以把N0的日子直接置为0占比对应的辐射计算结果也会被修正为合理小值。切记不能把这几天简单删掉否则站点年总量偏低的假象会误导后续分析。6.3 月份日界与日照计的“太阳日”差异日照时数观测在很多台站是以真太阳时为日界的而常规气象数据文件的日期是按北京时间日界归档的。这一点在大尺度统计上影响很小但在日出、日落时间极端的高纬地区冬季可能出现少数日期错位。表现为个别站点在日照时数为0的冬季日前后两天数值突变。处理思路不复杂如果你发现某站点冬季数据在日界附近有明显跳变可以尝试用三日滑动平均做平滑或者在月尺度上重新做聚合把这种边界噪声稀释掉。单独揪出某一天的精确对位成本太高对绝大多数工程需求没有必要。6.4 a、b系数“一招鲜”带来的区域失真早期我做全国批处理时图省事全量用同一组a、b系数结果西北地区估算值明显偏高四川盆地明显偏低。后来把全国分成几个气候区分别回归才把整体误差压下来。这提醒我一个道理模型的公式可以统一但参数绝不能出厂后就不动了。只要是做全国尺度的项目至少要按气候区或省级行政区给参数分档并且在成果里注明每个区域用的参数来源别让下游使用者在汇报时无据可查。这套代码跑通之后后续往里面加新站点只需要在stations表里增加经纬度和海拔重新跑一遍就出结果了。如果你需要进一步做月尺度、季尺度甚至年尺度分析对日结果按站点聚合即可pandas的groupby加resample一条链走下来顺滑得很。我实际做下来最大的体会是真正费时间的永远不是模型而是看不完的数据质量和没完没了的边界情况处理把前面这些硬骨头啃干净剩下的就是享受一套代码跑完八百多个站点的爽快感。本文还有配套的精品资源点击获取 日照时数转太阳辐射：840站点批量估算实战与避坑指南
简介面向全国840个气象站点的日照时数转日总太阳辐射Python代码专为气象科研人员、环境科学研究者及有一定Python基础的数据分析学习者设计。脚本覆盖数据读取、缺失值与异常值清洗、基于林格曼方程等经验公式的辐射换算、结果输出以及matplotlib可视化等完整流程并包含错误处理与性能优化可应用于气候变化研究、太阳能资源评估和气候模型构建等实际场景。压缩包共3个文件包括1个Python脚本和2个csv样例数据分别存放输入日照时数与转换后的输出结果便于对照验证整体仅275KB。已有6336人学习下载。该代码既可充当自动化处理840个站点日辐射数据的实用工具也可作为学习pandas数据操作、科学计算与可视化的实践范例帮助读者深入理解日照时数到日总太阳辐射的能量转换机制。 840个气象站点、几十年的日照时数记录要一口气转换成日总太阳辐射这个需求听起来就是把一个公式套进循环里跑一遍但真正动手时会发现数据不是攥在一个文件里的系数不是一个固定值通吃全国的算出来的数字也不是靠眼睛看就能判断靠不靠谱的。这篇文章我就把这套流程完整摊开讲包括数据预处理、天文辐射计算、Angström-Prescott模型的落地写法以及840个站点批量处理时我踩过的各种坑。适合做光伏资源评估、农业气候区划、建筑能耗模拟的朋友参考。1. 为什么要用日照时数来推太阳辐射1.1 直测辐射数据的缺口远比想象中大太阳总辐射的直接观测依赖日射强度计这类仪器本身不便宜而且对维护要求非常高石英罩要定期清洁、灵敏度会漂移、内部的干燥剂失效就得换一不留神数据质量就崩了。所以全国真正长期稳定发布总辐射观测的站点数量并不多要做省级网格、全国尺度的资源图谱光靠这些站点根本撑不起空间分辨率。相比之下日照时数的观测几乎每个国家级气象站都在做。日照计的记录相对皮实数据序列长、站点密度高全国840个站点的日照时数资料基本上可以把各个气候亚区都覆盖到。所以很多工程场景下日照时数成了太阳辐射的“代理变量”用经验模型把它换算成总辐射属于一种非常成熟的工程妥协。1.2 日照时数转辐射的经典模型日照时数和太阳辐射之间最经典的桥梁就是Angström-Prescott回归方程H H0 × (a b × n/N)其中H是目标日总太阳辐射H0是大气层顶的日天文辐射量也叫地外辐射n是当天实测日照时数N是当天理论可照时数也就是白昼时长a和b是经验系数。模型的核心逻辑很直观日照占比n/N越高云越少实际到达地面的辐射占天文辐射的比例就越高。这个模型从二十世纪二十年代提出至今经过无数区域修正依然在光伏、气象、农业领域大量使用。它最大的优点是只依赖日照时数一个常规观测要素就能计算而且只要标定好本地系数精度相当可观。业界也尝试过引入气温、湿度、降水做更复杂的机器学习模型但数据获取成本和可解释性都不如这个经典公式。840个站点要做统一批处理用这个模型是最合适的起点。2. 840个站点的原始数据先整理成能跑模型的样子2.1 数据文件长什么样读进来要注意什么国内数据集常见的形式是“一站一个文件”文件名就是区站号文件内是类似这样按空格或者固定宽度排列的文本57083 2015 01 01 0.0 57083 2015 01 02 8.2 57083 2015 01 03 -9999列的含义通常是区站号、年、月、日、日照时数单位小时。不同来源的数据缺失值编码不一样常见的有-9999、9999、32766这一类的占位符读进来之后第一步就是把它们替换成NaN。还有一个极其容易踩的细节有些数据集里日照时数会乘以0.1存储也就是说表里存的是8.2实际要除以10才是0.82小时不看元数据直接算结果会离谱到天际。先准备两份输入一份是站点元信息区站号、纬度、经度、海拔一份是长表观测记录。下面是我惯用的读法import pandas as pd import numpy as np # 站点信息 stations pd.read_csv( stations_840.csv, encodinggbk, dtype{station_id: str}, ) # 日照时数长表所有站点合并后 # 假设列名station_id, year, month, day, sunshine_hour obs pd.read_csv( sunshine_hour_all.csv, sepr\s, names[station_id, year, month, day, sunshine_hour], dtype{station_id: str}, ) obs[date] pd.to_datetime( obs[[year, month, day]] ) obs[sunshine_hour] obs[sunshine_hour].replace( [-9999, 9999, 32766], np.nan ) obs obs.dropna(subset[sunshine_hour]) obs obs[(obs[sunshine_hour] 0) (obs[sunshine_hour] 24)]最后这个范围过滤很关键日照时数不可能小于0也不应该超过24小时一些脏数据比如28.5、-3.2这时候就暴露了。2.2 合并站点信息和观测记录模型需要用到纬度和经度所以要先把站点元信息合并到观测长表里。合并之后还有一个动作不能省检查每个站点经纬度是否有缺失、是否有重复的日期记录。df obs.merge(stations, onstation_id, howleft) df df.dropna(subset[lat, lon]) # 检查同一个站点是否有重复日期 dupe df.duplicated(subset[station_id, date], keepFalse) if dupe.sum(): print(f发现 {dupe.sum()} 条重复日期记录请先人工核对) df df.drop_duplicates(subset[station_id, date])840个站点合并完行数通常是几百万级但这一步在pandas里也就是几秒的事。真正花时间的是检查那些“看起来成功但实际失败”的站点比如某些站点某一年文件缺失、某些站点经纬度为零值这些不查清楚后面算出来的辐射值就是一堆心里没底的数字。3. 辐射估算模型的核心计算一步一步拆给你看3.1 日地外辐射H0的计算Angström-Prescott里的H0不是一个常数它随着日期、纬度变化是天文位置决定的。计算H0需要三个量日地距离修正系数dr、太阳赤纬角δ、日落时角ωs。为了在840个站点上做向量化计算我建议全部用numpy批量处理不要写for循环逐步调用math库。先计算一年中的第几天然后算日角doy df[date].dt.dayofyear.to_numpy() # 日角弧度 gamma 2 * np.pi * (doy - 1) / 365.0 # 日地距离修正系数 dr 1 0.033 * np.cos(gamma) # 太阳赤纬角Spencer公式弧度 dec ( 0.006918 - 0.399912 * np.cos(gamma) 0.070257 * np.sin(gamma) - 0.006758 * np.cos(2 * gamma) 0.000907 * np.sin(2 * gamma) - 0.002697 * np.cos(3 * gamma) 0.00148 * np.sin(3 * gamma) )太阳赤纬角还可以用更简单的Cooper公式近似但Spencer公式的精度对日尺度计算来说更有保障尤其在春秋分附近两者差异能达到0.5度以上这对辐射值的影响不是可以忽略的。日地外辐射总量的标准公式是H0 (24×3600×Gsc/π) × dr × (ωs×sinφ×sinδ cosφ×cosδ×sinωs)其中Gsc是太阳常数工程计算里取1367 W/m²φ是纬度弧度。最终结果除以10^6转成MJ/m²/day。用numpy直接这样写lat_rad np.deg2rad(df[lat].to_numpy()) Gsc 1367.0 omega_0 np.arccos( np.clip(-np.tan(lat_rad) * np.tan(dec), -1.0, 1.0) ) H0 ( (24.0 * 3600.0 * Gsc / np.pi) * dr * ( omega_0 * np.sin(lat_rad) * np.sin(dec) np.cos(lat_rad) * np.cos(dec) * np.sin(omega_0) ) / 1e6 )这里有一个非常重要的细节-tan(lat)*tan(dec)在数学上可能落在[-1,1]区间之外高纬度地区尤其常见如果直接丢进arccos会返回NaN。用np.clip把输入夹到[-1,1]再计算能避免一整个站点因为冬季极端日期而全部报废。3.2 理论白昼时长N的计算白昼时长N本质上由日落时角直接换算卫星轨道上太阳中心在地平线上的持续时间单位是小时N (24/π) × ωs对应到代码就是daylight_hours 24.0 / np.pi * omega_0注意N的物理含义极昼情况下ωs等于π弧度N就是24小时极夜情况下ωs是0N就是0。由于前面已经做了clip这两个极端值在计算上天然成立不需要额外判断。但后续算n/N的时候N等于0的那一天会变成除零需要单独处理。3.3 a和b系数怎么取才不至于瞎算这是整个流程里最容易出争议的地方。Angström-Prescott方程的a、b系数并不是普适常数不同气候区、不同季节差异很大。全国840个站点横跨湿润、半干旱、干旱、高原等多个气候区用一个固定系数套全国误差在部分地区会非常难看。文献里的经典取值月尺度上很多研究用a0.18、b0.55这个组合在不少区域适用但并不全球通用。日尺度的估算通常推荐另一个量级。我做实际项目时习惯按气候区给几套参考系数然后再用靠近站点的辐射观测站做校验调整。适用区域ab备注华北平原、东北中部0.200.52春季风沙天误差偏大长江中下游0.180.50梅雨季容易低估四川盆地、云贵东部0.200.44多云地区日照-辐射关系弱青藏高原0.250.48建议再做海拔修正西北干旱区新疆、甘肃0.220.55晴天比例高相关性好如果你有站点附近辐射站的实测资料强烈建议做一次本地回归把自己地区的a、b直接标定出来。方法很简单把实测辐射除以H0得到y把日照占比n/N作为x做一元线性回归截距就是a斜率就是b。这套自标定流程在代码里也就是np.polyfit两行的事但效果远好于随便套文献值。4. 全国840个站点的批处理实现4.1 向量化计算告别840个for循环数据量在百万级时最直观的想法是for循环遍历每个站点然后内层再遍历每一天840个站点跑下来可能要等一阵子。实际上完全可以一步到位把所有站点所有日期的数据一次性丢进numpy数组做向量化计算几百万行也只需几秒到十几秒。把经纬度、日期、日照时数都合并到同一个DataFrame后按上面的方式算出H0和daylight_hours然后直接算估算辐射sunshine_hour df[sunshine_hour].to_numpy() # 日照占比注意极夜N0的情况 ratio np.divide( sunshine_hour, daylight_hours, outnp.zeros_like(sunshine_hour, dtypenp.float64), wheredaylight_hours 0, ) # 给每个站点挂上对应区域的a、b示例全国统一初值 a_val 0.20 b_val 0.52 ghi_estimate H0 * (a_val b_val * ratio) # 物理量约束辐射值不为负 ghi_estimate np.maximum(ghi_estimate, 0.0) df[ghi_estimate] ghi_estimatenp.divide里的where参数是个容易被忽略但是极其关键的操作它只允许白天时长大于0的日期参与除法剩下的地方直接用0填充省去一次mask处理。这样写出来的代码没有循环也没有额外的内存爆炸风险。4.2 输出结果的组织方式算完的结果建议不要只憋在一个超大CSV里后续按站点、按日期切片取数会很痛苦。我习惯按站点拆分输出一个站点一个文件或者整体存成parquet格式。下面这个写法实用且不容易出错import os os.makedirs(output/ghi_daily, exist_okTrue) for station_id, g in df.groupby(station_id): out g[[date, ghi_estimate]].sort_values(date) out.to_csv( foutput/ghi_daily/{station_id}.csv, indexFalse, float_format%.2f, )如果你后续想要直接做全国面板分析推荐用df.to_parquet(ghi_daily.parquet)一次性落盘Polars、Dask读取都比CSV快一个量级。我个人一般是两种都存CSV用来给业务同事看parquet留给自己做二次计算。5. 结果验算不要让看起来合理的数字骗了你5.1 用邻近辐射观测站做基准对比算完辐射值之后最忌直接拿去画图甚至写报告。840个站点里可能就有那么几十个站在数据质量上有隐藏问题必须做验证。最有效的办法是找到那些同时具备日照时数和实测总辐射的站点把模型估算值和实测值放在一起看偏差。抽样选十几个覆盖不同气候区的站点计算平均绝对误差MAE、均方根误差RMSE和相对偏差from sklearn.metrics import mean_absolute_error, mean_squared_error # est: 模型估算值, obs: 实测辐射值 mae mean_absolute_error(obs, est) rmse mean_squared_error(obs, est, squaredFalse) mape np.mean(np.abs((obs - est) / obs)) * 100 print(fMAE: {mae:.2f} MJ/m2/day) print(fRMSE: {rmse:.2f} MJ/m2/day) print(fMAPE: {mape:.1f}%)日尺度的Angström-Prescott模型误差到15%~25%都不算意外这是模型的正常波动除非你想做融资级别的精算否则不必太焦虑。但如果你发现某个站点MAPE超过了40%大概率不是模型问题而是这个站点的日照时数或辐射观测本身存在系统误差这时候要去查原始数据而不是硬调系数。5.2 通过物理自检快速淘汰异常在没有实测辐射站对照的区域至少要过一遍物理合理性检查。几个最实用的判断点日辐射估计值不应该超过同一天H0的值更不可能为负值。晴天日照时数接近N时日辐射估算值理论上应在H0的60%~75%区间如果出现接近H0甚至超过H0的数值说明系数设置或数据单位出了问题。同一站点相邻日期的辐射值不应出现剧烈震荡日照时数从0跳到12辐射值从2跳到25是合理的但日照时数都是10的连续两天辐射值却从20掉到5就要警惕。还可以按站点查看年总量均值。全国除四川盆地、贵州一带外年总辐射量多数在4500~6500 MJ/m²/year这个区间如果某站点算出来年总量只有2000或者冲到8000不要怀疑气候特殊先去检查数据。6. 我踩过的几个坑以及对应的解决办法6.1 arccos函数越界导致整年数据报废这是所有踩坑里影响面最大的一个。没加clip之前我把全国数据一次性丢进去算结果青藏高原和东北北部一大批冬季日期的H0变成了NaN一侧是坦坦荡荡的报错另一侧是几万个NaN在DataFrame里安静潜伏。如果后面直接groupby聚合这些NaN会自动被跳过你根本察觉不到直到画图时才发现一片空白。用np.clip把三角函数的输入限制在[-1,1]是唯一可靠的解法。数学上-tanφ·tanδ本来就该落在这个区间超出部分完全是由浮点误差和极端纬度造成的伪异常clip掉不会丢信息。6.2 极昼极夜地区除以零高纬度站点在冬至前后可能出现N0的日子此时不管实际日照时数是多少日照占比n/N都会变成inf或者NaN。用普通除法pandas会给出inf和警告后面再算辐射值就全部被污染。我习惯用np.divide的where参数处理原因前文提过它可以把N0的日子直接置为0占比对应的辐射计算结果也会被修正为合理小值。切记不能把这几天简单删掉否则站点年总量偏低的假象会误导后续分析。6.3 月份日界与日照计的“太阳日”差异日照时数观测在很多台站是以真太阳时为日界的而常规气象数据文件的日期是按北京时间日界归档的。这一点在大尺度统计上影响很小但在日出、日落时间极端的高纬地区冬季可能出现少数日期错位。表现为个别站点在日照时数为0的冬季日前后两天数值突变。处理思路不复杂如果你发现某站点冬季数据在日界附近有明显跳变可以尝试用三日滑动平均做平滑或者在月尺度上重新做聚合把这种边界噪声稀释掉。单独揪出某一天的精确对位成本太高对绝大多数工程需求没有必要。6.4 a、b系数“一招鲜”带来的区域失真早期我做全国批处理时图省事全量用同一组a、b系数结果西北地区估算值明显偏高四川盆地明显偏低。后来把全国分成几个气候区分别回归才把整体误差压下来。这提醒我一个道理模型的公式可以统一但参数绝不能出厂后就不动了。只要是做全国尺度的项目至少要按气候区或省级行政区给参数分档并且在成果里注明每个区域用的参数来源别让下游使用者在汇报时无据可查。这套代码跑通之后后续往里面加新站点只需要在stations表里增加经纬度和海拔重新跑一遍就出结果了。如果你需要进一步做月尺度、季尺度甚至年尺度分析对日结果按站点聚合即可pandas的groupby加resample一条链走下来顺滑得很。我实际做下来最大的体会是真正费时间的永远不是模型而是看不完的数据质量和没完没了的边界情况处理把前面这些硬骨头啃干净剩下的就是享受一套代码跑完八百多个站点的爽快感。本文还有配套的精品资源点击获取