简介一份面向科研人员与数据分析爱好者的Python代码包可将全国840个气象站点的日照时数SSD转换为日总太阳辐射SR适用于气候变化研究、太阳能资源评估与气候模型构建等场景。资源包为rar压缩包共3个文件包含1个Python脚本和2个CSV数据文件站点输入数据与辐射输出结果整体大小275KB轻量易用。脚本覆盖数据读取、预处理、基于林格曼方程或经验公式的转换计算、结果保存及可视化等完整流程并包含错误处理与优化。已有6343人学习下载适合需要批量处理气象站点辐射数据的科研人员也适合想学习pandas数据处理与自动化脚本的Python爱好者。通过这套代码读者可直接套用模板完成日照时数到日总太阳辐射的批量转换也可借此理解大规模气象数据的清洗、计算与输出思路。1. 840 个气象站点日照时数转日总太阳辐射一份可以直接运行的 Python 脚本某天我拿到一张全国 840 个气象站点的日照时数表字段只有站点编号、经纬度、日期和日照小时数任务是把它们转成日总太阳辐射。当时的直觉是“找个公式套进去不就行了”真正动手才发现坑比想象中多纬度要转弧度、日序不能按 365 天拍脑袋、天文辐射公式里那串三角函数写错一个符号输出就全是 NaN。这篇文章就把这套 Python 脚本的完整实现拆开讲清楚——算法怎么选、代码怎么写、840 个站点批量跑起来怎么不出错、拿到结果怎么验证。适合正在做太阳能资源评估、农业气象或水文模型的开发者这也是一个可以直接拿去跑、按自己数据改字段名就能用的脚本。2. 算法与数据结构从日照时数到天文辐射再到日总辐射的换算链路2.1 840 个站点数据集的真实形态先明确我们处理的是什么。全国 840 个气象站点的日照时数数据常见的落盘格式是长表 CSV一行代表一个站点某一天的观测记录典型字段是station_id、longitude、latitude、date、sunshine_hours。有些源数据会用宽表格式——每个站点一列、日期在行上这种先melt成长表再处理。sunshine_hours的物理含义要特别确认它是一天中太阳直接辐照度超过 120 W/m² 的累计时长单位通常是小时。有些数据集会把分钟数当小时数存如果不做单位校验后面算出来的日总辐射会整体偏大这类问题很难从单一数值上发现必须做整站总量复核。提示拿到数据第一步永远是df.describe()加df.isna().sum()先了解每一列的缺失情况和量纲再谈计算。2.2 日照时数转日总辐射的四步换算整个转换并不复杂核心是经验公式 Angstrom-Prescott日总太阳辐射 Rs 与太阳天文辐射 Ra 成线性关系斜率由日照百分率决定。Rs (as bs * n / N) * Ra其中n是实际日照时数N是当天最大可能日照时数as、bs是经验系数。这公式看起来简单麻烦的是N和Ra都要从天文参数推导由日序求太阳赤纬 delta进而求日落时角 ws由 ws 算出最大可能日照时数 N由日地距离修正因子和 ws 算出天文辐射 Ra最后套 Angstrom-Prescott 公式得到日总辐射。这四步对应的物理量必须一一算对。很多教程只给最终公式把太阳赤纬、时角这些中间量直接跳过导致新手把纬度、角度混着用算出的 Ra 数值离谱。下面逐个参数拆开讲。2.3 关键天文参数与 Angstrom-Prescott 系数的经验取值日序doy直接用date.dt.dayofyear获取注意闰年。太阳赤纬用 Cooper 近似公式角度制下是23.45 * sin(2π * (284 doy) / 365)这里极易踩坑的是 sin 内部的表达式用角度单位还是弧度单位。日地距离修正因子 dr 用1 0.033 * cos(2π * doy / 365)。日落时角 ws 是arccos(-tan(phi) * tan(delta))其中 phi 必须转成弧度。天文辐射 Ra 的单位是 MJ/m²/day这是后续所有计算的基准单位。太阳常数 Gsc 取值 0.0820 MJ/m²/min公式为24 * 60 / π * Gsc * dr * (ws * sin(phi) * sin(delta) cos(phi) * cos(delta) * sin(ws))。最大可能日照时数 N 用24 / π * ws换算成小时。as 和 bs 的经验值水文模拟里常用 as0.25、bs0.50两者之和接近 0.75代表晴空条件下大气透射率上限部分区域会用 as0.18、bs0.55适用于高原清洁大气。脚本里把 as、bs 作为命令行参数暴露换区域不必改代码。参数含义典型值 / 单位Gsc太阳常数0.0820 MJ/m²/mindr日地距离修正因子0.97~1.03delta太阳赤纬-23.45°~23.45°ws日落时角0~πN最大可能日照时数0~24 hasAngstrom 截距0.25bsAngstrom 斜率0.503. 核心代码拆解天文辐射 Ra 与日照百分率批量计算的实现3.1 读取数据与字段归一化不管源数据列名是什么第一步统一映射。我用一个RENAME_MAP字典把不同命名风格的字段归一到固定名称这样后续计算函数不用关心原始数据长什么样。import numpy as np import pandas as pd RENAME_MAP { station: station_id, site_id: station_id, lat: latitude, lon: longitude, 经度: longitude, 纬度: latitude, 日照时数: sunshine_hours, } def load_data(path): df pd.read_csv(path, encodingutf-8-sig) df.columns [RENAME_MAP.get(col, col) for col in df.columns] required [station_id, longitude, latitude, date, sunshine_hours] missing [col for col in required if col not in df.columns] if missing: raise ValueError(f缺乏必需字段: {missing}) df[date] pd.to_datetime(df[date]) df[sunshine_hours] pd.to_numeric(df[sunshine_hours], errorscoerce) return df逻辑说明RENAME_MAP负责把中英文混用的表头统一成内部字段名pd.to_datetime把日期字符串解析成 datetime 类型解析失败会自动变 NaT后面再统一处理。pd.to_numeric加errorscoerce很关键——不这么写一个脏字符就能让整列变成 object 类型后续运算直接报错。参数说明path是 CSV 路径返回的 DataFrame 保证含五个必需字段。3.2 一次性算出日序、赤纬、日落时角接下来用 numpy 向量化计算天文参数。840 个站点 × 365 天数据量在 30 万行级别用 for 循环逐行算在 Python 里是灾难numpy 列运算才是正解。def calc_solar_params(df): doy df[date].dt.dayofyear.to_numpy() lat_rad np.radians(df[latitude].to_numpy()) delta np.radians( 23.45 * np.sin(np.radians(360 * (284 doy) / 365)) ) dr 1 0.033 * np.cos(2 * np.pi * doy / 365) ws np.arccos( np.clip(-np.tan(lat_rad) * np.tan(delta), -1.0, 1.0) ) return doy, lat_rad, delta, dr, ws逻辑说明doy由dayofyear直接提取numpy 数组参与后续运算。delta的公式在外层转了np.radians保证sin和cos收到的参数是弧度如果漏掉这一步赤纬会被当成弧度参与计算结果偏差可达 20% 以上。ws的np.clip是对 arccos 定义域的兜底——高纬度地区冬季-tan(phi)*tan(delta)可能略小于 -1不 clip 的话会输出 NaN处理 840 个站点时这种边缘值一定会遇到。参数说明返回的doy范围是 1~366lat_rad单位弧度delta单位弧度ws单位弧度。这些中间量全部保持弧度制直到最后一步输出才转换单位避免混用。3.3 天文辐射 Ra 与日照百分率批量计算天文辐射 Ra 的计算是核心。公式24*60/π * Gsc * dr * (ws * sin(phi) * sin(delta) cos(phi) * cos(delta) * sin(ws))本质是对整个白天太阳辐照度的时间积分所以24*60是把每分钟的辐照度累加成一整天。GSC 0.0820 # MJ/m2/min def calc_radiation(df): doy, lat_rad, delta, dr, ws calc_solar_params(df) ra ( 24 * 60 / np.pi * GSC * dr * ( ws * np.sin(lat_rad) * np.sin(delta) np.cos(lat_rad) * np.cos(delta) * np.sin(ws) ) ) max_sunshine 24 / np.pi * ws sunshine_frac np.clip( df[sunshine_hours].to_numpy() / max_sunshine, 0.0, 1.0 ) rs (0.25 0.50 * sunshine_frac) * ra return pd.DataFrame({ station_id: df[station_id], date: df[date], longitude: df[longitude], latitude: df[latitude], sunshine_hours: df[sunshine_hours], day_total_radiation_MJ: np.round(rs, 3), day_total_radiation_kWh: np.round(rs / 3.6, 4), })逻辑说明ra的计算整列向量化840 个站点 30 万行数据在一秒内完成。max_sunshine是当天理论日照上限用24/π * ws把弧度折算成小时。sunshine_frac是实际日照时数与上限的比值clip 到 0~1 是因为仪器误差会导致个别记录超过理论值。最终rs (0.25 0.50 * sunshine_frac) * ra晴天sunshine_frac1时 Rs 等于 0.75 倍 Ra符合晴空辐射经验上限。参数说明输出同时给 MJ/m²/day 和 kWh/m²/day 两列光伏设计常用后者省得二次换算。rs / 3.6是 MJ 到 kWh 的换算关系因为 1 kWh 3.6 MJ。3.4 结果导出保留原始字段与编码处理输出 CSV 时把原始站点信息一并保留方便后续按站点分组统计。编码用utf-8-sig是为了让 Excel 打开时中文表头不乱码这是经常被忽略但实际很影响交付体验的细节。def main(input_path, output_path): df load_data(input_path) result calc_radiation(df) result.to_csv(output_path, indexFalse, encodingutf-8-sig) print(f处理完成: {len(result)} 条记录) print(f日总辐射范围: {result[day_total_radiation_MJ].min():.2f} f~ {result[day_total_radiation_MJ].max():.2f} MJ/m2/day)逻辑说明main函数把数据加载、辐射计算、导出串成完整流程print输出的最小值和最大辐射值是最快的 sanity check——如果最大值超过 40 MJ/m²/day大概率是纬度或日序计算出了问题这个阈值在中国区域基本不会误报。4. 运行与调试流程命令行、输出文件与结果验证4.1 环境准备与最小依赖这个脚本只依赖 pandas 和 numpy没有其他第三方库Python 3.9 以上直接跑。虚拟环境安装命令python -m venv .venv source .venv/bin/activate # Windows 下是 .venv\Scripts\activate pip install pandas numpypandas 负责读取和输出numpy 负责向量化数学运算。如果你的源数据是 Excel 格式再加一个pip install openpyxl就行读取时pd.read_excel替代pd.read_csv。整个环境装下来不到 200MB和动不动就装全套地理库的方案比这份脚本在普通办公笔记本上就能跑。4.2 命令行参数设计脚本入口支持四个关键参数列名映射在代码里改日常运行不需要动代码。python sunshine_to_radiation.py \ --input 全国840站点_日照时数_2023.csv \ --date-col date \ --lat-col latitude \ --sun-col sunshine_hours \ --output out/daily_radiation_2023.csv逻辑说明--input指向原始 CSV 路径--date-col、--lat-col、--sun-col指定源文件里的日期、纬度、日照时数列名适配不同来源的数据集--output是结果路径。这样做的好处是换一批数据只需要改命令行参数不需要打开脚本改代码——对于经常换数据源的气象数据处理场景来说这个设计能省下大量重复劳动。参数字段按argparse实现缺省值分别是date、latitude、sunshine_hours通常不用传。4.3 输出文件长什么样结果 CSV 每行是一个站点一天的数据列结构如下字段类型说明station_idstr站点编号datedatetime日期longitude / latitudefloat经纬度sunshine_hoursfloat原始日照时数day_total_radiation_MJfloat日总辐射MJ/m²/dayday_total_radiation_kWhfloat日总辐射kWh/m²/dayday_total_radiation_MJ是算法直接输出第二个 kWh 列是换算值。用 Excel 打开后先看最大值和最小值——中国区域内日总辐射的合理范围大致是 0~35 MJ/m²/day如果出现 50 甚至 100 的数字说明公式里某个环节出了问题。4.4 结果验证三步走拿到输出别急着用先做三件事验证。第一步单站抽检。挑一个纬度居中的站点把这一站全年的日总辐射求和得到年总量单位 MJ/m²应该在 4000~7000 之间。低于 3000 说明日照百分率被严重低估可能原始日照时数单位是分钟高于 8000 则要查天文辐射计算。第二步季节趋势检查。北半球站点 6 月的日均辐射应显著高于 12 月如果出现冬季反而高的倒挂现象大概率是日序公式里加减天数出了错。第三步与公开参考数据对照。选 3~5 个站点把多年平均的月总辐射和某再分析数据集或气候公报里的数值比对误差在 10% 以内算正常。气候态的 Angstrom 系数本身有区域差异超过 15% 的偏差优先怀疑系数选取是否适合本地。5. 避坑指南与常见问题排查单位、角度、缺失值与高原站点5.1 数据读取类问题NaN 与全零输出现象脚本跑完没有报错但输出文件里的辐射值全是 NaN或者某一站全为 0。原因最常见的是日期字段解析失败pd.to_datetime遇到2023/1/1 00:00这类半小时粒度时间戳时会产生 NaTdayofyear随之变成 NaN整行乘出来全是 NaN。另一个常见原因是日照时数列存在---或9999这类缺失值标记pd.to_numeric如果没加errorscoerce整个列会变成字符串类型运算结果全部是 NaN。解决load_data里强制pd.to_numeric(errorscoerce)并在读取后打印缺失统计。对于“全为 0”的情况直接检查原始 CSV 里该站点是否确实有观测部分站点冬季日照时数记录为 0 是正常的不一定有错。可以在输出前加一行print(df.groupby(station_id)[sunshine_hours].agg([count, min, max]))这个打印能快速发现哪些站点缺数据、哪些站点数值范围异常。5.2 公式计算类问题辐射值反常的大现象某几个站点的日总辐射超过 50 MJ/m²/day明显超出物理上限。原因纬度没有转弧度。np.tan(lat_rad)和np.sin(lat_rad)用的必须是弧度如果直接用df[latitude]的度数参与运算在纬度 30°~50° 时误差能达到 30% 以上。另一种可能是太阳赤纬公式里外层少了np.radiansdelta变成了一组很小的数值天文辐射 Ra 虚高。解决把calc_solar_params里所有三角函数输入统一打印出来检查量级。sin(delta)应该在 -0.4~0.4 之间ws应该在 0~3.2 之间如果ws出现复数或大于 π 的值基本就是角度和弧度混用了。我在代码里会把lat_rad、delta、ws三个量各打印一行 min/max跑完扫一眼心理就有底。5.3 高纬度站点的边界问题现象纬度高于 66° 的站点在冬季出现 NaN夏季辐射值却比邻近站点高出一截。原因日落时角公式ws arccos(-tan(phi) * tan(delta))在极昼极夜区域会出现定义域越界——冬天-tan(phi)*tan(delta)可能小于 -1arccos直接返回 NaN夏天它又可能大于 1。没有 clip 的版本会静默产出空值不会报错这最坑人。解决np.clip(-np.tan(lat_rad) * np.tan(delta), -1.0, 1.0)这行代码就是在兜这个底。极夜期间日照时数为 0辐射本身应该是 0clip 后 arccos 返回 πN 和 Ra 算出来是一个合理的小值再乘以 0 的日照百分率得到 0逻辑自洽。5.4 云南贵州等低日照区域系统性偏低现象某几个站点的年总辐射比气候公报低 20% 以上但不是所有站点都这样。原因Angstrom 系数 as0.25、bs0.50 是在特定区域的多年平均拟合结果对四川盆地、贵州这些多云地区并不适用。这些区域日照时数普遍偏低但实际散射辐射占比高总辐射未必像日照时数显示的那么低。经验公式在这种气候区先天吃亏不是代码 bug。解决把 as 和 bs 暴露成命令行参数--as-coef 0.26 --bs-coef 0.42用本地站的多年实测辐射率定系数。如果没有实测数据至少换用高原地带常用的 as0.18、bs0.55 试跑一遍对比差异。这类调整属于参数校准问题不该在代码里硬编码。5.5 840 个站点跑得极慢的问题现象脚本跑了几小时才出结果CPU 占用还不到 30%。原因用 for 循环逐行遍历 DataFrame 计算三角函数。30 万行数据在 Python 的纯循环里要跑 20 分钟以上而同样的计算用 numpy 向量化一秒钟内完成。解决脚本里calc_solar_params和calc_radiation全部用 numpy 数组整列运算没有任何一层 Python 循环。如果你拿到别的版本写的循环实现重构时只需要记住一条原则——所有涉及三角函数、幂运算、乘除的逐行计算一律先.to_numpy()再整列操作最后拼回 DataFrame。从apply改成向量化840 个站点一年的数据从 20 分钟压缩到 3 秒这个提升立竿见影。6. 二次开发思路把批量脚本改造成可配置的辐射计算工具项目源码只是个能跑通的最小实现实际开发里 840 个站点的数据往往是按月分文件存放还可能有几千个文件需要批量处理需要考虑把脚本封装成可扩展的工具形态。我的建议是把核心函数重构成一个类把 as/bs 系数和单位换算都作为实例属性这样二次开发时只需实例化不同配置就能服务多种数据源class RadiationConverter: def __init__(self, as_coef0.25, bs_coef0.50, output_unitMJ): self.as_coef as_coef self.bs_coef bs_coef self.output_unit output_unit def convert_file(self, input_path, output_path, date_coldate, lat_collatitude, sun_colsunshine_hours): df load_data(input_path) result calc_radiation(df, self.as_coef, self.bs_coef) if self.output_unit kWh: result[day_total_radiation] result[day_total_radiation_MJ] / 3.6 result.to_csv(output_path, indexFalse, encodingutf-8-sig)这个封装把硬编码的 0.25 和 0.50 变成了构造参数把年份、区域、数据源差异全部隔离在实例配置里。再往下可以加一个config.yaml描述各数据源的列名映射用yaml.safe_load读取后自动适配不同机构提供的数据省得每次拿到新数据都改代码。另一个值得投入的改造方向是增加质量诊断模式。现在脚本只管算不管数据质量实际跑气象数据时缺测、异常值、站点漂移才是真正的血泪来源。加一个--validate开关输出每个站点的有效记录数、日照时数年总量、平均日辐射跑完数据回头看一眼就能定位问题这比算完再检查要高效得多。我还要加一个功能多文件批量处理。840 个站点如果按站点分文件存放外层用一个for循环遍历目录每个文件计算完追加写出避免一次性加载到内存里卡死。最后再加个--skip-errors参数某个文件格式有问题时跳过并记录日志不让整个批次中断——这在实际生产中比任何算法优化都更实用。从那以后我每次拿到气象数据都会先跑一遍describe()和isna().sum()确认单位再开始算算完立刻看最大最小值范围不再急着把结果交给下游。这套流程救过我好几次希望帮到你。本文还有配套的精品资源点击获取 840个气象站点日照时数转日总太阳辐射的Python脚本
简介一份面向科研人员与数据分析爱好者的Python代码包可将全国840个气象站点的日照时数SSD转换为日总太阳辐射SR适用于气候变化研究、太阳能资源评估与气候模型构建等场景。资源包为rar压缩包共3个文件包含1个Python脚本和2个CSV数据文件站点输入数据与辐射输出结果整体大小275KB轻量易用。脚本覆盖数据读取、预处理、基于林格曼方程或经验公式的转换计算、结果保存及可视化等完整流程并包含错误处理与优化。已有6343人学习下载适合需要批量处理气象站点辐射数据的科研人员也适合想学习pandas数据处理与自动化脚本的Python爱好者。通过这套代码读者可直接套用模板完成日照时数到日总太阳辐射的批量转换也可借此理解大规模气象数据的清洗、计算与输出思路。1. 840 个气象站点日照时数转日总太阳辐射一份可以直接运行的 Python 脚本某天我拿到一张全国 840 个气象站点的日照时数表字段只有站点编号、经纬度、日期和日照小时数任务是把它们转成日总太阳辐射。当时的直觉是“找个公式套进去不就行了”真正动手才发现坑比想象中多纬度要转弧度、日序不能按 365 天拍脑袋、天文辐射公式里那串三角函数写错一个符号输出就全是 NaN。这篇文章就把这套 Python 脚本的完整实现拆开讲清楚——算法怎么选、代码怎么写、840 个站点批量跑起来怎么不出错、拿到结果怎么验证。适合正在做太阳能资源评估、农业气象或水文模型的开发者这也是一个可以直接拿去跑、按自己数据改字段名就能用的脚本。2. 算法与数据结构从日照时数到天文辐射再到日总辐射的换算链路2.1 840 个站点数据集的真实形态先明确我们处理的是什么。全国 840 个气象站点的日照时数数据常见的落盘格式是长表 CSV一行代表一个站点某一天的观测记录典型字段是station_id、longitude、latitude、date、sunshine_hours。有些源数据会用宽表格式——每个站点一列、日期在行上这种先melt成长表再处理。sunshine_hours的物理含义要特别确认它是一天中太阳直接辐照度超过 120 W/m² 的累计时长单位通常是小时。有些数据集会把分钟数当小时数存如果不做单位校验后面算出来的日总辐射会整体偏大这类问题很难从单一数值上发现必须做整站总量复核。提示拿到数据第一步永远是df.describe()加df.isna().sum()先了解每一列的缺失情况和量纲再谈计算。2.2 日照时数转日总辐射的四步换算整个转换并不复杂核心是经验公式 Angstrom-Prescott日总太阳辐射 Rs 与太阳天文辐射 Ra 成线性关系斜率由日照百分率决定。Rs (as bs * n / N) * Ra其中n是实际日照时数N是当天最大可能日照时数as、bs是经验系数。这公式看起来简单麻烦的是N和Ra都要从天文参数推导由日序求太阳赤纬 delta进而求日落时角 ws由 ws 算出最大可能日照时数 N由日地距离修正因子和 ws 算出天文辐射 Ra最后套 Angstrom-Prescott 公式得到日总辐射。这四步对应的物理量必须一一算对。很多教程只给最终公式把太阳赤纬、时角这些中间量直接跳过导致新手把纬度、角度混着用算出的 Ra 数值离谱。下面逐个参数拆开讲。2.3 关键天文参数与 Angstrom-Prescott 系数的经验取值日序doy直接用date.dt.dayofyear获取注意闰年。太阳赤纬用 Cooper 近似公式角度制下是23.45 * sin(2π * (284 doy) / 365)这里极易踩坑的是 sin 内部的表达式用角度单位还是弧度单位。日地距离修正因子 dr 用1 0.033 * cos(2π * doy / 365)。日落时角 ws 是arccos(-tan(phi) * tan(delta))其中 phi 必须转成弧度。天文辐射 Ra 的单位是 MJ/m²/day这是后续所有计算的基准单位。太阳常数 Gsc 取值 0.0820 MJ/m²/min公式为24 * 60 / π * Gsc * dr * (ws * sin(phi) * sin(delta) cos(phi) * cos(delta) * sin(ws))。最大可能日照时数 N 用24 / π * ws换算成小时。as 和 bs 的经验值水文模拟里常用 as0.25、bs0.50两者之和接近 0.75代表晴空条件下大气透射率上限部分区域会用 as0.18、bs0.55适用于高原清洁大气。脚本里把 as、bs 作为命令行参数暴露换区域不必改代码。参数含义典型值 / 单位Gsc太阳常数0.0820 MJ/m²/mindr日地距离修正因子0.97~1.03delta太阳赤纬-23.45°~23.45°ws日落时角0~πN最大可能日照时数0~24 hasAngstrom 截距0.25bsAngstrom 斜率0.503. 核心代码拆解天文辐射 Ra 与日照百分率批量计算的实现3.1 读取数据与字段归一化不管源数据列名是什么第一步统一映射。我用一个RENAME_MAP字典把不同命名风格的字段归一到固定名称这样后续计算函数不用关心原始数据长什么样。import numpy as np import pandas as pd RENAME_MAP { station: station_id, site_id: station_id, lat: latitude, lon: longitude, 经度: longitude, 纬度: latitude, 日照时数: sunshine_hours, } def load_data(path): df pd.read_csv(path, encodingutf-8-sig) df.columns [RENAME_MAP.get(col, col) for col in df.columns] required [station_id, longitude, latitude, date, sunshine_hours] missing [col for col in required if col not in df.columns] if missing: raise ValueError(f缺乏必需字段: {missing}) df[date] pd.to_datetime(df[date]) df[sunshine_hours] pd.to_numeric(df[sunshine_hours], errorscoerce) return df逻辑说明RENAME_MAP负责把中英文混用的表头统一成内部字段名pd.to_datetime把日期字符串解析成 datetime 类型解析失败会自动变 NaT后面再统一处理。pd.to_numeric加errorscoerce很关键——不这么写一个脏字符就能让整列变成 object 类型后续运算直接报错。参数说明path是 CSV 路径返回的 DataFrame 保证含五个必需字段。3.2 一次性算出日序、赤纬、日落时角接下来用 numpy 向量化计算天文参数。840 个站点 × 365 天数据量在 30 万行级别用 for 循环逐行算在 Python 里是灾难numpy 列运算才是正解。def calc_solar_params(df): doy df[date].dt.dayofyear.to_numpy() lat_rad np.radians(df[latitude].to_numpy()) delta np.radians( 23.45 * np.sin(np.radians(360 * (284 doy) / 365)) ) dr 1 0.033 * np.cos(2 * np.pi * doy / 365) ws np.arccos( np.clip(-np.tan(lat_rad) * np.tan(delta), -1.0, 1.0) ) return doy, lat_rad, delta, dr, ws逻辑说明doy由dayofyear直接提取numpy 数组参与后续运算。delta的公式在外层转了np.radians保证sin和cos收到的参数是弧度如果漏掉这一步赤纬会被当成弧度参与计算结果偏差可达 20% 以上。ws的np.clip是对 arccos 定义域的兜底——高纬度地区冬季-tan(phi)*tan(delta)可能略小于 -1不 clip 的话会输出 NaN处理 840 个站点时这种边缘值一定会遇到。参数说明返回的doy范围是 1~366lat_rad单位弧度delta单位弧度ws单位弧度。这些中间量全部保持弧度制直到最后一步输出才转换单位避免混用。3.3 天文辐射 Ra 与日照百分率批量计算天文辐射 Ra 的计算是核心。公式24*60/π * Gsc * dr * (ws * sin(phi) * sin(delta) cos(phi) * cos(delta) * sin(ws))本质是对整个白天太阳辐照度的时间积分所以24*60是把每分钟的辐照度累加成一整天。GSC 0.0820 # MJ/m2/min def calc_radiation(df): doy, lat_rad, delta, dr, ws calc_solar_params(df) ra ( 24 * 60 / np.pi * GSC * dr * ( ws * np.sin(lat_rad) * np.sin(delta) np.cos(lat_rad) * np.cos(delta) * np.sin(ws) ) ) max_sunshine 24 / np.pi * ws sunshine_frac np.clip( df[sunshine_hours].to_numpy() / max_sunshine, 0.0, 1.0 ) rs (0.25 0.50 * sunshine_frac) * ra return pd.DataFrame({ station_id: df[station_id], date: df[date], longitude: df[longitude], latitude: df[latitude], sunshine_hours: df[sunshine_hours], day_total_radiation_MJ: np.round(rs, 3), day_total_radiation_kWh: np.round(rs / 3.6, 4), })逻辑说明ra的计算整列向量化840 个站点 30 万行数据在一秒内完成。max_sunshine是当天理论日照上限用24/π * ws把弧度折算成小时。sunshine_frac是实际日照时数与上限的比值clip 到 0~1 是因为仪器误差会导致个别记录超过理论值。最终rs (0.25 0.50 * sunshine_frac) * ra晴天sunshine_frac1时 Rs 等于 0.75 倍 Ra符合晴空辐射经验上限。参数说明输出同时给 MJ/m²/day 和 kWh/m²/day 两列光伏设计常用后者省得二次换算。rs / 3.6是 MJ 到 kWh 的换算关系因为 1 kWh 3.6 MJ。3.4 结果导出保留原始字段与编码处理输出 CSV 时把原始站点信息一并保留方便后续按站点分组统计。编码用utf-8-sig是为了让 Excel 打开时中文表头不乱码这是经常被忽略但实际很影响交付体验的细节。def main(input_path, output_path): df load_data(input_path) result calc_radiation(df) result.to_csv(output_path, indexFalse, encodingutf-8-sig) print(f处理完成: {len(result)} 条记录) print(f日总辐射范围: {result[day_total_radiation_MJ].min():.2f} f~ {result[day_total_radiation_MJ].max():.2f} MJ/m2/day)逻辑说明main函数把数据加载、辐射计算、导出串成完整流程print输出的最小值和最大辐射值是最快的 sanity check——如果最大值超过 40 MJ/m²/day大概率是纬度或日序计算出了问题这个阈值在中国区域基本不会误报。4. 运行与调试流程命令行、输出文件与结果验证4.1 环境准备与最小依赖这个脚本只依赖 pandas 和 numpy没有其他第三方库Python 3.9 以上直接跑。虚拟环境安装命令python -m venv .venv source .venv/bin/activate # Windows 下是 .venv\Scripts\activate pip install pandas numpypandas 负责读取和输出numpy 负责向量化数学运算。如果你的源数据是 Excel 格式再加一个pip install openpyxl就行读取时pd.read_excel替代pd.read_csv。整个环境装下来不到 200MB和动不动就装全套地理库的方案比这份脚本在普通办公笔记本上就能跑。4.2 命令行参数设计脚本入口支持四个关键参数列名映射在代码里改日常运行不需要动代码。python sunshine_to_radiation.py \ --input 全国840站点_日照时数_2023.csv \ --date-col date \ --lat-col latitude \ --sun-col sunshine_hours \ --output out/daily_radiation_2023.csv逻辑说明--input指向原始 CSV 路径--date-col、--lat-col、--sun-col指定源文件里的日期、纬度、日照时数列名适配不同来源的数据集--output是结果路径。这样做的好处是换一批数据只需要改命令行参数不需要打开脚本改代码——对于经常换数据源的气象数据处理场景来说这个设计能省下大量重复劳动。参数字段按argparse实现缺省值分别是date、latitude、sunshine_hours通常不用传。4.3 输出文件长什么样结果 CSV 每行是一个站点一天的数据列结构如下字段类型说明station_idstr站点编号datedatetime日期longitude / latitudefloat经纬度sunshine_hoursfloat原始日照时数day_total_radiation_MJfloat日总辐射MJ/m²/dayday_total_radiation_kWhfloat日总辐射kWh/m²/dayday_total_radiation_MJ是算法直接输出第二个 kWh 列是换算值。用 Excel 打开后先看最大值和最小值——中国区域内日总辐射的合理范围大致是 0~35 MJ/m²/day如果出现 50 甚至 100 的数字说明公式里某个环节出了问题。4.4 结果验证三步走拿到输出别急着用先做三件事验证。第一步单站抽检。挑一个纬度居中的站点把这一站全年的日总辐射求和得到年总量单位 MJ/m²应该在 4000~7000 之间。低于 3000 说明日照百分率被严重低估可能原始日照时数单位是分钟高于 8000 则要查天文辐射计算。第二步季节趋势检查。北半球站点 6 月的日均辐射应显著高于 12 月如果出现冬季反而高的倒挂现象大概率是日序公式里加减天数出了错。第三步与公开参考数据对照。选 3~5 个站点把多年平均的月总辐射和某再分析数据集或气候公报里的数值比对误差在 10% 以内算正常。气候态的 Angstrom 系数本身有区域差异超过 15% 的偏差优先怀疑系数选取是否适合本地。5. 避坑指南与常见问题排查单位、角度、缺失值与高原站点5.1 数据读取类问题NaN 与全零输出现象脚本跑完没有报错但输出文件里的辐射值全是 NaN或者某一站全为 0。原因最常见的是日期字段解析失败pd.to_datetime遇到2023/1/1 00:00这类半小时粒度时间戳时会产生 NaTdayofyear随之变成 NaN整行乘出来全是 NaN。另一个常见原因是日照时数列存在---或9999这类缺失值标记pd.to_numeric如果没加errorscoerce整个列会变成字符串类型运算结果全部是 NaN。解决load_data里强制pd.to_numeric(errorscoerce)并在读取后打印缺失统计。对于“全为 0”的情况直接检查原始 CSV 里该站点是否确实有观测部分站点冬季日照时数记录为 0 是正常的不一定有错。可以在输出前加一行print(df.groupby(station_id)[sunshine_hours].agg([count, min, max]))这个打印能快速发现哪些站点缺数据、哪些站点数值范围异常。5.2 公式计算类问题辐射值反常的大现象某几个站点的日总辐射超过 50 MJ/m²/day明显超出物理上限。原因纬度没有转弧度。np.tan(lat_rad)和np.sin(lat_rad)用的必须是弧度如果直接用df[latitude]的度数参与运算在纬度 30°~50° 时误差能达到 30% 以上。另一种可能是太阳赤纬公式里外层少了np.radiansdelta变成了一组很小的数值天文辐射 Ra 虚高。解决把calc_solar_params里所有三角函数输入统一打印出来检查量级。sin(delta)应该在 -0.4~0.4 之间ws应该在 0~3.2 之间如果ws出现复数或大于 π 的值基本就是角度和弧度混用了。我在代码里会把lat_rad、delta、ws三个量各打印一行 min/max跑完扫一眼心理就有底。5.3 高纬度站点的边界问题现象纬度高于 66° 的站点在冬季出现 NaN夏季辐射值却比邻近站点高出一截。原因日落时角公式ws arccos(-tan(phi) * tan(delta))在极昼极夜区域会出现定义域越界——冬天-tan(phi)*tan(delta)可能小于 -1arccos直接返回 NaN夏天它又可能大于 1。没有 clip 的版本会静默产出空值不会报错这最坑人。解决np.clip(-np.tan(lat_rad) * np.tan(delta), -1.0, 1.0)这行代码就是在兜这个底。极夜期间日照时数为 0辐射本身应该是 0clip 后 arccos 返回 πN 和 Ra 算出来是一个合理的小值再乘以 0 的日照百分率得到 0逻辑自洽。5.4 云南贵州等低日照区域系统性偏低现象某几个站点的年总辐射比气候公报低 20% 以上但不是所有站点都这样。原因Angstrom 系数 as0.25、bs0.50 是在特定区域的多年平均拟合结果对四川盆地、贵州这些多云地区并不适用。这些区域日照时数普遍偏低但实际散射辐射占比高总辐射未必像日照时数显示的那么低。经验公式在这种气候区先天吃亏不是代码 bug。解决把 as 和 bs 暴露成命令行参数--as-coef 0.26 --bs-coef 0.42用本地站的多年实测辐射率定系数。如果没有实测数据至少换用高原地带常用的 as0.18、bs0.55 试跑一遍对比差异。这类调整属于参数校准问题不该在代码里硬编码。5.5 840 个站点跑得极慢的问题现象脚本跑了几小时才出结果CPU 占用还不到 30%。原因用 for 循环逐行遍历 DataFrame 计算三角函数。30 万行数据在 Python 的纯循环里要跑 20 分钟以上而同样的计算用 numpy 向量化一秒钟内完成。解决脚本里calc_solar_params和calc_radiation全部用 numpy 数组整列运算没有任何一层 Python 循环。如果你拿到别的版本写的循环实现重构时只需要记住一条原则——所有涉及三角函数、幂运算、乘除的逐行计算一律先.to_numpy()再整列操作最后拼回 DataFrame。从apply改成向量化840 个站点一年的数据从 20 分钟压缩到 3 秒这个提升立竿见影。6. 二次开发思路把批量脚本改造成可配置的辐射计算工具项目源码只是个能跑通的最小实现实际开发里 840 个站点的数据往往是按月分文件存放还可能有几千个文件需要批量处理需要考虑把脚本封装成可扩展的工具形态。我的建议是把核心函数重构成一个类把 as/bs 系数和单位换算都作为实例属性这样二次开发时只需实例化不同配置就能服务多种数据源class RadiationConverter: def __init__(self, as_coef0.25, bs_coef0.50, output_unitMJ): self.as_coef as_coef self.bs_coef bs_coef self.output_unit output_unit def convert_file(self, input_path, output_path, date_coldate, lat_collatitude, sun_colsunshine_hours): df load_data(input_path) result calc_radiation(df, self.as_coef, self.bs_coef) if self.output_unit kWh: result[day_total_radiation] result[day_total_radiation_MJ] / 3.6 result.to_csv(output_path, indexFalse, encodingutf-8-sig)这个封装把硬编码的 0.25 和 0.50 变成了构造参数把年份、区域、数据源差异全部隔离在实例配置里。再往下可以加一个config.yaml描述各数据源的列名映射用yaml.safe_load读取后自动适配不同机构提供的数据省得每次拿到新数据都改代码。另一个值得投入的改造方向是增加质量诊断模式。现在脚本只管算不管数据质量实际跑气象数据时缺测、异常值、站点漂移才是真正的血泪来源。加一个--validate开关输出每个站点的有效记录数、日照时数年总量、平均日辐射跑完数据回头看一眼就能定位问题这比算完再检查要高效得多。我还要加一个功能多文件批量处理。840 个站点如果按站点分文件存放外层用一个for循环遍历目录每个文件计算完追加写出避免一次性加载到内存里卡死。最后再加个--skip-errors参数某个文件格式有问题时跳过并记录日志不让整个批次中断——这在实际生产中比任何算法优化都更实用。从那以后我每次拿到气象数据都会先跑一遍describe()和isna().sum()确认单位再开始算算完立刻看最大最小值范围不再急着把结果交给下游。这套流程救过我好几次希望帮到你。本文还有配套的精品资源点击获取