企业官网建设流程全解析

简介：全国八百四十个气象站点日照时数转日总太阳辐射的python代码，面向气象科研人员、气候变化研究者以及数据分析学习者，用于解决将逐日日照时数批量换算为日总太阳辐射的重复计算问题，可服务于太阳能资源评估、气候模型构建和地表能量平衡研究。压缩包包含三个文件，两个表格数据文件和一个脚本文件，整体约二百七十五千字节；表格文件分别给出输入样例与转换结果，脚本文件则实现从数据读取、缺失值清洗、经验公式换算到结果保存的完整流程。脚本参考日照时数与总辐射的经验关系，结合站点经纬度、季节等参数完成计算，并借助pandas与matplotlib完成数据表处理和站点辐射分布图绘制，便于用户快速检查和比较不同区域结果。目前已有六千三百四十一人浏览学习，适合作为气象数据自动化处理的入门模板；非编程背景的研究者可直接替换数据运行，有编程基础的用户也能借此掌握批量气象数据的清洗、计算与可视化方法。

1. 用日照时数换太阳辐射：840个站点的批量解法

做光伏资源评估、作物模型模拟或者蒸散发计算时，日总太阳辐射是刚需输入，但国内辐射观测站点数量远比日照时数站点少。全国840个气象站都有日照时数记录，辐射站却只有百十来座，直接拿辐射数据做空间分析根本铺不开。日照时数转日总辐射的Angström-Prescott经验公式，是解决这个数据缺口最经典也最稳妥的路子，核心是用日照百分率拟合大气透过率，把日照时数换算成水平面日总辐射。这篇文章给一套能直接改的Python方案：从天文学公式到参数标定，再到840站批量处理与避坑，全程可以照着跑通。

2. 日照时数转辐射的核心原理：Angström-Prescott公式与两个天文参数

2.1 A-P公式：为什么只靠日照时数就能算辐射

日照时数本质上记录的是"太阳直接照射地面的累计小时数"，它和总辐射之间的关系，不是简单的线性倍增，而是通过大气对太阳辐射的衰减来联系的。Angström-Prescott公式把这个关系写成：

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

H是地面水平面日总辐射（MJ/m²），H0是大气层顶天文辐射日总量，S是实测日照时数，S0是可照时数（从日出到日落的理论时长），a和b是经验系数。S/S0叫日照百分率，a是全天阴天时大气透过率的近似，a+b是全天晴天的透过率近似，两者之差反映了云量对辐射的调制能力。

这个公式能成立的条件是：日照时数的累积效应和大气透过率之间存在统计稳定的关系。虽然真实大气里水汽、气溶胶、云类型都会影响辐射，但在日尺度上，日照相对时长已经携带了大部分信号。这也是为什么A-P公式从1924年提出至今，仍然是世界气象组织推荐的日照转辐射基础方法。相比需要输入云量、湿度、气压的物理模型，它需要的输入只有日照时数和站点纬度，非常适合大范围站点序列计算。

2.2 天文日照时数S0与天文辐射H0：两个必须先算对的量

S0和H0不是从观测数据里读出来的，是用天文公式算的。S0取决于纬度φ和太阳赤纬δ，赤纬又随年序日n变化。常用的Cooper公式足够精确：

δ = 23.45 × sin(360 × (284 + n) / 365)

然后算日落时角ωs：

cos(ωs) = -tan(φ) × tan(δ)

S0就是2 × ωs（度）除以15，单位小时。H0的计算要再引入日地距离修正因子E0和太阳常数Isc（取1367 W/m²）：

E0 = 1 + 0.033 × cos(2π × n / 365)

H0 = (86400 / π) × E0 × Isc × [ωs×sin(φ)×sin(δ) + cos(φ)×cos(δ)×sin(ωs)]

这里的ωs、φ、δ都要用弧度。算出来是J/m²，除以1e6转成MJ/m²。注意H0是大气层顶的辐射日总量，不是地面实测值，A-P公式里H0的作用是把辐射归一化，消除纬度、季节、日地距离的影响，让a和b在不同站点间有可比性。这也是为什么A-P公式能跨纬度使用的原因。

