简介一套面向气象科研与环境数据分析的Python处理脚本专门用于将全国840个气象站点的日照时数转换为日总太阳辐射适合需要批量处理辐射数据、开展太阳能评估或构建气候模型的研究者。压缩包仅275KB共3个文件其中py脚本承担核心转换流程两个csv文件分别提供输入站点的日照时数样例与转换后的辐射输出结果可直接替换数据运行。脚本实现覆盖数据读取、缺失值清洗、基于物理经验公式的辐射换算、结果导出及简单可视化并借助pandas提升处理效率便于扩展到更大规模站点数据。资源目前已有6343人学习下载对于想掌握气象数据自动化处理、理解日照与辐射换算原理的读者是一份可直接运行的实用模板。1. 日照时数转日总太阳辐射840 个站点批量算值在哪、坑在哪做光伏选址、作物模型或者生态水文模拟的从业者手头最不缺的就是全国气象站点逐日日照时数最缺的反而是太阳辐射数据。辐射观测站数量少、分布稀而日照站点密度高得多用日照时数反推日总太阳辐射是这个领域用了很多年的替代方案。原理一句话能讲完但真正把 840 个站点、逐日数据一次性跑完就会发现缺测编码、天文参数边界、区域系数选择每一环都能让结果偏掉两成以上。这套 Python 脚本做的事情就是把「日照时数 → 日总太阳辐射」的换算封装成可批处理的完整流程输入站点经纬度和逐日日照时数输出逐日总太阳辐射并保留天文辐射和可照时数等中间量供排查。适合正在做光伏资源评估、农业模型输入数据准备、径流与蒸散发模拟的人直接拿去改。2. 换算公式与天文参数为什么日总辐射能靠日照百分率算出来2.1 公式拆解地表辐射为什么和日照百分率近似线性日总太阳辐射的经典估算是日照百分率线性回归式Rs Ra × (a b × n/N)。Ra 是大气层顶天文辐射N 是当日可照时数n 是实测日照时数n/N 就是日照百分率。日照时数计记录的其实是当天直射辐射强度超过约 120 W/m² 的累计小时数它本身就是一个云量累计指标云多日照小时数短到达地面的辐射被削得厉害云少日照小时数接近可照时数辐射接近晴空水平。把透过率对日照百分率做线性近似a 就是 n/N0 即全天阴雨时的透过率ab 就是 n/N1 即全天晴朗时的透过率。默认的全国经验值通常取 a0.25、b0.50意味着最阴的天能透过约 25% 的天文辐射晴天能透过约 75%。这个线性关系在月尺度上很稳日尺度上会有噪声但逐日使用仍然是最常用的做法。2.2 三个天文参数赤纬、日地距离修正、日落时角要把 Ra 算准得先知道每天地球相对于太阳的几何位置。三个参数就够了太阳赤纬 δ 决定太阳直射点纬度日地距离修正因子 dr 修正轨道偏心导致的辐射通量差异日落时角 ωs 决定白天长短。这三者都只跟一年中的第几天有关写成函数非常短。import numpy as np def solar_declination(julian_day): 太阳赤纬单位弧度 return 0.409 * np.sin(2 * np.pi * julian_day / 365 - 1.39) def earth_sun_distance_factor(julian_day): 日地距离修正因子无量纲约 0.97~1.03 return 1 0.033 * np.cos(2 * np.pi * julian_day / 365) def sunset_hour_angle(lat_deg, declination): 日落时角单位弧度用 clip 兜住定义域边界 phi np.radians(lat_deg) arg -np.tan(phi) * np.tan(declination) return np.arccos(np.clip(arg, -1.0, 1.0))赤纬和日地距离因子采用的都是通用的近似展开式误差在 1% 以内对辐射估算完全够用。日落时角里那句np.clip是关键高纬度地区冬天 tan(φ)·tan(δ) 会超出 [-1, 1]不 clip 的话arccos直接返回 NaN后面全盘报销。2.3 天文辐射 Ra 与可照时数 N 的完整实现有了上面三个参数Ra 和 N 就是一步套用def extraterrestrial_radiation(julian_day, lat_deg): 大气层顶日总辐射单位 MJ/(m2·d) Gsc 0.0820 # 太阳常数单位 MJ/(m2·min) phi np.radians(lat_deg) dr earth_sun_distance_factor(julian_day) delta solar_declination(julian_day) ws sunset_hour_angle(lat_deg, delta) Ra (24 * 60 / np.pi) * Gsc * dr * ( ws * np.sin(phi) * np.sin(delta) np.cos(phi) * np.cos(delta) * np.sin(ws) ) return Ra def day_length(julian_day, lat_deg): 可照时数单位小时 delta solar_declination(julian_day) ws sunset_hour_angle(lat_deg, delta) return 24 * ws / np.piGsc 取 0.0820 MJ/(m²·min)这是太阳常数换算成每分钟能量通量的常用值前面的24 * 60 / np.pi把每分钟的量积成每日总量。Ra 输出单位是 MJ/m²如果后面要跟实测辐射对比注意对方用的单位。可照时数就是日落时角换算成小时一天最长不会超过 24高纬度夏季能到 20 小时以上。这套天文计算我建议直接抄别自己重新推导容易在系数上翻车。3. 脚本落地840 个站点的批处理流程与关键函数3.1 输入数据组织站点表加逐日观测表代码的输入按两张表组织。第一张是站点元数据包含站点编号、纬度、经度、海拔第二张是逐日观测包含日期、站点编号、日照时数。日照时数单位是小时正常范围 0 到 24。stn_id,lat,lon,elev 50136,53.47,122.35,433.0 50468,51.72,126.65,242.0date,stn_id,sunshine_hours 2020-01-01,50136,7.2 2020-01-02,50136,0.0 2020-01-03,50136,11.5拿到原始数据后第一件事不是算辐射而是清洗。气象站点的日照时数文件里常见 32744、32766、999999 这类特征缺测码直接参与计算会把结果拉飞。清洗函数要把缺测码、负值、超过 24 的值统一转成 NaNdef clean_sunshine(s): 把日照时数里的缺测编码和异常值统一转 NaN s s.replace([32700, 32744, 32766, 999998, 999999], np.nan) s s.mask(s 0, np.nan) s s.mask(s 24, np.nan) return s注意顺序先编码替换再做范围检查。有些站点会把微量负值记录成 -1 表示仪器未工作不能简单取绝对值直接转 NaN 更稳妥。3.2 单站换算从一个日期到一列结果单站的换算逻辑是批处理的基础先写成纯函数方便调试。输入日照时数、纬度、儒略日序号和系数输出当天的 Rsdef sunshine_to_daily_radiation(sunshine_hours, lat_deg, julian_day, a0.25, b0.50): 日照时数(小时) - 日总太阳辐射(MJ/m2) Ra extraterrestrial_radiation(julian_day, lat_deg) N day_length(julian_day, lat_deg) if N 0.1: return 0.0 # 极夜或可照时间极短日照百分率失去意义 ratio float(np.clip(sunshine_hours / N, 0.0, 1.0)) return (a b * ratio) * Ra这个函数把 n/N 强制限制在 0 到 1 之间这是必须的自动站日照计偶尔会把强散射光误判为直射导致 n 略大于 N如果日照百分率超过 1算出来的 Rs 可能超过晴空值物理上讲不通。极夜那天 N 接近于 0日照百分率没意义直接给 0.0 而不是报错。3.3 840 站批处理用向量化而不是循环套循环单站函数可以用于调试但 840 个站乘 365 天大约是 30 万行数据用apply逐行调 Python 函数会慢到让人怀疑人生。批处理要按站点分组把天文参数全部向量化计算def batch_sunshine_to_radiation(meta, obs, a0.25, b0.50): meta: 站点表; obs: 逐日日照表; 返回带 Rs/Ra/N 的结果表 frames [] for _, m in meta.iterrows(): sub obs[obs[stn_id] m[stn_id]].copy() if sub.empty: continue jd pd.to_datetime(sub[date]).dt.dayofyear.values.astype(float) phi np.radians(m[lat]) dr 1 0.033 * np.cos(2 * np.pi * jd / 365) delta 0.409 * np.sin(2 * np.pi * jd / 365 - 1.39) ws np.arccos(np.clip(-np.tan(phi) * np.tan(delta), -1.0, 1.0)) Ra (24 * 60 / np.pi) * 0.0820 * dr * ( ws * np.sin(phi) * np.sin(delta) np.cos(phi) * np.cos(delta) * np.sin(ws) ) N 24 * ws / np.pi n sub[sunshine_hours].values with np.errstate(divideignore, invalidignore): ratio np.where(N 0.1, n / N, 0.0) ratio np.clip(ratio, 0.0, 1.0) sub[Ra] Ra sub[N] N sub[Rs] (a b * ratio) * Ra frames.append(sub) return pd.concat(frames, ignore_indexTrue)外层按站点循环内层全部 numpy 数组运算30 万行在普通笔记本上几秒跑完。np.where(N 0.1, n / N, 0.0)负责极夜日的除零保护np.errstate抑制中间过程的警告但不影响结果。输出里同时保留 Ra 和 N是因为后面验证结果时还要用它们算晴空指数。4. 区域系数与边界条件高原、湿润区、高纬度站点的修正4.1 分区域系数参考全国统一 0.25/0.50 只是起点a、b 的经验值在不同气候区差异明显。全国默认 0.25/0.50 是温带中纬度地区的平均结果拿到高原或者南方湿润区直接用系统性偏差会很难看。常用来做分区域参考的一组取值如下区域类型ab适用场景全国默认0.250.50快速估算、无分区信息北方干旱半干旱0.22 ~ 0.250.50 ~ 0.55云量少、晴天透过率高南方湿润区0.18 ~ 0.220.55 ~ 0.62水汽多、阴天透过率低高原站点0.27 ~ 0.300.38 ~ 0.45海拔 3000 米以上、大气干洁物理逻辑是这样的湿润区水汽和气溶胶多阴天时云层更厚所以 a 偏小但偶尔放晴时空气通透ab 偏大。高原地区大气稀薄干洁散射比例低日照百分率和透过率的关系更陡所以 b 明显偏小、a 略大。如果没有本地率定条件至少按这个表给站点分组而不是全图一个系数。4.2 极夜与高纬度边界可照时数接近零的日子单独处理纬度高于大约 65 度的站点冬季会出现可照时数极短甚至全天无日照的日子。这时候 n/N 失去了物理意义分子分母都趋近于零任何微小的观测误差都会被放大成荒谬的百分比。批量函数里N 0.1这个阈值就是为这类日子准备的。N 小于 0.1 小时的天数结果直接置 0.0同时建议单独输出一个标记列方便后续统计时把这些点排除掉。不要试图在这些天强行计算日照百分率算出来的数字只会污染月均值。4.3 单位换算MJ/m² 与 kWh/m² 别记反光伏场景常见的是 kWh/m²而公式输出是 MJ/m²换算系数是 1 kWh/m² 3.6 MJ/m²。也就是说除以 3.6 就能从兆焦换成千瓦时sub[Rs_kwh] sub[Rs] / 3.6做光伏发电量估算时kWh/m² 在数值上就约等于峰值日照小时数可以直接乘上系统效率估算发电量。这个换算我在代码注释里每次都写清楚因为太容易在单位上栽跟头有人把 Ra 的 MJ 当成 kWh 直接用结果冬季辐射偏高 20% 以上查了半天才发现是单位问题。5. 常见问题与避坑日照资料里最容易翻车的五个细节5.1 算出来辐射比天文辐射还高现象某站点夏季晴天计算 Rs 接近甚至超过 Ra明显不物理。原因原始日照时数文件里有 32744、999999 这类缺测编码没被清洗被当成真实日照值参与计算另外自动站偶发记录 n 超过 N 的情况。解决清洗函数里先替换缺测编码再检查范围计算时强制np.clip(n/N, 0, 1)同时统计一下被 clip 的行数如果占比超过 1%回去查原始数据的质量标记。5.2 高纬度站点冬天整列 NaN现象程序在 12 月份的数据上崩掉或者输出整片 NaN而低纬度站点正常。原因tan(φ)·tan(δ) 超出 [-1, 1]arccos定义域越界极夜天数 N 等于 0除零产生 NaN。解决日落时角计算里加np.clip(arg, -1.0, 1.0)日照百分率用np.where(N 0.1, n/N, 0.0)而不是直接相除。这两处是整套代码里最容易被忽略但影响最大的边界保护。5.3 相邻两个站点结果系统性差 20%现象同一气候区、相距不到一百公里的两个站点月均太阳辐射差异明显不合理。原因站点表里纬度精度不一致有的站点经纬度只有小数点后一位有的保留到两位纬度差 0.5 度在冬季就会让 Ra 差出好几个百分点。另外可能有站点经纬度坐标被误填成度分格式但没有换算成十进制度。解决统一把站点经纬度转成十进制度保留至少五位小数批处理前先画一张站点分布图做目检看有没有站点位置明显偏离。5.4 辐射曲线的相位整体偏移现象夏季辐射偏低、冬季偏高或者峰值出现的时间比邻近站早一个月。原因日照时数的统计时次问题。气象站册子上的日照时数常常按前一日 20 时至当日 20 时统计和自然日错开了半天部分自动站又按北京时 0 点归日两种口径混在一起就会造成相位差。解决先读数据说明确认统计时次如果不一致把日期整体偏移一天再做dayofyear。这个坑最隐蔽因为不画逐日曲线根本看不出来。5.5 高原站点冬季系统性偏高现象全国统一用 0.25/0.50 跑完海拔 3000 米以上的站点冬季 Rs 比同纬度低海拔站点高 30% 以上。原因高原大气干洁、云量少日照百分率本就偏高加上 b 取 0.50 对高原而言偏大导致晴天透过率被高估。解决至少按第 4 章的分区域表给高原站点换系数如果手头有短期的实测辐射数据用最小二乘对Rs/Ra和n/N做回归率定本地 a、b。没有实测数据时宁可保守取 a0.28、b0.42也不要全国一刀切。6. 验证技巧用晴空指数和邻近站点给计算结果把脉6.1 晴空指数 Kt先给结果做一次快速体检Rs 算完不等于结束我习惯先算晴空指数 Kt Rs / Ra。Kt 的物理含义是大气总透过率正常范围应该在 0 到 1 之间。逐日情况下阴雨天 Kt 通常落在 0.1 到 0.3晴天在 0.6 到 0.8超过 0.85 基本可以断定不物理低于 0.05 多半是观测或清洗出了问题。result[kt] result[Rs] / result[Ra] bad result[(result[kt] 0.85) | (result[kt] 0.05)] print(fKt 异常占比: {len(bad) / len(result):.3%})正常数据的 Kt 异常占比应该远低于 1%。如果异常比例超过 2%先别急着分析结果回头查清洗和系数而不是继续往下走。这一条几乎每次都能拦住问题数据。6.2 邻近站点互证与月均值对比Kt 体检通过后再做两个交叉验证。第一个是按月统计各站 Rs 均值同一气候区的相邻站点月均值差异一般不超过 15%如果某站连续几个月明显偏离邻近站多半是站点经纬度或系数分组出了问题。第二个是把结果和气候学常识对一下全国大部分地区夏季月均日辐射在 15 到 25 MJ/m²冬季在 5 到 10 MJ/m²高原站夏季能到 25 以上但冬季不会低得离谱。明显越界的数值哪怕通过了 Kt 检查也要再确认一次站点表。monthly result.groupby( [result[date].dt.month, result[stn_id]] )[Rs].mean().unstack()从那以后我每次拿到一批新的日照数据都会先抽头 7 天跑一遍小样本把 Ra、N、Kt 三个中间量打出来看一遍确认曲线形态正常再全量执行最后才做月均对比。这套流程看起来多花十分钟实际上省掉的是后面反复排查数据的大半天。希望帮到你。本文还有配套的精品资源点击获取