企业官网建设流程全解析

简介：面向气象、环境与新能源领域科研人员的Python实用工具，用于将全国840个气象站点的日照时数批量转换为日总太阳辐射，解决太阳能资源评估与气候研究中的基础数据换算问题。压缩包共3个文件，包括1个py脚本和2个csv数据文件，总体积约275KB；py脚本覆盖数据读取、清洗、公式转换、结果保存等环节，csv文件提供输入样例与转换后的输出结果，便于对照验证。已有6336人学习/下载。脚本基于常见经验公式，并结合站点经纬度与季节因素，使用者只需准备含日期、经纬度、日照时数的表格数据即可运行，也可参考代码改造为更大规模站点的处理流程。对初学者而言，这是一份能直接上手的自动化气象数据处理模板；对有一定编程基础的用户，可通过阅读代码掌握pandas数据操作、批量循环计算、异常值处理及结果导出等实用技巧，有助于深化对地表能量平衡和太阳能资源评估的理解。 做水文模型、生态模型、光伏资源评估的人，十有八九都遇到过这个尴尬：模型核心输入要日总太阳辐射（MJ/m²·d），但手里翻来覆去只有日照时数。地面辐射观测站全国就100多个，而日照时数观测站有840个，而且序列普遍长得多。用日照时数反推太阳辐射，业内最经典的就是Angström-Prescott公式（中文圈常叫A-P公式或埃斯特龙公式），原理不复杂，代码量也不大。这篇就把全国840个气象站点批量转换的完整思路和Python代码拆开讲清楚，从天文辐射计算、可照时数推导，到a、b系数的拟合和批量跑批避坑，一次性说透。适合做气候分析、农业模型、太阳能资源评估的研究生和工程师参考，小白也能照着跑通。

1. 思路与原理：为什么日照时数能推算总辐射

1.1 Angström-Prescott公式的基本形态

先理清楚一个概念：到达地面的太阳总辐射，由直接辐射和散射辐射组成，日照时数本身只是“太阳直射持续的时间长短”，它跟辐射能量之间并不是线性关系，但存在强相关性。1924年Angström提出用日照百分率估算辐射，后来Prescott改成地外辐射做基准，变成现在最常用的形式：

H / H0 = a + b × (S / S0)

拆开看：

• H：实际地表日总太阳辐射（MJ/m²·d），这就是我们要算的目标。
• H0：大气层顶的天文辐射日总量（MJ/m²·d），也叫地外辐射，理论上可以精确计算。
• S：实测日照时数（小时），来自气象站观测数据。
• S0：可照时数（小时），即理论白昼长度，同样可以计算。
• a、b：经验系数，反映大气对辐射的削弱程度，跟地区云量、气溶胶、水汽有关。

这个公式把比值关系拆成两部分：先算“天上该有多少辐射”，再用日照百分率去折算“实际到达地面多少”。所以代码的核心就落到两个计算上：一个是天文辐射H0，一个是可照时数S0，最后套系数。

1.2 天文参数与可照时数的计算要点

天文辐射取决于三个因素：日地距离、太阳赤纬、纬度。计算过程绕不开几个关键参数：

• 日序J：1月1日为1，12月31日为365（闰年366），代码里要把日期转成这个序列。
• 日地距离修正dr：地球绕太阳轨道是椭圆的，每一天实际接收的辐射有正负3.3%的波动。
• 太阳赤纬δ：太阳直射点纬度随季节在南北纬23.5°之间摆动。
• 时角ωs：决定日出日落时刻，也直接决定可照时数。

计算式如下（单位统一说明，角度全部用弧度制）：

dr = 1 + 0.033 × cos(2πJ / 365) δ = 0.409 × sin(2πJ / 365 − 1.39) ωs = arccos(−tan(纬度) × tan(δ))

求到ωs后，可照时数S0 = 24 × ωs / π（单位小时）；天文辐射日总量H0 = (24×3600×Gsc/π) × dr × (ωs×sinφ×sinδ + cosφ×cosδ×sinωs) / 10⁶，其中Gsc为太阳常数1367 W/m²，φ为纬度（弧度制）。

这里最容易栽跟头的是弧度换算和arccos溢出。纬度超过66.5°有极昼极夜，tan值相乘可能超过1，arccos直接返回NaN。国内站点基本在53°N以南，但保险起见还是要做截断，否则批量跑批的时候遇到边界站点就断掉了。

2. 数据准备：站点信息与日照数据的组织方式

2.1 两份基础数据怎么整理

处理840个站点，首先要把输入数据分成两个文件，建议都用CSV：

站点信息表（station_info.csv），必含四列：

• station_id：站点号，字符串类型，别用数字，不然像“58027”这种带前导零的会被pandas读成数字丢信息。
• lat：纬度，单位度，北纬为正。
• lon：经度，单位度（算地方时修正用得上，但A-P公式本身只用纬度）。
• alt：海拔高度，备用。

日照数据表（sunshine_hours.csv），长表格式：