2.3 参数a、b：全国统一还是分站标定

a和b没有物理常数级别的标准值，它们是经验回归系数。同一个公式在不同气候区，a大约在0.18~0.25，b在0.45~0.55，a+b在0.65~0.75之间。干燥的高原地区a偏小、b偏大，因为晴天透过率高；湿润的南方地区a偏大、b偏小，因为阴天也有较多散射辐射。直接从文献抄一套全国统一参数是最快的方法，但误差通常有10%~20%，做精细项目时不够用。

更可靠的做法是用辐射观测站的数据标定：某个辐射站如果有同步的实测辐射H_obs和日照时数S，就先用天文公式算该站的H0和S0，然后对H_obs/H0和S/S0做一元线性回归，截距就是a，斜率就是b。这个做法在时间序列上要求至少有一整年的逐日数据，覆盖晴、阴、雨各种天气状况，得到的参数才稳定。全国840个站里只有少数站有同步辐射观测，标定后可以把参数用空间插值推广到全部站点，这个进阶流程放到第6章讲。

3. 数据准备与参数标定：先把手里的气象数据整理成能算的样子

3.1 站点元数据与日照数据的标准格式

批量处理的第一步，是把散乱的原始数据规整成统一格式。我一般会准备两个文件：一个是站点元数据表，包含站点ID、纬度（十进制度，北纬为正）、经度、海拔；另一个是日照时数数据表，长表格式，每行包含站点ID、日期（YYYYMMDD或YYYY-MM-DD）、日照时数值。原始数据从观测网拿下来通常是按年分文件的文本格式，站点ID占列、日期占列、值可能带奇怪的小数位和缺测标记，需要先做清洗。

缺测值和异常值的标记在各数据源里不统一，常见的有32766、9999、空白字符串。清洗时要统一转成NaN，否则后面计算会带着脏数据跑。单位也要盯紧，有的数据集把日照时数放大了10倍存储（即记录值除以10才是小时），有的则直接是小时，建议在清洗后打一个简单的数据概览：按站点统计最大值、小于0的个数、大于24的个数，一眼看出单位问题。

3.2 用辐射站点标定a、b的回归代码

标定流程是整条链路里唯一需要"数据质量洁癖"的环节。先选一个有同步辐射和日照观测的站点，要求至少连续一年逐日数据，剔除缺测和明显异常（辐射大于天文辐射、日照为负等），然后做线性回归。以下是标定代码的完整实现，注释里标了每个步骤的用途。

import pandas as pd import numpy as np def solar_geometry(n, lat_deg): """计算年序日n对应赤纬、日落时角、可照时数和天文辐射日总量。 输入n为1-366，lat_deg为纬度(度,北纬为正)。""" # Cooper公式计算赤纬，单位先转弧度 decl_deg = 23.45 * np.sin(np.radians(360 * (284 + n) / 365)) decl_rad = np.radians(decl_deg) lat_rad = np.radians(lat_deg) # 日落时角余弦，注意acos定义域要截断到[-1,1] cos_omega = -np.tan(lat_rad) * np.tan(decl_rad) cos_omega = np.clip(cos_omega, -1, 1) omega_rad = np.arccos(cos_omega) # 可照时数，单位小时 S0 = 2 * np.degrees(omega_rad) / 15 # 日地距离修正因子E0 E0 = 1 + 0.033 * np.cos(2 * np.pi * n / 365) # 天文辐射日总量，Isc=1367 W/m2，结果从J/m2转MJ/m2 Isc = 1367 H0 = (86400 / np.pi) * E0 * Isc * ( omega_rad * np.sin(lat_rad) * np.sin(decl_rad) + np.cos(lat_rad) * np.cos(decl_rad) * np.sin(omega_rad) ) / 1e6 return S0, H0 # 读辐射站数据：列含date, sunshine_h, radiation_MJ df = pd.read_csv("rad_station.csv", parse_dates=["date"]) df["n"] = df["date"].dt.dayofyear # 计算该站逐日天文参数 lat = df["lat"].iloc[0] # 站点纬度，假设表里有 S0_list, H0_list = [], [] for n in df["n"]: s0, h0 = solar_geometry(n, lat) S0_list.append(s0) H0_list.append(h0) df["S0"] = S0_list df["H0"] = H0_list # 剔除辐射超天文辐射的物理不可能数据 df = df[df["radiation_MJ"] <= df["H0"] * 1.02] df = df.dropna(subset=["sunshine_h", "radiation_MJ"]) # 构造回归变量 df["y"] = df["radiation_MJ"] / df["H0"] # 归一化辐射 df["x"] = df["sunshine_h"] / df["S0"] # 日照百分率 # 期望截距a在0.1~0.3, 斜率b在0.4~0.6, a+b约0.7 a, b = np.polyfit(df["x"], df["y"], 1) print(f"标定结果: a={a:.3f}, b={b:.3f}, a+b={a+b:.3f}")

