日照时数转日总太阳辐射：840站批量换算的工程实现与避坑指南

简介：一套面向气象科研与环境数据分析的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/N=0 即全天阴雨时的透过率，a+b 就是 n/N=1 即全天晴朗时的透过率。默认的全国经验值通常取 a=0.25、b=0.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.pi

Gsc 取 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.0

date,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 的值统一转成 NaN：

def 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 单站换算：从一个日期到一列结果

单站的换算逻辑是批处理的基础，先写成纯函数方便调试。输入日照时数、纬度、儒略日序号和系数，输出当天的 Rs：

def sunshine_to_daily_radiation(sunshine_hours, lat_deg, julian_day, a=0.25, b=0.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, a=0.25, b=0.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(divide='ignore', invalid='ignore'): 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_index=True)

外层按站点循环，内层全部 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 是温带中纬度地区的平均结果，拿到高原或者南方湿润区直接用，系统性偏差会很难看。常用来做分区域参考的一组取值如下：

物理逻辑是这样的：湿润区水汽和气溶胶多，阴天时云层更厚，所以 a 偏小；但偶尔放晴时空气通透，a+b 偏大。高原地区大气稀薄干洁，散射比例低，日照百分率和透过率的关系更陡，所以 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。没有实测数据时，宁可保守取 a=0.28、b=0.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(f"Kt 异常占比: {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 三个中间量打出来看一遍，确认曲线形态正常再全量执行，最后才做月均对比。这套流程看起来多花十分钟，实际上省掉的是后面反复排查数据的大半天。希望帮到你。

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