• date：日期，统一成YYYY-MM-DD字符串。
• station_id：站点号，跟站点表一一对应。
• sunshine_hours：当日日照时数，单位小时。

长表是pandas最擅长的格式，别用“一个站点一列”的宽表，后面groupby和merge都不顺手。实测下来，气象数据站导出Excel常有合并单元格、标题行错位的问题，清洗阶段建议直接用pd.read_csv(encoding='gbk')，国内气象数据GBK编码很常见。

2.2 日期转日序的技巧

日期转日序J，最直接的是datetime模块的timetuple().tm_yday，但840个站点、多年逐日数据动辄几十万行，用Python原生的datetime逐行循环会慢得让人怀疑人生。正确姿势是用pandas向量化：

import pandas as pd sunshine['date'] = pd.to_datetime(sunshine['date']) sunshine['year'] = sunshine['date'].dt.year sunshine['doy'] = sunshine['date'].dt.dayofyear # 自动处理闰年

dayofyear拿到的就是1到365/366的日序，pandas底层是C实现的，速度比for循环快一到两个数量级。就算数据量到百万级，也就一两秒的事。注意闰年会出现日序366，代码里cos公式用的365不受影响，可以放心用。

3. 核心代码实现：从单站函数到全国批量跑批

3.1 先写一个单站的辐射计算函数

我不建议一上来就写840站的循环，先把单个站点的函数调通，再向量化批量，出错好定位。下面这个函数就是整套代码的核心计算单元：

import numpy as np import pandas as pd def calc_daily_radiation(lat_deg, doy, sunshine_hours, a=0.18, b=0.55): """ 计算日总太阳辐射(MJ/m2/d) lat_deg: 纬度, 单位度 doy: 日序, 1-366 sunshine_hours: 当日日照时数, 小时 a, b: Angstrom-Prescott经验系数 """ phi = np.radians(lat_deg) # 纬度转弧度 Gsc = 1367 # 太阳常数 W/m2 # 日地距离修正 dr = 1 + 0.033 * np.cos(2 * np.pi * doy / 365) # 太阳赤纬 delta = 0.409 * np.sin(2 * np.pi * doy / 365 - 1.39) # 时角，注意截断避免arccos溢出 cos_ws = -np.tan(phi) * np.tan(delta) cos_ws = np.clip(cos_ws, -1.0, 1.0) ws = np.arccos(cos_ws) # 可照时数(小时) S0 = 24 * ws / np.pi # 天文辐射日总量(MJ/m2/d) H0 = (24 * 3600 * Gsc / np.pi) * dr * ( ws * np.sin(phi) * np.sin(delta) + np.cos(phi) * np.cos(delta) * np.sin(ws) ) / 1e6 # 日照百分率, 避免日照时数超过可照时数 ratio = np.clip(sunshine_hours / S0, 0.0, 1.0) # 实际太阳辐射 H = H0 * (a + b * ratio) return H

这个函数有几处值得说。一是np.clip处理arccos溢出，等于给极昼极夜情况兜了底；二是日照百分率也做了截断，因为实际观测中偶尔会出现日照时数略大于理论可照时数（仪器灵敏度、太阳视半径等综合原因），不截断就会算出超过物理上限的离谱值；三是一切运算都放在numpy里，方便后面向量化。

3.2 有实测辐射数据时，怎么拟合a、b系数

很多地区的站点是有太阳辐射实测值的，那就别拍脑袋取系数，直接用最小二乘拟合。做法是把公式变形成一元线性回归：

y = H / H0，x = S / S0，y = a + b × x

写起来很简单：

from scipy import stats # 假设有实测辐射数据 df = pd.DataFrame({ 'H_obs': [12.5, 15.1, 18.3, ...], # 实测日总辐射 'H0': [30.2, 32.8, 35.1, ...], # 对应日期的天文辐射 'S': [6.2, 8.5, 10.1, ...], # 实测日照时数 'S0': [11.3, 11.9, 12.5, ...] # 对应日期的可照时数 }) df['y'] = df['H_obs'] / df['H0'] df['x'] = df['S'] / df['S0'] slope, intercept, r_value, p_value, std_err = stats.linregress(df['x'], df['y']) b_fit = slope a_fit = intercept print(f"拟合结果: a={a_fit:.3f}, b={b_fit:.3f}, R2={r_value**2:.3f}")

拟合的时候要注意：用月尺度数据拟合和用日尺度数据拟合，系数会有差别，日尺度拟合的R²普遍低一些，这是正常的。站点如果四季分明，还建议按月或按季节分组拟合，能得到12组a、b系数，精度会明显提升。

如果没有实测辐射数据，就看站点所在气候分区，参考邻近有辐射观测站点的系数。全国尺度上，常用的默认值a≈0.18、b≈0.55，但这是个粗糙的均值，西北干旱区b值能到0.60以上，湿润地区a值偏大、b值偏小，直接套默认值做研究会被审稿人挑战。

3.3 全国840站批量跑批的完整代码