这段代码的逻辑：先逐日算天文参数，再构造归一化变量做一元回归。np.polyfit返回的第二个值是斜率b，第一个是截距a。标定结果做合理性检查：a小于0、b大于0.8或者a+b超出0.6~0.8，都要回头查数据，大概率是辐射和日照记录不同步或者单位出了问题。回归样本量最好大于300天，少于一整年得到的参数在季节外推时不可靠。

3.3 没有辐射站时用文献分区参数

全国840个站不可能每个都有同步辐射观测，实际操作中是先标定有辐射资料的站，再向周边无辐射站推广。如果项目精度要求不高，或者只是做宏观趋势分析，直接用分区参考值也能跑。常见的做法是按气候区给参数：湿润区a取0.22~0.24、b取0.46~0.50，半干旱区a取0.20~0.23、b取0.48~0.53，干旱区a取0.18~0.21、b取0.52~0.57。具体数值可以查所在区域已有文献的报道，但要注意文献里H0的计算方式是否和你一致。

有个隐藏的坑：不同文献对H0的定义不一样，有的用天文辐射日总量，有的用晴天辐射（即大气上界辐射经过理想大气衰减后的值）。用晴天辐射做基准的A-P公式，a、b数值和用天文辐射做基准的不能混用。我见过有人把文献里的b=0.62直接套进天文辐射基准的公式，结果辐射估算值普遍偏高。判断基准的方法很简单：如果a+b接近0.75~0.8，大概率是基于晴天辐射；基于天文辐射的a+b通常不超过0.75。

4. 840个站点批量转换的Python实现：按站循环与向量化结合

4.1 主流程设计：读取元数据、逐站计算、统一输出

批量处理的核心是流程组织。840个站点的日照数据如果分散在多个文件里，按站循环逐文件处理最直观，但文件IO开销大；如果所有站点数据在一个大表里，用pandas的groupby一次算完更高效。我建议先按数据组织方式选方案：日照数据按"站点×日期"的宽表存储时，用groupby；按站分开存储时，用循环处理。无论哪种，都要先把站点纬度合并进日照数据表，因为计算S0和H0必须知道纬度。

主流程分三步：加载站点元数据、合并纬度、逐站（或分组）计算并追加结果。840站的日序列如果每站有20年数据，总行数约600万，pandas处理没有任何压力，不需要上数据库或者分布式。计算完成后统一输出长表CSV，列包括站点ID、日期、日照时数、S0、H0、估算辐射，后续做统计和GIS分析都方便。

4.2 核心函数：一次性算出S0、H0和日辐射

单个站点的日辐射估算函数可以写成独立模块，后续无论是循环还是groupby都复用同一套逻辑。下面是完整函数，跟前一版的区别是增加了DataFrame批量处理和异常状态标记。

