简介本资源是一套面向水文、气象及农业工程领域科研人员与高校师生的潜在蒸散发ET₀计算Python工具集覆盖25种主流经验与半理论公式包括Penman-Monteith、Hargreaves-Samani、Thornthwaite、Priestley-Taylor等经典方法适用于不同数据可得性条件下的区域蒸散发估算与模型比选。压缩包共11个文件含5个核心Python源码.py与6个已编译字节码.pyc其中et_methods.py封装全部算法逻辑utils.py与converter.py提供单位转换与输入校验支持global_variables.py统一管理参数常量结构清晰、模块解耦便于二次开发与方法替换。资源体积仅117KB轻量易集成已获1014人学习下载。使用者可直接调用各方法函数进行批量计算快速对比不同模型在本地气候数据下的适用性显著降低从公式推导到代码实现的技术门槛。1. 项目概述从气象公式到一行代码搞气象、水文、农业或者生态研究的朋友对“潜在蒸散发”这个概念肯定不陌生。简单说它指的是在水分供应充足条件下地表可能蒸发和植物可能蒸腾的总水量。听起来是个纯理论的气象学参数对吧但它的实际应用场景广得吓人——从判断一个地区是不是干旱缺水到给农田灌溉定个科学的浇水计划再到预测气候变化对水资源的影响都离不开它。以前算这玩意儿要么依赖昂贵专业的气象软件要么就得抱着厚厚的公式手册手动按计算器过程繁琐还容易出错。但现在不一样了我们有了Python。今天要聊的就是怎么用Python把那些复杂的Penman-Monteith、Priestley-Taylor公式变成几行清晰、可复用的代码。这不仅仅是“写个脚本”而是搭建一个属于自己的、灵活强大的水文气象分析工具链的开端。无论你是刚接触科研的学生还是需要快速处理数据的工程师这套代码都能让你从重复劳动中解放出来把精力真正放在数据分析和问题洞察上。2. 核心原理与公式选择为什么是Penman-Monteith在动手写代码之前得先搞清楚我们要算的是什么以及为什么选这个公式。潜在蒸散发ET0的估算方法很多从简单的温度法到复杂的能量平衡法选择哪种直接决定了代码的复杂度和结果的可靠性。2.1 主流计算方法对比为了有个直观认识我们先看一个简单的对比方法核心原理所需数据优点缺点适用场景温度法如Hargreaves基于日均温与温差的经验公式日均温、最高温、最低温数据需求少计算简单精度相对较低区域性较强数据稀缺的初步估算、历史气候分析辐射法如Priestley-Taylor基于净辐射的能量平衡简化版净辐射、气温物理基础较好无需风速湿度在干旱地区可能低估包含经验系数湿润地区、下垫面均匀如茂密草地综合法Penman-Monteith结合能量平衡与空气动力学过程气温、湿度、风速、辐射、气压物理机制最完备精度高国际推荐标准数据需求多计算复杂高精度研究、业务化应用、国际对比2.2 为什么首选FAO Penman-Monteith公式联合国粮农组织FAO推荐的Penman-Monteith公式是目前国际上公认的标准方法。它不是一个黑箱其核心思想是同时考虑了能量供给太阳辐射提供蒸发的动力和空气干燥力风速和湿度差带走水汽的能力这两个关键过程。公式看起来有点吓人但拆解后很好理解ET0 [Δ·(Rn - G) γ·(900/(T273))·u2·(es - ea)] / [Δ γ·(1 0.34·u2)]其中Δ饱和水汽压曲线斜率。它随温度变化温度越高空气能容纳的水汽越多斜率越大蒸发潜力越强。这部分是“能量项”的关键系数。Rn地表净辐射。这是蒸发能量的根本来源等于太阳短波辐射收入减去地面长波辐射支出。G土壤热通量。对于日尺度计算通常可以忽略或用一个简单比例估算。γ干湿表常数。和大气压有关海拔越高气压越低γ值略增。T, u2, es, ea分别是气温℃、2米高风速m/s、饱和水汽压kPa、实际水汽压kPa。(es - ea)就是饱和差直接体现了空气的“干燥力”。注意公式里的常数900和0.34是FAO针对“参考作物”高度0.12m的绿草地校准后的结果。如果你的研究对象是其他植被如森林、作物可能需要调整这些参数这就是所谓的“作物系数”Kc校正实际蒸散发 ETc Kc * ET0。选择实现这个公式意味着我们的代码将具备专业级的精度和广泛的认可度。虽然需要的气象要素多一点但如今很多公开数据集如NASA POWER ERA5都能提供这使得复现标准流程成为可能。3. 代码架构设计与数据准备直接写一个几百行的函数把所有东西塞进去是初学者的做法。好的工程代码应该是模块化的、清晰的、易于调试和扩展的。我们的计算程序可以分成几个核心模块。3.1 模块化设计思路数据输入与校验模块负责读取原始数据如CSV、Excel、NetCDF检查数据完整性有无缺测处理异常值并统一单位如辐射从W/m²换算为MJ/m²/day。气象参数计算模块这是核心。将原始观测数据气温、露点温度、风速、辐射转换为公式所需的中间变量如饱和水汽压es、实际水汽压ea、饱和水汽压曲线斜率Δ、干湿表常数γ等。每个变量一个函数。辐射计算模块如果无法直接获得净辐射Rn则需要根据日照时数、纬度、日期等计算太阳辐射Ra和净短波辐射Rns、净长波辐射Rnl最后得到Rn。这部分涉及天文计算独立出来逻辑更清晰。主计算模块整合以上所有中间变量按照FAO Penman-Monteith公式完成最终ET0的计算。结果输出与可视化模块将计算结果保存为文件并生成时间序列图、月统计图等直观展示分析结果。3.2 关键数据获取与预处理“垃圾进垃圾出。”数据的质量直接决定结果的可靠性。你需要准备至少日尺度的以下数据最高气温Tmax与最低气温Tmin用于计算日均温Tmean (TmaxTmin)/2以及饱和水汽压。相对湿度RH或露点温度Tdew用于计算实际水汽压ea。如果只有平均相对湿度RH_mean则ea (es(Tmin) es(Tmax))/2 * RH_mean / 100。使用露点温度直接计算ea更准确。风速u通常需要2米高度的风速u2。如果风速仪高度不同需按对数风速廓线公式换算到2米高。太阳辐射Rs或日照时数n这是计算净辐射Rn的关键。如果有实测太阳辐射Rs单位MJ/m²/day最好。如果没有FAO提供了通过日照时数和天文辐射Ra估算Rs的公式。站点纬度lat和海拔alt用于计算太阳辐射、大气压等。实操心得数据源方面NASA POWER数据集是免费且易用的宝藏它提供了全球任意位置基于再分析的气象数据包括直接可用的太阳辐射格式规整非常适合研究和初步应用。对于中国区域国家气象科学数据中心提供的地面气象站数据更精确但需要申请和处理。在代码中务必为每个输入参数添加详细的注释说明其单位和来源并编写数据有效性检查如温度范围是否合理湿度是否在0-100之间这能避免很多隐蔽的错误。4. 核心函数实现与代码逐行解析接下来我们进入核心环节用Python实现FAO-56 Penman-Monteith公式。我们将遵循模块化原则逐个击破。4.1 基础气象参数计算函数这些函数是构建公式的“砖块”。我们假设输入数据都是日尺度的NumPy数组或Pandas Series。import numpy as np import pandas as pd def mean_saturation_vapor_pressure(Tmin, Tmax): 计算日平均饱和水汽压 (es, kPa)。 根据FAO采用日最高温和最低温对应的饱和水汽压的平均值。 # 饱和水汽压计算公式 (Tetens公式) es_Tmin 0.6108 * np.exp((17.27 * Tmin) / (Tmin 237.3)) es_Tmax 0.6108 * np.exp((17.27 * Tmax) / (Tmax 237.3)) es (es_Tmin es_Tmax) / 2.0 return es def actual_vapor_pressure_from_rh(Tmin, Tmax, RH_mean): 通过平均相对湿度计算实际水汽压 (ea, kPa)。 RH_mean: 平均相对湿度 (%) es mean_saturation_vapor_pressure(Tmin, Tmax) ea es * (RH_mean / 100.0) return ea def actual_vapor_pressure_from_dewpoint(Tdew): 通过露点温度计算实际水汽压 (ea, kPa)。更推荐的方法。 ea 0.6108 * np.exp((17.27 * Tdew) / (Tdew 237.3)) return ea def slope_of_saturation_vapor_pressure(Tmean): 计算饱和水汽压曲线斜率 (Δ, kPa/°C)。 Tmean: 日均温 (°C) numerator 4098 * (0.6108 * np.exp((17.27 * Tmean) / (Tmean 237.3))) denominator (Tmean 237.3) ** 2 delta numerator / denominator return delta def atmospheric_pressure(altitude): 计算大气压 (P, kPa)。 altitude: 海拔高度 (m) P 101.3 * ((293 - 0.0065 * altitude) / 293) ** 5.26 return P def psychrometric_constant(P): 计算干湿表常数 (γ, kPa/°C)。 P: 大气压 (kPa) gamma 0.665e-3 * P # 常数0.000665是比热容等参数的组合 return gamma4.2 辐射计算模块这是难点也是容易出错的地方。我们实现一个完整的、由日照时数推算净辐射的流程。def extraterrestrial_radiation(lat, doy): 计算日地外辐射 (Ra, MJ/m²/day)。 lat: 纬度 (度北纬为正) doy: 年日序 (1-365/366) lat_rad np.radians(lat) # 转换为弧度 # 太阳磁偏角 (δ) delta 0.409 * np.sin((2 * np.pi / 365) * doy - 1.39) # 日地距离倒数 (dr) dr 1 0.033 * np.cos(2 * np.pi / 365 * doy) # 日落时角 (ωs) omega_s np.arccos(-np.tan(lat_rad) * np.tan(delta)) # 地外辐射 Ra (公式简化版) Ra (24 * 60 / np.pi) * 0.0820 * dr * ( omega_s * np.sin(lat_rad) * np.sin(delta) np.cos(lat_rad) * np.cos(delta) * np.sin(omega_s) ) return Ra def net_radiation(Tmin, Tmax, Rs, lat, doy, albedo0.23): 计算地表净辐射 (Rn, MJ/m²/day)。 Rs: 入射太阳短波辐射 (MJ/m²/day)。若无可通过日照时数估算。 albedo: 地表反照率参考草地下默认0.23。 # 1. 计算净短波辐射 (Rns) Rns (1 - albedo) * Rs # 2. 计算净长波辐射 (Rnl) # 斯蒂芬-玻尔兹曼常数 (MJ/K⁴/m²/day) sigma 4.903e-9 # 通过最高最低温估算实际水汽压 (简化处理) ea actual_vapor_pressure_from_rh(Tmin, Tmax, (TminTmax)/2) # 此处用平均温估算湿度不精确仅示例 # 绝对温度 Tmax_K Tmax 273.16 Tmin_K Tmin 273.16 # 净长波辐射公式 Rnl sigma * ((Tmax_K**4 Tmin_K**4)/2) * (0.34 - 0.14 * np.sqrt(ea)) * (1.35 * (Rs / extraterrestrial_radiation(lat, doy)) - 0.35) # 确保Rnl为负值或零能量损失 Rnl np.where(Rnl 0, 0, Rnl) # 3. 净辐射 Rn Rn Rns Rnl return Rn重要提示net_radiation函数中的长波辐射计算部分对ea很敏感。示例中用一个粗略估算来演示流程。在实际应用中务必使用更准确的实际水汽压数据如来自露点温度或相对湿度。此外(1.35 * (Rs/Ra) - 0.35)是云量影响的修正因子当 Rs/Ra日照百分率数据质量差时这个估算误差会被放大。4.3 主计算函数整合现在我们把所有“砖块”垒起来建成计算ET0的“房子”。def calculate_et0_fao56(Tmin, Tmax, RH_mean, wind_speed, Rs, lat, doy, altitude): FAO-56 Penman-Monteith 日潜在蒸散发计算主函数。 参数均为日尺度序列NumPy数组或标量。 返回ET0 (mm/day)。 # 1. 计算中间变量 Tmean (Tmin Tmax) / 2.0 es mean_saturation_vapor_pressure(Tmin, Tmax) ea actual_vapor_pressure_from_rh(Tmin, Tmax, RH_mean) # 建议替换为更精确的ea计算函数 delta slope_of_saturation_vapor_pressure(Tmean) P atmospheric_pressure(altitude) gamma psychrometric_constant(P) Rn net_radiation(Tmin, Tmax, Rs, lat, doy) G 0 # 日尺度土壤热通量通常忽略 # 2. 确保风速为2米高这里假设输入已是u2 u2 wind_speed # 3. FAO Penman-Monteith 公式 numerator_part1 delta * (Rn - G) numerator_part2 gamma * (900 / (Tmean 273)) * u2 * (es - ea) denominator delta gamma * (1 0.34 * u2) ET0 (numerator_part1 numerator_part2) / denominator return ET05. 完整工作流示例与结果分析有了核心函数我们需要一个完整的数据处理流程来驱动它。这里用一个模拟的CSV数据文件为例。5.1 数据读取与预处理实战假设我们有一个weather_data.csv文件包含以下列Date,Tmax_C,Tmin_C,RH_mean,WindSpeed_ms,Sunshine_hours,Latitude,Altitude_m。import pandas as pd import matplotlib.pyplot as plt # 1. 读取数据 df pd.read_csv(weather_data.csv, parse_dates[Date]) df[DOY] df[Date].dt.dayofyear # 计算年日序 # 2. 数据清洗与检查 print(df.describe()) print(df.isnull().sum()) # 处理缺失值对于短时间缺测可以用前后均值插值长时间缺失需谨慎。 df.fillna(methodffill, inplaceTrue) # 前向填充根据情况选择 # 3. 单位转换与衍生变量计算 # 假设风速已是2米高风速。将日照时数转换为太阳辐射 (Rs)。 # 使用FAO的日照时数-辐射转换公式 lat df[Latitude].iloc[0] # 假设站点纬度固定 alt df[Altitude_m].iloc[0] # 假设站点海拔固定 def sunshine_to_radiation(n, lat, doy): 将日照时数 (n, 小时) 转换为太阳辐射 (Rs, MJ/m²/day) Ra extraterrestrial_radiation(lat, doy) # 理论最大日照时数 (N) 计算略复杂此处简化使用FAO近似公式 # 更精确的N计算需要日落时角为简化假设一个近似值或调用完整函数 # 这里演示直接使用一个简化转换系数不精确仅示意 # 实际应用中应使用完整的Angstrom公式或直接获取Rs数据 Rs (0.25 0.5 * (n / 12)) * Ra # 这是一个非常粗略的估算 return Rs # 注意上述 sunshine_to_radiation 函数是高度简化的仅用于演示流程。 # 强烈建议使用实测Rs或更精确的模型如使用extraterrestrial_radiation和日落时角计算N再用Angstrom公式。 # 本例中我们假设数据中已有Rs_MJ列或者我们用其他方式获得了准确的Rs。 # 为了继续演示我们假设df已经有了正确的Rs_MJ列。 # 4. 调用主计算函数 # 假设我们已经有了Rs_MJ列 df[ET0_mm_day] calculate_et0_fao56( Tmindf[Tmin_C].values, Tmaxdf[Tmax_C].values, RH_meandf[RH_mean].values, wind_speeddf[WindSpeed_ms].values, Rsdf[Rs_MJ].values, # 使用假设存在的太阳辐射列 latlat, doydf[DOY].values, altitudealt ) # 5. 查看结果 print(df[[Date, ET0_mm_day]].head()) print(f年均ET0: {df[ET0_mm_day].mean():.2f} mm/day) print(f月总ET0: \n{df.groupby(df[Date].dt.month)[ET0_mm_day].sum()})5.2 结果可视化与解读计算出的ET0是一个时间序列可视化能帮助我们快速发现规律和异常。# 绘制ET0时间序列 plt.figure(figsize(14, 6)) plt.plot(df[Date], df[ET0_mm_day], b-, linewidth0.8, labelDaily ET0) # 绘制月平均线 df[YearMonth] df[Date].dt.to_period(M) monthly_avg df.groupby(YearMonth)[ET0_mm_day].mean() # 需要将Period索引转换为绘图可用的日期 monthly_avg.index monthly_avg.index.to_timestamp() plt.plot(monthly_avg.index, monthly_avg.values, r-, linewidth2, labelMonthly Avg) plt.xlabel(Date) plt.ylabel(ET0 (mm/day)) plt.title(Potential Evapotranspiration (FAO-56 PM) Time Series) plt.legend() plt.grid(True, alpha0.3) plt.tight_layout() plt.savefig(et0_time_series.png, dpi300) plt.show() # 绘制月累计ET0柱状图 monthly_total df.groupby(df[Date].dt.month)[ET0_mm_day].sum() months [Jan, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, Nov, Dec] plt.figure(figsize(10, 5)) plt.bar(months, monthly_total) plt.xlabel(Month) plt.ylabel(Cumulative ET0 (mm)) plt.title(Monthly Total Potential Evapotranspiration) plt.tight_layout() plt.savefig(et0_monthly_total.png, dpi300) plt.show()通过图表你可以清晰地看到ET0的季节性变化夏季高冬季低以及由天气波动如阴雨天辐射低导致ET0骤降引起的日际变化。将月累计ET0与降水量对比就能初步评估该地区的干湿状况。6. 常见陷阱、调试技巧与进阶优化即使公式和代码都正确在实际运行中你还是会遇到各种问题。下面是一些我踩过的坑和解决方案。6.1 数据质量导致的典型问题ET0出现负值或异常高值检查辐射数据这是最常见的原因。确保你的太阳辐射Rs单位是MJ/m²/day而不是 W/m²。1 W/m² 0.0864 MJ/m²/day。单位错误会导致Rn计算差两个数量级。检查温度范围确认Tmax和Tmin没有单位错误如华氏度当摄氏度用。检查是否有非物理值如Tmin Tmax。检查湿度数据实际水汽压ea不能大于饱和水汽压es。如果使用RH计算确保RH在0-100%之间。如果(es - ea)为负会导致公式第二部分为负可能产生不合理结果。结果序列存在NaN检查分母为零公式分母是Δ γ * (1 0.34 * u2)。Δ和γ通常为正但需检查输入数据中是否有导致Δ计算异常的温度值如极端低温。检查辐射计算中的除零在net_radiation函数中Rs / Ra可能导致除零虽然Ra理论上不为零但需确保doy和lat参数合理。6.2 代码调试与验证策略分步验证不要一次性运行整个流程。单独测试每个辅助函数。输入几个已知值手动计算饱和水汽压es、斜率Δ与代码输出对比。找一个计算器或权威软件如FAO提供的ET0 Calculator Excel版作为基准。准备一套标准输入数据FAO文档的示例运行你的代码逐项对比中间变量Rn, es, ea, Δ, γ和最终ET0结果。确保误差在可接受范围通常相对误差1%。敏感性分析改变某个输入参数如温度增加1°C风速增加0.5 m/s观察ET0的变化量是否与理论预期相符温度对ET0影响最大其次是辐射和湿度风速影响相对较小。这能帮你理解模型行为和发现潜在错误。6.3 性能优化与工程化建议当需要处理多年、多站点的数据时效率就很重要了。向量化操作我们上面的代码已经使用了NumPy数组操作天然是向量化的比用for循环快几个数量级。确保传入函数的都是数组而不是在函数内部循环。处理大型网格数据如果你有NetCDF格式的再分析数据如ERA5可以使用xarray库它能够优雅地处理多维数组和坐标并支持分块计算避免内存溢出。import xarray as xr # 假设ds是一个包含所有气象变量的xarray Dataset # 你可以将我们的函数改写为支持xarray的ufunc或者直接对每个网格点应用并行计算对于成百上千个站点的计算可以使用multiprocessing或joblib库进行并行处理。构建命令行工具或Web应用使用argparse库封装你的脚本使其可以通过命令行参数指定输入文件、输出路径和站点信息。或者用Streamlit或Dash快速构建一个交互式Web应用让不熟悉代码的同事也能上传数据并查看结果。6.4 扩展思考从ET0到实际用水管理计算ET0只是第一步。真正的价值在于应用。作物需水量ETc记住ETc Kc * ET0。你需要找到或校准目标作物在不同生长阶段的作物系数Kc。FAO-56文档附录提供了大量作物的参考Kc值。灌溉调度结合土壤水分平衡模型根据ETc、有效降雨和土壤有效持水量就能制定科学的灌溉计划告诉你“什么时候浇”和“浇多少”。干旱监测计算标准化降水蒸散指数SPEI它同时考虑了降水和潜在蒸散发是比SPI更综合的干旱指标。把这段Python代码作为起点你已经拥有了一个核心的水文气象分析引擎。围绕它构建数据管道、可视化界面和应用模型就能解决许多实际的科研和工程问题。代码本身不难难的是对物理过程的理解和对数据质量的把控。多验证多思考这个工具会成为你研究工作中非常得力的助手。本文还有配套的精品资源点击获取 
长河君 · 2026/9/3 5:10:32 · 阅读 1.2k Python实现FAO Penman-Monteith公式：精准计算潜在蒸散发（ET0）
简介本资源是一套面向水文、气象及农业工程领域科研人员与高校师生的潜在蒸散发ET₀计算Python工具集覆盖25种主流经验与半理论公式包括Penman-Monteith、Hargreaves-Samani、Thornthwaite、Priestley-Taylor等经典方法适用于不同数据可得性条件下的区域蒸散发估算与模型比选。压缩包共11个文件含5个核心Python源码.py与6个已编译字节码.pyc其中et_methods.py封装全部算法逻辑utils.py与converter.py提供单位转换与输入校验支持global_variables.py统一管理参数常量结构清晰、模块解耦便于二次开发与方法替换。资源体积仅117KB轻量易集成已获1014人学习下载。使用者可直接调用各方法函数进行批量计算快速对比不同模型在本地气候数据下的适用性显著降低从公式推导到代码实现的技术门槛。1. 项目概述从气象公式到一行代码搞气象、水文、农业或者生态研究的朋友对“潜在蒸散发”这个概念肯定不陌生。简单说它指的是在水分供应充足条件下地表可能蒸发和植物可能蒸腾的总水量。听起来是个纯理论的气象学参数对吧但它的实际应用场景广得吓人——从判断一个地区是不是干旱缺水到给农田灌溉定个科学的浇水计划再到预测气候变化对水资源的影响都离不开它。以前算这玩意儿要么依赖昂贵专业的气象软件要么就得抱着厚厚的公式手册手动按计算器过程繁琐还容易出错。但现在不一样了我们有了Python。今天要聊的就是怎么用Python把那些复杂的Penman-Monteith、Priestley-Taylor公式变成几行清晰、可复用的代码。这不仅仅是“写个脚本”而是搭建一个属于自己的、灵活强大的水文气象分析工具链的开端。无论你是刚接触科研的学生还是需要快速处理数据的工程师这套代码都能让你从重复劳动中解放出来把精力真正放在数据分析和问题洞察上。2. 核心原理与公式选择为什么是Penman-Monteith在动手写代码之前得先搞清楚我们要算的是什么以及为什么选这个公式。潜在蒸散发ET0的估算方法很多从简单的温度法到复杂的能量平衡法选择哪种直接决定了代码的复杂度和结果的可靠性。2.1 主流计算方法对比为了有个直观认识我们先看一个简单的对比方法核心原理所需数据优点缺点适用场景温度法如Hargreaves基于日均温与温差的经验公式日均温、最高温、最低温数据需求少计算简单精度相对较低区域性较强数据稀缺的初步估算、历史气候分析辐射法如Priestley-Taylor基于净辐射的能量平衡简化版净辐射、气温物理基础较好无需风速湿度在干旱地区可能低估包含经验系数湿润地区、下垫面均匀如茂密草地综合法Penman-Monteith结合能量平衡与空气动力学过程气温、湿度、风速、辐射、气压物理机制最完备精度高国际推荐标准数据需求多计算复杂高精度研究、业务化应用、国际对比2.2 为什么首选FAO Penman-Monteith公式联合国粮农组织FAO推荐的Penman-Monteith公式是目前国际上公认的标准方法。它不是一个黑箱其核心思想是同时考虑了能量供给太阳辐射提供蒸发的动力和空气干燥力风速和湿度差带走水汽的能力这两个关键过程。公式看起来有点吓人但拆解后很好理解ET0 [Δ·(Rn - G) γ·(900/(T273))·u2·(es - ea)] / [Δ γ·(1 0.34·u2)]其中Δ饱和水汽压曲线斜率。它随温度变化温度越高空气能容纳的水汽越多斜率越大蒸发潜力越强。这部分是“能量项”的关键系数。Rn地表净辐射。这是蒸发能量的根本来源等于太阳短波辐射收入减去地面长波辐射支出。G土壤热通量。对于日尺度计算通常可以忽略或用一个简单比例估算。γ干湿表常数。和大气压有关海拔越高气压越低γ值略增。T, u2, es, ea分别是气温℃、2米高风速m/s、饱和水汽压kPa、实际水汽压kPa。(es - ea)就是饱和差直接体现了空气的“干燥力”。注意公式里的常数900和0.34是FAO针对“参考作物”高度0.12m的绿草地校准后的结果。如果你的研究对象是其他植被如森林、作物可能需要调整这些参数这就是所谓的“作物系数”Kc校正实际蒸散发 ETc Kc * ET0。选择实现这个公式意味着我们的代码将具备专业级的精度和广泛的认可度。虽然需要的气象要素多一点但如今很多公开数据集如NASA POWER ERA5都能提供这使得复现标准流程成为可能。3. 代码架构设计与数据准备直接写一个几百行的函数把所有东西塞进去是初学者的做法。好的工程代码应该是模块化的、清晰的、易于调试和扩展的。我们的计算程序可以分成几个核心模块。3.1 模块化设计思路数据输入与校验模块负责读取原始数据如CSV、Excel、NetCDF检查数据完整性有无缺测处理异常值并统一单位如辐射从W/m²换算为MJ/m²/day。气象参数计算模块这是核心。将原始观测数据气温、露点温度、风速、辐射转换为公式所需的中间变量如饱和水汽压es、实际水汽压ea、饱和水汽压曲线斜率Δ、干湿表常数γ等。每个变量一个函数。辐射计算模块如果无法直接获得净辐射Rn则需要根据日照时数、纬度、日期等计算太阳辐射Ra和净短波辐射Rns、净长波辐射Rnl最后得到Rn。这部分涉及天文计算独立出来逻辑更清晰。主计算模块整合以上所有中间变量按照FAO Penman-Monteith公式完成最终ET0的计算。结果输出与可视化模块将计算结果保存为文件并生成时间序列图、月统计图等直观展示分析结果。3.2 关键数据获取与预处理“垃圾进垃圾出。”数据的质量直接决定结果的可靠性。你需要准备至少日尺度的以下数据最高气温Tmax与最低气温Tmin用于计算日均温Tmean (TmaxTmin)/2以及饱和水汽压。相对湿度RH或露点温度Tdew用于计算实际水汽压ea。如果只有平均相对湿度RH_mean则ea (es(Tmin) es(Tmax))/2 * RH_mean / 100。使用露点温度直接计算ea更准确。风速u通常需要2米高度的风速u2。如果风速仪高度不同需按对数风速廓线公式换算到2米高。太阳辐射Rs或日照时数n这是计算净辐射Rn的关键。如果有实测太阳辐射Rs单位MJ/m²/day最好。如果没有FAO提供了通过日照时数和天文辐射Ra估算Rs的公式。站点纬度lat和海拔alt用于计算太阳辐射、大气压等。实操心得数据源方面NASA POWER数据集是免费且易用的宝藏它提供了全球任意位置基于再分析的气象数据包括直接可用的太阳辐射格式规整非常适合研究和初步应用。对于中国区域国家气象科学数据中心提供的地面气象站数据更精确但需要申请和处理。在代码中务必为每个输入参数添加详细的注释说明其单位和来源并编写数据有效性检查如温度范围是否合理湿度是否在0-100之间这能避免很多隐蔽的错误。4. 核心函数实现与代码逐行解析接下来我们进入核心环节用Python实现FAO-56 Penman-Monteith公式。我们将遵循模块化原则逐个击破。4.1 基础气象参数计算函数这些函数是构建公式的“砖块”。我们假设输入数据都是日尺度的NumPy数组或Pandas Series。import numpy as np import pandas as pd def mean_saturation_vapor_pressure(Tmin, Tmax): 计算日平均饱和水汽压 (es, kPa)。 根据FAO采用日最高温和最低温对应的饱和水汽压的平均值。 # 饱和水汽压计算公式 (Tetens公式) es_Tmin 0.6108 * np.exp((17.27 * Tmin) / (Tmin 237.3)) es_Tmax 0.6108 * np.exp((17.27 * Tmax) / (Tmax 237.3)) es (es_Tmin es_Tmax) / 2.0 return es def actual_vapor_pressure_from_rh(Tmin, Tmax, RH_mean): 通过平均相对湿度计算实际水汽压 (ea, kPa)。 RH_mean: 平均相对湿度 (%) es mean_saturation_vapor_pressure(Tmin, Tmax) ea es * (RH_mean / 100.0) return ea def actual_vapor_pressure_from_dewpoint(Tdew): 通过露点温度计算实际水汽压 (ea, kPa)。更推荐的方法。 ea 0.6108 * np.exp((17.27 * Tdew) / (Tdew 237.3)) return ea def slope_of_saturation_vapor_pressure(Tmean): 计算饱和水汽压曲线斜率 (Δ, kPa/°C)。 Tmean: 日均温 (°C) numerator 4098 * (0.6108 * np.exp((17.27 * Tmean) / (Tmean 237.3))) denominator (Tmean 237.3) ** 2 delta numerator / denominator return delta def atmospheric_pressure(altitude): 计算大气压 (P, kPa)。 altitude: 海拔高度 (m) P 101.3 * ((293 - 0.0065 * altitude) / 293) ** 5.26 return P def psychrometric_constant(P): 计算干湿表常数 (γ, kPa/°C)。 P: 大气压 (kPa) gamma 0.665e-3 * P # 常数0.000665是比热容等参数的组合 return gamma4.2 辐射计算模块这是难点也是容易出错的地方。我们实现一个完整的、由日照时数推算净辐射的流程。def extraterrestrial_radiation(lat, doy): 计算日地外辐射 (Ra, MJ/m²/day)。 lat: 纬度 (度北纬为正) doy: 年日序 (1-365/366) lat_rad np.radians(lat) # 转换为弧度 # 太阳磁偏角 (δ) delta 0.409 * np.sin((2 * np.pi / 365) * doy - 1.39) # 日地距离倒数 (dr) dr 1 0.033 * np.cos(2 * np.pi / 365 * doy) # 日落时角 (ωs) omega_s np.arccos(-np.tan(lat_rad) * np.tan(delta)) # 地外辐射 Ra (公式简化版) Ra (24 * 60 / np.pi) * 0.0820 * dr * ( omega_s * np.sin(lat_rad) * np.sin(delta) np.cos(lat_rad) * np.cos(delta) * np.sin(omega_s) ) return Ra def net_radiation(Tmin, Tmax, Rs, lat, doy, albedo0.23): 计算地表净辐射 (Rn, MJ/m²/day)。 Rs: 入射太阳短波辐射 (MJ/m²/day)。若无可通过日照时数估算。 albedo: 地表反照率参考草地下默认0.23。 # 1. 计算净短波辐射 (Rns) Rns (1 - albedo) * Rs # 2. 计算净长波辐射 (Rnl) # 斯蒂芬-玻尔兹曼常数 (MJ/K⁴/m²/day) sigma 4.903e-9 # 通过最高最低温估算实际水汽压 (简化处理) ea actual_vapor_pressure_from_rh(Tmin, Tmax, (TminTmax)/2) # 此处用平均温估算湿度不精确仅示例 # 绝对温度 Tmax_K Tmax 273.16 Tmin_K Tmin 273.16 # 净长波辐射公式 Rnl sigma * ((Tmax_K**4 Tmin_K**4)/2) * (0.34 - 0.14 * np.sqrt(ea)) * (1.35 * (Rs / extraterrestrial_radiation(lat, doy)) - 0.35) # 确保Rnl为负值或零能量损失 Rnl np.where(Rnl 0, 0, Rnl) # 3. 净辐射 Rn Rn Rns Rnl return Rn重要提示net_radiation函数中的长波辐射计算部分对ea很敏感。示例中用一个粗略估算来演示流程。在实际应用中务必使用更准确的实际水汽压数据如来自露点温度或相对湿度。此外(1.35 * (Rs/Ra) - 0.35)是云量影响的修正因子当 Rs/Ra日照百分率数据质量差时这个估算误差会被放大。4.3 主计算函数整合现在我们把所有“砖块”垒起来建成计算ET0的“房子”。def calculate_et0_fao56(Tmin, Tmax, RH_mean, wind_speed, Rs, lat, doy, altitude): FAO-56 Penman-Monteith 日潜在蒸散发计算主函数。 参数均为日尺度序列NumPy数组或标量。 返回ET0 (mm/day)。 # 1. 计算中间变量 Tmean (Tmin Tmax) / 2.0 es mean_saturation_vapor_pressure(Tmin, Tmax) ea actual_vapor_pressure_from_rh(Tmin, Tmax, RH_mean) # 建议替换为更精确的ea计算函数 delta slope_of_saturation_vapor_pressure(Tmean) P atmospheric_pressure(altitude) gamma psychrometric_constant(P) Rn net_radiation(Tmin, Tmax, Rs, lat, doy) G 0 # 日尺度土壤热通量通常忽略 # 2. 确保风速为2米高这里假设输入已是u2 u2 wind_speed # 3. FAO Penman-Monteith 公式 numerator_part1 delta * (Rn - G) numerator_part2 gamma * (900 / (Tmean 273)) * u2 * (es - ea) denominator delta gamma * (1 0.34 * u2) ET0 (numerator_part1 numerator_part2) / denominator return ET05. 完整工作流示例与结果分析有了核心函数我们需要一个完整的数据处理流程来驱动它。这里用一个模拟的CSV数据文件为例。5.1 数据读取与预处理实战假设我们有一个weather_data.csv文件包含以下列Date,Tmax_C,Tmin_C,RH_mean,WindSpeed_ms,Sunshine_hours,Latitude,Altitude_m。import pandas as pd import matplotlib.pyplot as plt # 1. 读取数据 df pd.read_csv(weather_data.csv, parse_dates[Date]) df[DOY] df[Date].dt.dayofyear # 计算年日序 # 2. 数据清洗与检查 print(df.describe()) print(df.isnull().sum()) # 处理缺失值对于短时间缺测可以用前后均值插值长时间缺失需谨慎。 df.fillna(methodffill, inplaceTrue) # 前向填充根据情况选择 # 3. 单位转换与衍生变量计算 # 假设风速已是2米高风速。将日照时数转换为太阳辐射 (Rs)。 # 使用FAO的日照时数-辐射转换公式 lat df[Latitude].iloc[0] # 假设站点纬度固定 alt df[Altitude_m].iloc[0] # 假设站点海拔固定 def sunshine_to_radiation(n, lat, doy): 将日照时数 (n, 小时) 转换为太阳辐射 (Rs, MJ/m²/day) Ra extraterrestrial_radiation(lat, doy) # 理论最大日照时数 (N) 计算略复杂此处简化使用FAO近似公式 # 更精确的N计算需要日落时角为简化假设一个近似值或调用完整函数 # 这里演示直接使用一个简化转换系数不精确仅示意 # 实际应用中应使用完整的Angstrom公式或直接获取Rs数据 Rs (0.25 0.5 * (n / 12)) * Ra # 这是一个非常粗略的估算 return Rs # 注意上述 sunshine_to_radiation 函数是高度简化的仅用于演示流程。 # 强烈建议使用实测Rs或更精确的模型如使用extraterrestrial_radiation和日落时角计算N再用Angstrom公式。 # 本例中我们假设数据中已有Rs_MJ列或者我们用其他方式获得了准确的Rs。 # 为了继续演示我们假设df已经有了正确的Rs_MJ列。 # 4. 调用主计算函数 # 假设我们已经有了Rs_MJ列 df[ET0_mm_day] calculate_et0_fao56( Tmindf[Tmin_C].values, Tmaxdf[Tmax_C].values, RH_meandf[RH_mean].values, wind_speeddf[WindSpeed_ms].values, Rsdf[Rs_MJ].values, # 使用假设存在的太阳辐射列 latlat, doydf[DOY].values, altitudealt ) # 5. 查看结果 print(df[[Date, ET0_mm_day]].head()) print(f年均ET0: {df[ET0_mm_day].mean():.2f} mm/day) print(f月总ET0: \n{df.groupby(df[Date].dt.month)[ET0_mm_day].sum()})5.2 结果可视化与解读计算出的ET0是一个时间序列可视化能帮助我们快速发现规律和异常。# 绘制ET0时间序列 plt.figure(figsize(14, 6)) plt.plot(df[Date], df[ET0_mm_day], b-, linewidth0.8, labelDaily ET0) # 绘制月平均线 df[YearMonth] df[Date].dt.to_period(M) monthly_avg df.groupby(YearMonth)[ET0_mm_day].mean() # 需要将Period索引转换为绘图可用的日期 monthly_avg.index monthly_avg.index.to_timestamp() plt.plot(monthly_avg.index, monthly_avg.values, r-, linewidth2, labelMonthly Avg) plt.xlabel(Date) plt.ylabel(ET0 (mm/day)) plt.title(Potential Evapotranspiration (FAO-56 PM) Time Series) plt.legend() plt.grid(True, alpha0.3) plt.tight_layout() plt.savefig(et0_time_series.png, dpi300) plt.show() # 绘制月累计ET0柱状图 monthly_total df.groupby(df[Date].dt.month)[ET0_mm_day].sum() months [Jan, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, Nov, Dec] plt.figure(figsize(10, 5)) plt.bar(months, monthly_total) plt.xlabel(Month) plt.ylabel(Cumulative ET0 (mm)) plt.title(Monthly Total Potential Evapotranspiration) plt.tight_layout() plt.savefig(et0_monthly_total.png, dpi300) plt.show()通过图表你可以清晰地看到ET0的季节性变化夏季高冬季低以及由天气波动如阴雨天辐射低导致ET0骤降引起的日际变化。将月累计ET0与降水量对比就能初步评估该地区的干湿状况。6. 常见陷阱、调试技巧与进阶优化即使公式和代码都正确在实际运行中你还是会遇到各种问题。下面是一些我踩过的坑和解决方案。6.1 数据质量导致的典型问题ET0出现负值或异常高值检查辐射数据这是最常见的原因。确保你的太阳辐射Rs单位是MJ/m²/day而不是 W/m²。1 W/m² 0.0864 MJ/m²/day。单位错误会导致Rn计算差两个数量级。检查温度范围确认Tmax和Tmin没有单位错误如华氏度当摄氏度用。检查是否有非物理值如Tmin Tmax。检查湿度数据实际水汽压ea不能大于饱和水汽压es。如果使用RH计算确保RH在0-100%之间。如果(es - ea)为负会导致公式第二部分为负可能产生不合理结果。结果序列存在NaN检查分母为零公式分母是Δ γ * (1 0.34 * u2)。Δ和γ通常为正但需检查输入数据中是否有导致Δ计算异常的温度值如极端低温。检查辐射计算中的除零在net_radiation函数中Rs / Ra可能导致除零虽然Ra理论上不为零但需确保doy和lat参数合理。6.2 代码调试与验证策略分步验证不要一次性运行整个流程。单独测试每个辅助函数。输入几个已知值手动计算饱和水汽压es、斜率Δ与代码输出对比。找一个计算器或权威软件如FAO提供的ET0 Calculator Excel版作为基准。准备一套标准输入数据FAO文档的示例运行你的代码逐项对比中间变量Rn, es, ea, Δ, γ和最终ET0结果。确保误差在可接受范围通常相对误差1%。敏感性分析改变某个输入参数如温度增加1°C风速增加0.5 m/s观察ET0的变化量是否与理论预期相符温度对ET0影响最大其次是辐射和湿度风速影响相对较小。这能帮你理解模型行为和发现潜在错误。6.3 性能优化与工程化建议当需要处理多年、多站点的数据时效率就很重要了。向量化操作我们上面的代码已经使用了NumPy数组操作天然是向量化的比用for循环快几个数量级。确保传入函数的都是数组而不是在函数内部循环。处理大型网格数据如果你有NetCDF格式的再分析数据如ERA5可以使用xarray库它能够优雅地处理多维数组和坐标并支持分块计算避免内存溢出。import xarray as xr # 假设ds是一个包含所有气象变量的xarray Dataset # 你可以将我们的函数改写为支持xarray的ufunc或者直接对每个网格点应用并行计算对于成百上千个站点的计算可以使用multiprocessing或joblib库进行并行处理。构建命令行工具或Web应用使用argparse库封装你的脚本使其可以通过命令行参数指定输入文件、输出路径和站点信息。或者用Streamlit或Dash快速构建一个交互式Web应用让不熟悉代码的同事也能上传数据并查看结果。6.4 扩展思考从ET0到实际用水管理计算ET0只是第一步。真正的价值在于应用。作物需水量ETc记住ETc Kc * ET0。你需要找到或校准目标作物在不同生长阶段的作物系数Kc。FAO-56文档附录提供了大量作物的参考Kc值。灌溉调度结合土壤水分平衡模型根据ETc、有效降雨和土壤有效持水量就能制定科学的灌溉计划告诉你“什么时候浇”和“浇多少”。干旱监测计算标准化降水蒸散指数SPEI它同时考虑了降水和潜在蒸散发是比SPI更综合的干旱指标。把这段Python代码作为起点你已经拥有了一个核心的水文气象分析引擎。围绕它构建数据管道、可视化界面和应用模型就能解决许多实际的科研和工程问题。代码本身不难难的是对物理过程的理解和对数据质量的把控。多验证多思考这个工具会成为你研究工作中非常得力的助手。本文还有配套的精品资源点击获取