840站点日照时数转日总太阳辐射的Python脚本实战

简介：面向气象、环境科学与太阳能资源评估方向研究者的Python工具，围绕全国840个气象站点日照时数（SSD）转日总太阳辐射（SR）这一常见数据处理任务整理。压缩包共3个文件，由2个CSV数据表和1个Python脚本组成，整体约275KB；CSV分别承载站点输入数据与计算结果，脚本负责读取、清洗、公式换算和批量输出，结构紧凑，便于套用到同类观测数据。已有6348人学习。脚本以经纬度、日期和日照时数为基础，结合林格曼方程或经验公式估算日总太阳辐射，并在计算前处理缺失值、异常值与日期格式不一致等问题，减少手工整理成本。对科研人员，它提供了可复核的批量计算模板，可替换站点数据与参数后用于气候变化、能源利用和气候模型的前期分析；对数据分析学习者，则可借pandas读写CSV、逐站计算与结果落盘的完整流程，理解气象指标换算和脚本组织方式，也方便在此基础上扩展统计汇总或可视化。

1. 840 个站点的日照时数摆在一起，先别急着套公式

前几年接手一个跨区域蒸散估算的活，手头只有一张 840 个气象站点的日值表：站号、日期、日照时数，三列。第一反应是 Ångström-Prescott，a 取 0.25、b 取 0.50，脚本半小时写完，跑完一看，高纬那几个站夏季算出来的日总太阳辐射比天文辐射还高；还有一批站的辐射与天文辐射之比全年是一条水平线。公式本身没错，问题出在输入变量 n 上：单位是 0.1 小时，没有除以 10；缺测码 32766 被当成了 3276.6 小时；可照时数 N 这一项干脆被漏掉了。这套 840 站点日照时数转日总太阳辐射的 python 脚本真正价值不在公式，而在于把单位、缺测、坐标、可照时数这几道坎提前处理掉了。凡是手上有气象站日值、要做蒸散、作物模型或者光伏资源初筛的，都可以照着往下走。

2. 日照时数怎么变成 MJ/m²：Ångström-Prescott 的四个输入量

Ångström-Prescott 关系本身只有一行：Rs/Ra = a + b·(n/N)。很多人就卡在这行式子里，以为把 n 代进去、乘一个系数就完了。实际上这个式子有四个输入量，n 只是其中之一，另外三个——天文辐射 Ra、可照时数 N、经验系数 a 和 b——每一个都能独立把结果带偏。下面按“先算天、再定系数、最后对齐数据”的顺序拆开。

2.1 天文辐射 Ra 与可照时数 N 的逐日算法

Ra 是大气层顶水平面上一天的累积辐射量，只跟纬度和年积日有关，跟云、气溶胶、海拔全都没关系。它之所以重要，是因为它给了 Rs 一个物理上限：地面日总辐射不可能显著超过 Ra。任何一站的 Rs/Ra 超过 1，一定是数据或代码出了问题，这就是后面排查时最有用的一条硬性限制。

计算链是这样的：先由年积日得到日地距离订正系数 dr 和太阳赤纬 δ，再由纬度和赤纬算出日落时角 ωs，最后把这三者代进 Ra 的积分式。可照时数 N 则是 ωs 换算成小时的结果，ωs 是弧度，乘 24/π 才是小时。这一步是第一个翻车点，很多人把 ωs 直接当小时用，N 就只有 1.5 左右，于是 n/N 动不动大于 1。

import numpy as np import pandas as pd GSC = 0.0820 # 太阳常数, MJ m-2 min-1 DAYS_IN_YEAR = 365.0 # 固定 365, 与通用手册口径保持一致 def ra_and_n(doy, lat_deg): """ doy : 年积日, 标量或 ndarray, 1-366 lat_deg : 纬度, 北纬为正, 十进制度 返回 : Ra (MJ m-2 d-1), N (h) """ doy = np.asarray(doy, dtype=float) phi = np.radians(np.asarray(lat_deg, dtype=float)) dr = 1.0 + 0.033 * np.cos(2 * np.pi * doy / DAYS_IN_YEAR) # 日地距离订正 delta = 0.409 * np.sin(2 * np.pi * doy / DAYS_IN_YEAR - 1.39) # 太阳赤纬, rad x = -np.tan(phi) * np.tan(delta) ws = np.arccos(np.clip(x, -1.0, 1.0)) # 日落时角, rad Ra = (24 * 60 / np.pi) * GSC * dr * ( ws * np.sin(phi) * np.sin(delta) + np.cos(phi) * np.cos(delta) * np.sin(ws) ) N = 24.0 / np.pi * ws # 可照时数, h return Ra, N

逐项说明一下。dr的振幅只有 ±3.3%，看着不起眼，但它决定了近日点和远日点的差异，在年尺度统计上会被放大。delta用的是通用近似式，误差在 0.1° 量级，对 Ra 的影响小于 0.2%，够用。np.clip(x, -1, 1)是防止 arccos 定义域越界的保底措施；中国陆域纬度在 3°N 到 54°N 之间，实际不会出现极昼极夜，但站点表里一旦混进纬度 90 或者符号写错的记录，np.tan会让 x 跑到几千，从而得到 NaN，而 NaN 会一路传到结果表里，很难追查。先把它限制在定义域内，再单独做一轮纬度范围检查，比事后查找 NaN 要省事。

参数方面，doy直接用dt.dayofyear，闰年自动是 366，但分母我固定用 365。这不是疏忽：通用手册的公式就是以 365 为分母，换成 366 会让你的结果和别人的结果对不上，差异在 0.1% 量级，没有意义。lat_deg必须是十进制度，北纬为正，南纬为负，如果你的站点表用度分整数存储，先别急着往里扔。

2.2 系数 a、b：默认值

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