import pandas as pd import numpy as np def calc_daily_radiation(df, lat_deg, a, b): """输入df包含date和sunshine_h两列, lat_deg为站点纬度。 返回原df复制并新增S0, H0, radiation_MJ三列。""" out = df.copy() n = out["date"].dt.dayofyear.to_numpy() lat_rad = np.radians(lat_deg) # 赤纬(弧度) decl_rad = np.radians(23.45 * np.sin(np.radians(360 * (284 + n) / 365))) # 日落时角(弧度), 截断保证acos有定义 cos_omega = -np.tan(lat_rad) * np.tan(decl_rad) cos_omega = np.clip(cos_omega, -1, 1) omega_rad = np.arccos(cos_omega) # 可照时数(小时) out["S0"] = 2 * np.degrees(omega_rad) / 15 # 日地距离修正因子和天文辐射(MJ/m2) E0 = 1 + 0.033 * np.cos(2 * np.pi * n / 365) out["H0"] = (86400 / np.pi) * E0 * 1367 * ( omega_rad * np.sin(lat_rad) * np.sin(decl_rad) + np.cos(lat_rad) * np.cos(decl_rad) * np.sin(omega_rad) ) / 1e6 # A-P公式估算地面日总辐射 sunshine = out["sunshine_h"].to_numpy() valid_mask = (sunshine >= 0) & (out["S0"] > 0) out["radiation_MJ"] = np.nan out.loc[valid_mask, "radiation_MJ"] = ( out.loc[valid_mask, "H0"] * (a + b * sunshine[valid_mask] / out.loc[valid_mask, "S0"]) ) return out

这段代码用numpy数组操作代替逐日循环，840站全量计算时速度快一个量级。S0等于0的极夜情况必须单独处理，否则除以0会出现inf，然后污染整列数据。valid_mask里把S0>0作为条件，极夜天的辐射直接留NaN，后续统计时用skipna忽略即可。这里a和b作为参数传入，方便替换成逐站标定值或插值结果。

4.3 批量跑完840个站：输出结果的组织方式

实际批量运行我推荐用一个配置字典管理所有路径和参数，避免脚本里到处是硬编码路径。下面是主脚本的核心片段，采用按站循环方式，适合数据文件按站存储的场景。

import os import pandas as pd from pathlib import Path # 配置项：路径与参数 CONFIG = { "meta": "data/station_meta.csv", # 站点元数据: id, lat, lon, elev "sun_dir": "data/sunshine/", # 按站存储的日照文件, 文件名含站号 "out": "output/radiation_daily.csv", # 输出文件 "a": 0.22, "b": 0.50, # 默认参数, 可用标定值替换 } meta = pd.read_csv(CONFIG["meta"]) results = [] for row in meta.itertuples(): fpath = Path(CONFIG["sun_dir"]) / f"sun_{row.id}.csv" if not fpath.exists(): print(f"缺失: {row.id}") continue sun = pd.read_csv(fpath, parse_dates=["date"]) sun = sun.dropna(subset=["sunshine_h"]) # 逐站计算, 传入该站纬度和统一参数 rad = calc_daily_radiation(sun, row.lat, CONFIG["a"], CONFIG["b"]) rad["station_id"] = row.id results.append(rad) # 合并并输出, 按站和日期排序 all_rad = pd.concat(results, ignore_index=True) all_rad = all_rad.sort_values(["station_id", "date"]) all_rad.to_csv(CONFIG["out"], index=False, float_format="%.3f") print(f"完成, 共{len(all_rad)}行, 输出至{CONFIG['out']}")

循环里加了一个文件存在性检查，某个站文件缺失时不中断整体流程，而是打印提示后继续。这在处理真实批量数据时很实用，840个站里总有那么一两个文件损坏或者命名不规范，一个try-except包住单站处理逻辑更稳。输出用float_format限制小数位，文件体积更小，下游工具读起来也更快。

5. 避坑：日照数据里最容易翻车的5个细节

5.1 现象：日照时数大于可照时数

某站点冬季某天日照记录为9.2小时，但可照时数S0只有8.7小时，导致S/S0大于1，辐射估算值偏高。原因通常是仪器过灵敏或者日出日落时刻的弱辐照被误判，也可能是数据清洗时把前一天的尾数串了行。解决：在计算辐射前加一个截断，把日照时数超过S0的值压到S0。气象观测规范里日照时数确实允许略超理论可照时数（大气折射效应），但超出量一般不超过0.1~0.2小时，超太多就是数据问题。代码里加一行sunshine = np.minimum(sunshine, S0)即可。

5.2 现象：辐射估算结果整体偏高或偏低