单站函数调通后，批量处理就是数据合并和分组应用的事。强烈建议用numpy的向量化而不是一个站点一个站点for循环——840个站点，每天一条数据，10年就是300多万行，for循环分组长了也会慢。

# 读入数据 station_info = pd.read_csv('station_info.csv', encoding='gbk') sunshine = pd.read_csv('sunshine_hours.csv', encoding='gbk') # 合并站点纬度 df = pd.merge(sunshine, station_info[['station_id', 'lat']], on='station_id', how='left') # 日期转日序 df['date'] = pd.to_datetime(df['date']) df['doy'] = df['date'].dt.dayofyear # 计算可照时数和天文辐射 phi = np.radians(df['lat'].values) doy = df['doy'].values dr = 1 + 0.033 * np.cos(2 * np.pi * doy / 365) delta = 0.409 * np.sin(2 * np.pi * doy / 365 - 1.39) cos_ws = np.clip(-np.tan(phi) * np.tan(delta), -1.0, 1.0) ws = np.arccos(cos_ws) S0 = 24 * ws / np.pi H0 = (24 * 3600 * 1367 / np.pi) * dr * ( ws * np.sin(phi) * np.sin(delta) + np.cos(phi) * np.cos(delta) * np.sin(ws) ) / 1e6 # 按站点给不同的a、b系数: 用字典映射 ab_params = {'58027': (0.17, 0.56), '58362': (0.18, 0.55)} # 示例，实际按站点给 a_arr = df['station_id'].map(lambda x: ab_params.get(x, (0.18, 0.55))[0]).values b_arr = df['station_id'].map(lambda x: ab_params.get(x, (0.18, 0.55))[1]).values ratio = np.clip(df['sunshine_hours'].values / S0, 0.0, 1.0) df['solar_radiation'] = H0 * (a_arr + b_arr * ratio) # 输出结果 df[['station_id', 'date', 'solar_radiation']].to_csv('daily_solar_radiation.csv', index=False, encoding='utf-8')

这里有个细节容易被忽略：merge以后如果某些站点在station_info里缺失，lat会变成NaN，tan(NaN)会静默算出NaN，最后结果整行变空。所以merge之后一定要跑一句检查：

assert df['lat'].notna().all(), "存在站点未匹配到纬度信息，请检查station_id"

跑批之前先做数据体检，永远比跑完再排查省时间。

4. 实战中的坑与排查经验

4.1 高频问题速查表

第一个arccos报错我再多说一点。国内有个别站点在高海拔接近40°N的位置，虽然不容易触发极圈条件，但纬度精确到小数点后两位时，tan(phi)×tan(delta)在秋分前后可能非常接近1，浮点误差就足以让arccos超界。所以np.clip这行我从来不省。

4.2 几条实操心得

第一，别把a、b系数当常数。全国840个站用一个a=0.18、b=0.55，属于“能用但不够严谨”的层次。我处理全国数据时，通常先按气候区或纬度带分组：35°N以上、35°N到25°N、25°N以南，各拟合一套系数，再对个别有实测辐射的站点单独拟合，效果比统一系数好不少。

第二，注意数据单位。气象站的日照时数单位是小时，但有的时候一些数据库导出来会给分钟，比如把“7.5小时”记成“450”。这个错误特别隐蔽，因为数值一直为正，算出来结果也不会报错，就是整体偏大。批量跑之前先抽查几十条记录，跟站点天气现象记录做个对拍，几秒钟的事能避免返工。

第三，如果要做月尺度或者年尺度分析，别把日值算完再平均，直接先算H0和S0的月均值，再套A-P公式。因为公式里H0和S0是日变量，非线性叠加会引入偏差，虽然不大，但作为科研数据会被细究。

第四，多站点批量输出的时候建议带上太阳辐射的flag标识，比如把S/S0小于0.1的日值标记为“云量过大，估算值仅供参考”。总辐射估算本身在中高纬度阴雨天的误差最大，做应用的时候心里有数。

我在实际处理这批840站数据时，还发现一个细节：站点纬度如果是手工录入的，偶尔会有整度数和整十分度数的差别，有的是度分秒格式混入，读进来以后要统一做一次check，否则个别站点算出来S0跟当地日出日落时间对不上。用一个已知站点的任意一天倒推验证一下：比如北京1月15日，可照时数应该在9小时45分钟左右，算出来偏差超过10分钟就该回头查纬度了。

最后再分享一个小技巧：把840站的转换结果按年和月做一次空间可视化，把异常值直接画到地图上看，比如冬天华南站却算出了高于西北的辐射，那就是纬度合并配错了。眼睛看图找错，常常比在表格里刷数字更快。

整套代码跑下来，几百万条数据几分钟出结果。日照时数转总辐射这个需求，在气象、农业、能源领域会一直存在，A-P公式虽然不是精度最高的方法，但它胜在数据需求低、物理意义清晰、可复现性强，作为批量生产的工程方案非常靠谱。参数标定那步真的值得多花点功夫，系数准了，后面所有应用都省心。

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