同样的参数，同一批代码，A地区的估算辐射比实测平均高25%，B地区却正常。这类现象十有八九是参数标定的基准不匹配，也就是第3章提到的H0定义问题。解决：重新核对所用a、b是从哪种H0基准下拟合出来的。如果是基于晴天辐射的，要么把参数换算回天文辐射基准，要么把公式里的H0换成你的晴天辐射模型。这个坑隐蔽在参数文献里，很多论文不会明确写基准，只能通过a+b的数值反推。

5.3 现象：高纬度站点冬季出现NaN或inf

北纬50度以上的站在冬至前后，tan(φ)×tan(δ)绝对值大于1，acos输入越界，程序直接报错或者产生NaN。原因就是极夜期间不存在日出日落，公式本身失效。解决：用np.clip(cos_omega, -1, 1)把输入截断，S0算出来是0或接近0，再通过valid_mask把辐射置为NaN。不要试图用负数时角或者负日照去"修正"，物理上极夜没有日照也没有辐射，NaN是正确结果。

5.4 现象：缺测值被当成0值参与计算

日照时数文件里缺测标记是32766，清洗时只做了dropna(subset=["sunshine_h"])，但pandas读进来时32766是整数，不会自动变成NaN。于是这天的日照被当成3276.6小时参与计算，辐射估算值巨大。原因是对数据源的缺测编码不熟悉。解决：读取数据时加参数na_values=[32766, 9999, -9999]，或者在清洗阶段先做值域检查，日照时数超过24的直接置NaN。push到正式流程前，先对每个站做describe()看最大值和缺失数。

5.5 现象：a+b的合理性检查被忽略

标定回归得到的a=0.1、b=0.75，a+b=0.85，远超常见范围0.6~0.8，但回归R²很高，让人误以为结果可靠。原因可能是H0计算错误，比如赤纬公式里角度和弧度混用，或者太阳常数取错。解决：建立硬性校验规则，a小于0或大于0.4、b小于0.3或大于0.7、a+b小于0.6或大于0.8的标定结果直接标记不通过，返回检查天文计算函数。这个规则不绝对，但能拦住绝大部分低级错误。

6. 进阶：把结果从"能算"做到"可信"：交叉验证与参数空间化

批量跑通只是起点，真正要对外出数还需要做两件事：验证精度、做参数空间化。验证用留一站交叉验证：全国有辐射观测的站点里，每次留出一个站不参与标定，用其余站标定的参数（或插值参数）去估算留出站的辐射，然后和实测值对比。统计指标用MBE（平均偏差，单位MJ/m²）、RMSE（均方根误差）和R²。我一般要求日尺度RMSE控制在2.5~3.5 MJ/m²以内，月尺度RMSE小于1.5 MJ/m²，超过这个范围就说明参数推广的假设有问题。

参数空间化是提升840站整体精度的关键手段。具体做法是：把标定好的各辐射站点a、b值连同站点经纬度做成散点，用反距离权重插值或普通克里金插值到全部840个站点所在位置，使每个站都有自己的a、b，而不是共用一套全国参数。插值前先检查标定值的空间自相关，如果东南西北乱跳，说明标定样本太少或数据质量问题，插值结果不可靠。插值后做平滑处理，并限制a、b的取值范围，避免在数据稀疏区出现极端组合。

进阶方向还包括把日总量拆分到小时。光伏模拟需要小时辐射时，可以基于日照时数的日变化过程假设，用Collares-Pereira模型把日总辐射拆成逐时辐射，只需要输入日总量和可照时数。这个拆分的精度依赖云的时间分布假设，雨季误差偏大，但用于长期平均或气候态分析足够。

最后说个我自己的教训：最早做这套流程时，我用全国统一参数跑完840站，对比某省实测辐射站数据，发现夏季系统性偏高。排查了一周，最后定位到是日照数据的单位问题——手头数据里日照时数实际是0.1h记录，我当成小时直接用了。从那以后，任何新数据源我都会先画三张图：日照时数的月均值曲线、辐射估算的月均值曲线、和实测辐射的散点图，图比任何统计指标都更能暴露问题。数据驱动的活，每一步都得留个心眼。希望帮到你。

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