简介本资源是一套面向水文、气象及农业科研工作者的潜在蒸散发ET计算Python工具集覆盖25种主流经验与物理模型方法如Penman-Monteith、Hargreaves-Samani、Thornthwaite、Priestley-Taylor等可支撑流域水文模拟、灌溉需水量估算、气候响应分析等实际研究任务。压缩包共11个文件含5个可读可修改的.py源码模块如et_methods.py、utils.py、converter.py与6个已编译的.pyc文件总大小仅117KB轻量易集成适合Python中高级用户快速调用或二次开发。已有1014人学习下载资源结构清晰核心算法封装于et_methods模块单位换算与数据预处理由converter和utils支持global_variables统一管理参数init.py实现包级导入便于按需调用单个模型或批量对比不同方法结果。1. 从气象数据到代码实现为什么我们需要计算潜在蒸散发如果你在农业、水文、生态或者气候变化研究领域工作过那么“潜在蒸散发”这个词一定不陌生。它听起来有点学术但说白了就是在一个理想条件下一片完全被植被覆盖、水分供应充足的地面单位时间内能够蒸发和蒸腾到大气中的总水量。这个“理想条件”很关键它排除了土壤水分不足、植被类型差异等限制因素反映的是纯粹由气象条件比如太阳辐射、温度、湿度、风速决定的蒸发能力。为什么这个“理想值”如此重要因为它是一个基准。在实际应用中比如农田灌溉管理我们知道了潜在蒸散发量再结合土壤墒情和作物系数就能估算出作物实际需要多少水从而制定精准的灌溉计划避免水资源的浪费。在水文模型中它是计算流域水量平衡、预测径流的关键输入参数。在气候变化研究中潜在蒸散发的长期趋势分析是评估干旱风险、生态系统响应的重要指标。然而计算潜在蒸散发并不是一个简单的加减乘除。它背后是复杂的物理过程涉及能量平衡和空气动力学原理。历史上科学家们提出了多种经验或半经验公式来估算它比如彭曼公式、彭曼-蒙蒂斯公式、哈格里夫斯公式、索恩思韦特公式等。这些公式各有优劣有的需要的数据多但精度高有的数据要求低但适用区域有限。过去这些计算往往依赖于专业的商业软件或需要手动查表、套公式过程繁琐且容易出错。而现在Python以其强大的科学计算库和简洁的语法成为了处理这类问题的利器。我们可以用几行代码就自动化地完成从原始气象数据读取、质量检查、公式计算到结果可视化的全过程。这不仅大大提高了工作效率也让研究方法更加透明和可复现。今天我就结合自己处理农业气象数据的经验带你一步步用Python实现几种主流的潜在蒸散发计算方法并分享一些实操中容易踩的坑和优化技巧。2. 核心公式选型面对一堆气象数据我该用哪个公式当你拿到一组气象数据准备计算潜在蒸散发时第一个问题往往是该用哪个公式这不是拍脑袋决定的而是由你手头数据的完整性和精度要求共同决定的。选择不当要么巧妇难为无米之炊要么就是杀鸡用牛刀浪费了高质量数据。下面我梳理了几个最常用的公式及其数据需求你可以对号入座。2.1 彭曼-蒙蒂斯公式当之无愧的“金标准”如果你拥有相对完整的气象站数据那么彭曼-蒙蒂斯公式通常是首选。它被联合国粮农组织推荐是目前理论上最完备、应用最广泛的公式。它综合考虑了净辐射能量项和空气干燥度、风速空气动力项。所需核心数据日均气温最高温、最低温。日照时数或太阳辐射用于计算净辐射。相对湿度或露点温度反映空气湿度。风速通常指2米高处的风速。站点经纬度和海拔用于计算太阳常数、大气压等。为什么选它因为它物理基础扎实在全球多数地区表现稳定。但它的“娇贵”之处在于对数据质量要求高尤其是辐射数据。如果辐射数据缺失或不准计算结果误差会很大。2.2 哈格里夫斯公式数据匮乏时的“救星”在很多情况下尤其是历史数据或偏远地区我们可能只有温度数据。这时哈格里夫斯公式就派上大用场了。它只需要日均最高温、最低温以及站点纬度通过一个经验系数来估算太阳辐射进而计算潜在蒸散发。所需核心数据日均气温最高温、最低温。站点纬度用于估算地外辐射。它的优势与局限优势极其明显——数据需求极简。我在处理一些上世纪的气象数据时它几乎是唯一可行的选择。但它的精度通常低于彭曼-蒙蒂斯公式在非常潮湿或非常干燥的地区可能需要本地化校正系数。2.3 普里斯特利-泰勒公式湿润地区的简化方案这个公式基于能量平衡假设空气动力项可以用一个常数比例α通常取1.26与净辐射项关联。它适用于水分供应充足、大面积均匀的湿润表面比如茂密的森林或灌溉充分的农田。所需核心数据净辐射这是核心输入。气温用于计算饱和水汽压曲线斜率等热力学参数。适用场景当你关注的是能量限制为主的蒸发过程且风速、湿度数据不可靠时可以考虑它。但在干旱半干旱地区它会显著高估蒸散发。为了更直观地对比我整理了一个选型决策表公式名称核心数据需求计算复杂度适用场景主要局限彭曼-蒙蒂斯温度、湿度、风速、辐射、站点信息高数据齐全追求高精度全球多数地区数据要求高计算步骤多哈格里夫斯最高/最低温、纬度低只有温度数据历史数据或偏远地区精度相对较低需地区校准普里斯特利-泰勒净辐射、温度中湿润下垫面能量限制为主干旱区误差大依赖净辐射精度我的经验之谈在实际项目中我通常会做一个“数据审计”。先列出所有可用数据字段及其缺失率。如果辐射、风速数据质量尚可毫不犹豫用彭曼-蒙蒂斯。如果只有温度就用哈格里夫斯先跑出一个基准结果并在报告中明确说明其不确定性。有时我甚至会并行计算两种方法通过对比结果来交叉验证数据的合理性。3. 环境搭建与数据准备别在第一步就掉进坑里工欲善其事必先利其器。一个稳定、可复现的Python环境是后续所有工作的基础。很多人觉得安装包很简单但恰恰是这里隐藏着最多的版本冲突和依赖问题。3.1 创建独立的Python环境我强烈建议使用conda或venv为这个项目创建一个独立的虚拟环境。这能确保你的库版本不会干扰其他项目。# 使用 conda (假设你安装了Anaconda或Miniconda) conda create -n pet_calc python3.9 conda activate pet_calc # 或者使用 venv python -m venv pet_env # Windows pet_env\Scripts\activate # Linux/Mac source pet_env/bin/activate3.2 安装核心计算库在我们的环境中需要安装几个核心库pip install numpy pandas matplotlibnumpy: 数值计算的基石所有公式中的数组运算都靠它。pandas: 数据处理的瑞士军刀读取CSV/Excel、处理时间序列、数据清洗离不开它。matplotlib: 基础绘图库用于可视化结果和输入数据。一个关键的坑scipy的隐式依赖。在计算饱和水汽压、斜率等时我们可能会用到一些数学函数。虽然彭曼-蒙蒂斯公式可以手动实现但为了稳健和方便我们通常会用到scipy中的常量或优化函数。但请注意scipy在某些系统上安装可能因为编译依赖而失败。一个更轻量级的替代是使用metpy或pyet这类气象专用库它们封装了这些计算。这里我们先以纯手工计算为例确保通用性。如果需要可以后续安装pip install scipy3.3 气象数据的读取与清洗假设你有一份名为weather_data.csv的日尺度气象数据其格式可能如下DateTmax_CTmin_CRH_meanWindSpeed_2mSunshine_hoursLatitudeLongitudeElevation_m2023-07-0132.520.165.22.19.540.0116.5502023-07-0234.021.560.82.510.240.0116.550第一步用pandas加载数据import pandas as pd # 读取数据并指定日期列为索引 df pd.read_csv(weather_data.csv, parse_dates[Date], index_colDate) print(df.head()) print(df.info()) # 查看数据概要和缺失值第二步处理缺失值与异常值这是最耗时但也最重要的一步。气象数据常有缺失或明显错误如温度超过60°C。# 1. 简单查看缺失情况 print(df.isnull().sum()) # 2. 对于少量缺失可以用前后值插补需谨慎 # 例如用前一天的湿度填充当天的缺失值仅适用于连续缺失少的情况 df[RH_mean].fillna(methodffill, inplaceTrue) # 3. 对于明显异常值可以设定合理范围进行过滤或标记 # 例如假设研究地点在中国东部日最高温超过45°C或低于-20°C视为异常 df.loc[(df[Tmax_C] 45) | (df[Tmax_C] -20), Tmax_C] pd.NA df.loc[(df[Tmin_C] 35) | (df[Tmin_C] -30), Tmin_C] pd.NA # 4. 删除仍存在关键数据缺失的行例如温度数据缺失就无法计算 # 这里假设Tmax和Tmin是必须的 df.dropna(subset[Tmax_C, Tmin_C], inplaceTrue)注意插补方法需要根据数据特点选择。对于气象序列时间序列插值如线性插值、样条插值可能比简单的前向填充更合理。可以使用df.interpolate(methodtime)。但务必记录下你的处理步骤这在科研中至关重要。第三步单位检查与转换公式计算对单位非常敏感。确保你的数据单位与公式要求一致。常见转换包括温度公式通常使用摄氏度°C。如果你的数据是开尔文K需减去273.15。风速彭曼-蒙蒂斯公式通常需要2米高处的风速m/s。如果你的风速是10米高处或单位是km/h需要转换。辐射最易出错的地方。公式需要的是每日净辐射或太阳短波辐射单位是MJ/m²/day。如果你的日照时数单位是小时需要转换为辐射能量。4. 手把手实现三种主流公式的Python代码详解数据准备好了环境也搭好了现在让我们进入核心环节——编码实现。我将分别实现哈格里夫斯、普里斯特利-泰勒和彭曼-蒙蒂斯公式并解释每一步的物理意义和计算细节。4.1 哈格里夫斯公式实现简约而不简单哈格里夫斯公式的原始形式如下ET0 0.0023 * Ra * (Tmean 17.8) * sqrt(Tmax - Tmin)其中Ra是地外辐射MJ/m²/dayTmean是日均温Tmax和Tmin是日最高/最低温。这里的关键是计算Ra它取决于一年中的第几天DOY和站点纬度。import numpy as np import pandas as pd from math import pi, sin, cos, asin def calculate_ra(doy, lat_deg): 计算地外辐射 Ra (MJ/m²/day) Args: doy (int): 年积日1-365/366 lat_deg (float): 站点纬度度北纬为正南纬为负 Returns: float: 地外辐射 Ra # 1. 将纬度转换为弧度 lat_rad lat_deg * pi / 180.0 # 2. 计算日地相对距离倒数 dr 和太阳磁偏角 delta # 日角 theta 2 * pi * doy / 365.0 # 日地距离倒数 dr 1 0.033 * np.cos(theta) # 太阳磁偏角弧度 delta 0.409 * np.sin(theta - 1.39) # 3. 计算日落时角 omega_s (弧度) # 当 tan(lat)*tan(delta) -1 或 1 时表示极昼或极夜需要特殊处理 tan_term -np.tan(lat_rad) * np.tan(delta) # 限制值在[-1,1]之间避免数学错误 tan_term np.clip(tan_term, -1.0, 1.0) omega_s np.arccos(tan_term) # 4. 计算 Ra # 太阳常数 Gsc 0.0820 MJ/m²/min Gsc 0.0820 # 一天中的分钟数 day_minutes 24 * 60 Ra (day_minutes / pi) * Gsc * dr * ( omega_s * np.sin(lat_rad) * np.sin(delta) np.cos(lat_rad) * np.cos(delta) * np.sin(omega_s) ) return Ra def et0_hargreaves(tmax, tmin, lat, doy): 计算哈格里夫斯潜在蒸散发 ET0 (mm/day) Args: tmax (float): 日最高温 (°C) tmin (float): 日最低温 (°C) lat (float): 纬度 (度) doy (int): 年积日 Returns: float: ET0 (mm/day) tmean (tmax tmin) / 2.0 Ra calculate_ra(doy, lat) # 原始哈格里夫斯公式 et0 0.0023 * Ra * (tmean 17.8) * np.sqrt(tmax - tmin) return et0 # 应用到整个DataFrame # 首先我们需要为每一行计算年积日Day of Year df[DOY] df.index.dayofyear # 假设纬度存储在列Latitude中且所有行纬度相同 lat df[Latitude].iloc[0] # 使用apply函数逐行计算注意sqrt里温差可能为负需处理 df[ET0_Hargreaves] df.apply( lambda row: et0_hargreaves(row[Tmax_C], row[Tmin_C], lat, row[DOY]), axis1 )实操心得np.sqrt(tmax - tmin)这里有个隐患。在极少数情况下如某些海洋性气候tmax可能略低于tmin可能是数据错误或四舍五入导致导致对负数开方报错。一个稳健的做法是加上一个很小的数或取绝对值np.sqrt(np.abs(tmax - tmin) 1e-10)。另外FAO后来推荐了修正的哈格里夫斯公式系数不同如果需要更高精度可以查阅FAO-56手册。4.2 普里斯特利-泰勒公式实现抓住能量核心普里斯特利-泰勒公式ET0 alpha * (delta / (delta gamma)) * (Rn - G) / lambda其中alpha经验系数通常取1.26湿润地区。delta饱和水汽压曲线斜率kPa/°C。gamma干湿表常数kPa/°C。Rn地表净辐射MJ/m²/day。G土壤热通量MJ/m²/day对于日尺度通常近似为0。lambda水的汽化潜热约2.45 MJ/kg。def calculate_delta(tmean): 计算饱和水汽压曲线斜率 delta (kPa/°C) FAO-56 推荐公式 # 计算在温度Tmean下的饱和水汽压 e_s_t 0.6108 * np.exp((17.27 * tmean) / (tmean 237.3)) delta (4098 * e_s_t) / ((tmean 237.3) ** 2) return delta def calculate_gamma(pressure): 计算干湿表常数 gamma (kPa/°C) Args: pressure (float): 大气压 (kPa) Cp 1.013e-3 # 空气定压比热 (MJ/kg/°C) epsilon 0.622 # 水汽与干空气分子量之比 lambda_val 2.45 # 汽化潜热 (MJ/kg) gamma (Cp * pressure) / (epsilon * lambda_val) return gamma def et0_priestley_taylor(tmean, rn, pressure, alpha1.26, g0): 计算普里斯特利-泰勒潜在蒸散发 ET0 (mm/day) Args: tmean (float): 日均温 (°C) rn (float): 地表净辐射 (MJ/m²/day) pressure (float): 大气压 (kPa) alpha (float): 系数默认1.26 g (float): 土壤热通量 (MJ/m²/day)日尺度常取0 Returns: float: ET0 (mm/day) delta calculate_delta(tmean) gamma calculate_gamma(pressure) lambda_val 2.45 # 公式计算结果单位转换 (MJ/m²/day) / (MJ/kg) * 1000 mm/day et0 alpha * (delta / (delta gamma)) * (rn - g) / lambda_val return et0 # 应用到DataFrame # 假设我们有净辐射列 Rn_MJ 和大气压列 Pressure_kPa (可通过海拔估算) df[Tmean] (df[Tmax_C] df[Tmin_C]) / 2.0 # 估算大气压简化公式海拔单位米 df[Pressure_kPa_est] 101.3 * ((293 - 0.0065 * df[Elevation_m]) / 293) ** 5.26 df[ET0_PT] df.apply( lambda row: et0_priestley_taylor(row[Tmean], row[Rn_MJ], row[Pressure_kPa_est]), axis1 )注意净辐射Rn的计算本身就是一个复杂过程需要太阳辐射、反射率、长波辐射等数据。如果你的数据只有日照时数需要先通过安格斯-普雷斯科特等公式估算太阳辐射再计算净辐射。这往往是误差的主要来源。4.3 彭曼-蒙蒂斯公式实现挑战与细节这是最复杂的一个。我们将严格按照FAO-56手册的步骤来实现。公式如下ET0 (0.408 * delta * (Rn - G) gamma * (900/(T273)) * u2 * (es - ea)) / (delta gamma * (1 0.34 * u2))其中新增了u2: 2米高处的风速m/s。es: 饱和水汽压kPa。ea: 实际水汽压kPa。T: 日均温°C。def calculate_es(tmean): 计算饱和水汽压 es (kPa) es 0.6108 * np.exp((17.27 * tmean) / (tmean 237.3)) return es def calculate_ea(tmean, rh_mean): 根据平均相对湿度计算实际水汽压 ea (kPa) es calculate_es(tmean) ea es * (rh_mean / 100.0) return ea def calculate_rn_from_sunshine(tmax, tmin, sunshine_hours, lat, doy, elevation): 一个简化的净辐射估算示例基于日照时数。 实际应用应使用更可靠的辐射数据或完整公式。 这里仅作演示计算净短波辐射 Rns。 # 1. 计算地外辐射 Ra (复用之前的函数) Ra calculate_ra(doy, lat) # 2. 根据日照时数估算太阳辐射 Rs (FAO-56 公式) # 假设最大可能日照时数 N 可以粗略估算这里简化 # 实际应使用更精确的日照时角计算N N 2 * np.arccos(-np.tan(lat*pi/180) * np.tan(calculate_solar_declination(doy))) * 24 / (2*pi) N np.maximum(N, 0.1) # 避免除零 Rs (0.25 0.5 * (sunshine_hours / N)) * Ra # 简化公式 # 3. 计算净短波辐射 Rns (假设反照率0.23适用于参考作物) albedo 0.23 Rns (1 - albedo) * Rs # 4. 净长波辐射 Rnl 计算非常复杂依赖温度、水汽等此处大幅简化 # 仅作示意实际项目务必使用完整公式 Rnl 0.0 # 简化假设 # 5. 净辐射 Rn Rns - Rnl Rn Rns - Rnl return Rn def calculate_solar_declination(doy): 计算太阳磁偏角 delta (弧度) theta 2 * pi * doy / 365.0 delta 0.409 * np.sin(theta - 1.39) return delta def et0_fao56(tmax, tmin, rh_mean, wind_speed, sunshine_hours, lat, doy, elevation): FAO-56 彭曼-蒙蒂斯公式计算 ET0 (mm/day) 这是一个简化版本净辐射计算不完整仅用于演示流程。 tmean (tmax tmin) / 2.0 # 1. 计算饱和水汽压 es 和实际水汽压 ea es calculate_es(tmean) ea calculate_ea(tmean, rh_mean) # 2. 计算饱和水汽压曲线斜率 delta delta calculate_delta(tmean) # 3. 计算干湿表常数 gamma # 估算大气压 pressure 101.3 * ((293 - 0.0065 * elevation) / 293) ** 5.26 gamma calculate_gamma(pressure) # 4. 估算净辐射 Rn (这里调用简化函数实际应用需替换) Rn calculate_rn_from_sunshine(tmax, tmin, sunshine_hours, lat, doy, elevation) G 0 # 日尺度土壤热通量 # 5. 空气动力项计算 # 确保风速是2米高处如果不是需要转换 u2 wind_speed # 假设数据已是2米风速 numerator_wind gamma * (900 / (tmean 273)) * u2 * (es - ea) # 6. 能量项计算 numerator_energy 0.408 * delta * (Rn - G) # 7. 分母计算 denominator delta gamma * (1 0.34 * u2) # 8. 计算 ET0 et0 (numerator_energy numerator_wind) / denominator return et0 # 应用到DataFrame df[ET0_FAO56] df.apply( lambda row: et0_fao56( row[Tmax_C], row[Tmin_C], row[RH_mean], row[WindSpeed_2m], row[Sunshine_hours], row[Latitude], row[DOY], row[Elevation_m] ), axis1 )踩坑实录彭曼-蒙蒂斯公式的实现90%的问题出在净辐射Rn的计算上。上面的calculate_rn_from_sunshine函数是一个极度简化的示例绝对不能用于严肃的科研或业务计算。完整的净辐射计算需要分别计算入射短波辐射、反射短波辐射、入射长波辐射和射出长波辐射其中涉及云量、水汽压、 Stefan-Boltzmann常数等。在实际项目中我强烈建议直接使用可靠的辐射观测数据。如果必须估算使用FAO-56或ASCE标准中完整的净辐射计算模块。考虑使用成熟的第三方库如Python的pyet或refet库它们已经实现了经过严格测试的FAO-56公式。另一个常见错误是单位。确保风速是m/s温度是°C辐射是MJ/m²/day压力是kPa。一个单位错误会导致结果差一个数量级。5. 结果分析与可视化让数据自己说话计算完成后我们得到了三列或更多潜在蒸散发数据。如何验证它们的合理性并从中提取信息5.1 数据合理性检查首先进行简单的统计和逻辑检查# 查看基本统计信息 print(df[[ET0_Hargreaves, ET0_PT, ET0_FAO56]].describe()) # 检查是否存在负值或异常大的值ET0一般0-15 mm/day print((df[[ET0_Hargreaves, ET0_PT, ET0_FAO56]] 0).sum()) print((df[[ET0_Hargreaves, ET0_PT, ET0_FAO56]] 20).sum())负值通常不合理可能是计算错误如辐射为负且绝对值过大或输入数据异常如温差为负。过大值在极端炎热干燥多风的天气可能出现但超过20 mm/day需要谨慎核查输入数据特别是辐射和风速。季节性ET0应有明显的季节变化夏季高冬季低。如果曲线平坦可能有问题。5.2 时间序列可视化绘制全年的ET0变化曲线对比不同方法的结果。import matplotlib.pyplot as plt import matplotlib.dates as mdates plt.figure(figsize(14, 7)) plt.plot(df.index, df[ET0_Hargreaves], labelHargreaves, alpha0.7, linewidth1) plt.plot(df.index, df[ET0_PT], labelPriestley-Taylor, alpha0.7, linewidth1) plt.plot(df.index, df[ET0_FAO56], labelFAO-56 Penman-Monteith, alpha0.7, linewidth1) plt.xlabel(Date) plt.ylabel(Potential Evapotranspiration (mm/day)) plt.title(Daily Potential Evapotranspiration Calculated by Different Methods) plt.legend() plt.grid(True, alpha0.3) # 优化x轴日期显示 plt.gca().xaxis.set_major_formatter(mdates.DateFormatter(%Y-%m)) plt.gca().xaxis.set_major_locator(mdates.MonthLocator()) plt.gcf().autofmt_xdate() plt.tight_layout() plt.show()通过看图你可以直观地判断趋势是否一致三种方法应该表现出相似的年变化趋势。量级差异哈格里夫斯和彭曼-蒙蒂斯结果可能接近普里斯特利-泰勒在干旱季节可能偏低因为它忽略了空气动力项。异常波动某一天某个方法出现尖峰或低谷需要回去检查那天的原始数据如异常高温、大风或辐射数据。5.3 散点图与相关性分析定量比较不同方法之间的一致性。fig, axes plt.subplots(1, 2, figsize(12, 5)) # Hargreaves vs FAO-56 ax1 axes[0] ax1.scatter(df[ET0_FAO56], df[ET0_Hargreaves], alpha0.5, s10) # 添加1:1线 lims [0, max(df[ET0_FAO56].max(), df[ET0_Hargreaves].max())] ax1.plot(lims, lims, k--, alpha0.75, label1:1 Line) ax1.set_xlabel(FAO-56 PM ET0 (mm/day)) ax1.set_ylabel(Hargreaves ET0 (mm/day)) ax1.set_title(Hargreaves vs FAO-56 PM) ax1.legend() ax1.grid(True, alpha0.3) # Priestley-Taylor vs FAO-56 ax2 axes[1] ax2.scatter(df[ET0_FAO56], df[ET0_PT], alpha0.5, s10, colororange) ax2.plot(lims, lims, k--, alpha0.75, label1:1 Line) ax2.set_xlabel(FAO-56 PM ET0 (mm/day)) ax2.set_ylabel(Priestley-Taylor ET0 (mm/day)) ax2.set_title(Priestley-Taylor vs FAO-56 PM) ax2.legend() ax2.grid(True, alpha0.3) plt.tight_layout() plt.show() # 计算相关系数 corr_h_vs_pm df[ET0_Hargreaves].corr(df[ET0_FAO56]) corr_pt_vs_pm df[ET0_PT].corr(df[ET0_FAO56]) print(fCorrelation (Hargreaves vs PM): {corr_h_vs_pm:.3f}) print(fCorrelation (PT vs PM): {corr_pt_vs_pm:.3f})如果散点紧密分布在1:1线两侧说明两种方法一致性高。如果出现系统性的偏离如PT法点全部在1:1线下方则说明该方法在你的研究区存在系统偏差可能需要校准。6. 性能优化与工程化思考从脚本到工具当数据量很大比如全国站点数十年逐日数据时直接使用DataFrame.apply逐行计算可能会比较慢。我们可以利用numpy的向量化运算进行优化。6.1 向量化计算改造以哈格里夫斯公式为例我们可以重写函数使其直接接受数组输入def et0_hargreaves_vectorized(tmax_arr, tmin_arr, lat, doy_arr): 向量化版本的哈格里夫斯公式计算 Args: tmax_arr, tmin_arr, doy_arr: 一维numpy数组 lat: 标量纬度 tmean_arr (tmax_arr tmin_arr) / 2.0 # 计算Ra也需要向量化这里假设doy_arr是数组 # 我们需要一个向量化的calculate_ra def calculate_ra_vectorized(doy_arr, lat): lat_rad lat * np.pi / 180.0 theta 2 * np.pi * doy_arr / 365.0 dr 1 0.033 * np.cos(theta) delta 0.409 * np.sin(theta - 1.39) tan_term -np.tan(lat_rad) * np.tan(delta) tan_term np.clip(tan_term, -1.0, 1.0) omega_s np.arccos(tan_term) Gsc 0.0820 day_minutes 24 * 60 Ra (day_minutes / np.pi) * Gsc * dr * ( omega_s * np.sin(lat_rad) * np.sin(delta) np.cos(lat_rad) * np.cos(delta) * np.sin(omega_s) ) return Ra Ra_arr calculate_ra_vectorized(doy_arr, lat) # 处理温差可能为负的情况 delta_t tmax_arr - tmin_arr delta_t np.where(delta_t 0, 0, delta_t) # 将负温差设为0 et0_arr 0.0023 * Ra_arr * (tmean_arr 17.8) * np.sqrt(delta_t) return et0_arr # 使用向量化函数速度大幅提升 tmax_values df[Tmax_C].to_numpy() tmin_values df[Tmin_C].to_numpy() doy_values df[DOY].to_numpy() df[ET0_Hargreaves_vec] et0_hargreaves_vectorized(tmax_values, tmin_values, lat, doy_values)对于彭曼-蒙蒂斯等复杂公式向量化能带来数量级的性能提升。6.2 模块化与封装为了让代码更易用、易维护我们可以将相关函数组织成一个模块.py文件。例如创建一个名为pet_calculator.py的文件# pet_calculator.py import numpy as np import pandas as pd class PETCalculator: 潜在蒸散发计算器 def __init__(self, latitude, elevation): self.lat latitude self.elev elevation def calculate_ra(self, doy): # ... 实现计算Ra的代码 pass def calculate_delta(self, tmean): # ... 实现计算delta的代码 pass def et0_hargreaves(self, df): # 接收DataFrame返回计算好的Series pass def et0_fao56(self, df, use_column_names{tmax:Tmax_C, ...}): # 通过字典映射列名增加灵活性 pass # 使用时 from pet_calculator import PETCalculator calc PETCalculator(latitude40.0, elevation50) df[ET0] calc.et0_fao56(df)这样主程序会变得非常简洁并且计算逻辑可以复用。6.3 处理大规模数据与并行计算对于全国站点数据可以按站点分组利用multiprocessing或joblib库进行并行计算。from multiprocessing import Pool import pandas as pd def calculate_site_et0(site_data): 处理单个站点数据的函数 site_id, df_site site_data # 假设df_site包含该站点所有数据 calc PETCalculator(latitudedf_site[Latitude].iloc[0], elevationdf_site[Elevation_m].iloc[0]) df_site[ET0] calc.et0_fao56(df_site) return site_id, df_site # 假设all_data是一个字典键为站点ID值为该站点的DataFrame all_data {...} with Pool(processes4) as pool: # 使用4个进程 results pool.map(calculate_site_et0, all_data.items()) # 将结果合并7. 常见问题排查与经验分享即使按照步骤操作你也可能会遇到一些奇怪的问题。这里分享几个我踩过的坑和解决办法。7.1 结果全是NaN或无穷大可能原因1数据中存在NaN或inf。在计算过程中如果输入数据包含NaN大部分numpy运算结果也会是NaN。使用df.isnull().sum()和np.isinf(df).sum()检查数据。可能原因2数学域错误。例如计算sqrt(tmax - tmin)时遇到负数或者计算arccos(x)时x不在[-1,1]区间内。务必在函数中加入np.clip进行数值保护。可能原因3单位错误导致数值极端。例如风速单位是km/h但被当作m/s会导致空气动力项巨大。仔细检查所有输入数据的单位。7.2 计算结果与已知文献或软件结果对不上第一步核对公式版本。彭曼-蒙蒂斯公式有FAO-56、ASCE等多个版本系数略有不同。确保你实现的公式与对比来源一致。第二步逐步验证中间变量。不要只比较最终ET0。将你的中间变量如Ra,delta,es,ea,Rn与可靠来源如FAO-56手册附录的算例进行对比。这是定位问题最有效的方法。第三步检查常数取值。例如干湿表常数中的汽化潜热lambda是2.45 MJ/kg还是2.5太阳常数用0.0820还是0.0864 MJ/m²/min这些细微差别都会影响结果。第四步时区与日界。你的数据日期是当地时间还是世界时日平均值是从当地0点到24点吗辐射数据是日总量吗时间不一致会导致与基于不同时间基准的计算结果产生偏差。7.3 季节性曲线出现不合理的“锯齿”或突变检查原始数据质量直接绘制Tmax,Tmin,Sunshine_hours等原始数据的时序图。突变往往源于原始数据的错误或缺失值插补不当。辐射计算是重灾区如果使用日照时数估算辐射那么Sunshine_hours数据的质量直接决定了Rn和最终ET0的平滑度。日照数据本身可能波动很大。风速的影响在干燥季节风速对彭曼-蒙蒂斯公式结果影响显著。一个异常的风速峰值会导致ET0出现一个尖峰。7.4 如何选择最终结果如果计算了多种方法该信哪个以数据最全的方法为基准。如果你有完整的辐射、湿度、风速数据那么彭曼-蒙蒂斯FAO-56的结果通常最可靠。考虑研究区域。在干旱区普里斯特利-泰勒公式可能系统性低估。哈格里夫斯公式在有的地区需要本地化校准。进行交叉验证。如果有可能找到你研究区域内其他已发表的ET0数据或可靠的模型输出与你的结果进行对比。敏感性分析。可以稍微改变某个输入参数如温度±1°C看ET0的变化是否在合理范围内。这有助于理解结果的不确定性。最后记住文档和注释的重要性。在你的代码中清晰地注明每个公式的来源如FAO-56, ASCE Standardized Reference Evapotranspiration Equation等记录下所有的单位转换和假设。这不仅是良好的编程习惯更是科学研究可复现性的基本要求。当你半年后回头看这段代码或者其他人要使用它时详细的注释能节省大量时间。本文还有配套的精品资源点击获取 Python实现潜在蒸散发计算：从彭曼公式到数据可视化
简介本资源是一套面向水文、气象及农业科研工作者的潜在蒸散发ET计算Python工具集覆盖25种主流经验与物理模型方法如Penman-Monteith、Hargreaves-Samani、Thornthwaite、Priestley-Taylor等可支撑流域水文模拟、灌溉需水量估算、气候响应分析等实际研究任务。压缩包共11个文件含5个可读可修改的.py源码模块如et_methods.py、utils.py、converter.py与6个已编译的.pyc文件总大小仅117KB轻量易集成适合Python中高级用户快速调用或二次开发。已有1014人学习下载资源结构清晰核心算法封装于et_methods模块单位换算与数据预处理由converter和utils支持global_variables统一管理参数init.py实现包级导入便于按需调用单个模型或批量对比不同方法结果。1. 从气象数据到代码实现为什么我们需要计算潜在蒸散发如果你在农业、水文、生态或者气候变化研究领域工作过那么“潜在蒸散发”这个词一定不陌生。它听起来有点学术但说白了就是在一个理想条件下一片完全被植被覆盖、水分供应充足的地面单位时间内能够蒸发和蒸腾到大气中的总水量。这个“理想条件”很关键它排除了土壤水分不足、植被类型差异等限制因素反映的是纯粹由气象条件比如太阳辐射、温度、湿度、风速决定的蒸发能力。为什么这个“理想值”如此重要因为它是一个基准。在实际应用中比如农田灌溉管理我们知道了潜在蒸散发量再结合土壤墒情和作物系数就能估算出作物实际需要多少水从而制定精准的灌溉计划避免水资源的浪费。在水文模型中它是计算流域水量平衡、预测径流的关键输入参数。在气候变化研究中潜在蒸散发的长期趋势分析是评估干旱风险、生态系统响应的重要指标。然而计算潜在蒸散发并不是一个简单的加减乘除。它背后是复杂的物理过程涉及能量平衡和空气动力学原理。历史上科学家们提出了多种经验或半经验公式来估算它比如彭曼公式、彭曼-蒙蒂斯公式、哈格里夫斯公式、索恩思韦特公式等。这些公式各有优劣有的需要的数据多但精度高有的数据要求低但适用区域有限。过去这些计算往往依赖于专业的商业软件或需要手动查表、套公式过程繁琐且容易出错。而现在Python以其强大的科学计算库和简洁的语法成为了处理这类问题的利器。我们可以用几行代码就自动化地完成从原始气象数据读取、质量检查、公式计算到结果可视化的全过程。这不仅大大提高了工作效率也让研究方法更加透明和可复现。今天我就结合自己处理农业气象数据的经验带你一步步用Python实现几种主流的潜在蒸散发计算方法并分享一些实操中容易踩的坑和优化技巧。2. 核心公式选型面对一堆气象数据我该用哪个公式当你拿到一组气象数据准备计算潜在蒸散发时第一个问题往往是该用哪个公式这不是拍脑袋决定的而是由你手头数据的完整性和精度要求共同决定的。选择不当要么巧妇难为无米之炊要么就是杀鸡用牛刀浪费了高质量数据。下面我梳理了几个最常用的公式及其数据需求你可以对号入座。2.1 彭曼-蒙蒂斯公式当之无愧的“金标准”如果你拥有相对完整的气象站数据那么彭曼-蒙蒂斯公式通常是首选。它被联合国粮农组织推荐是目前理论上最完备、应用最广泛的公式。它综合考虑了净辐射能量项和空气干燥度、风速空气动力项。所需核心数据日均气温最高温、最低温。日照时数或太阳辐射用于计算净辐射。相对湿度或露点温度反映空气湿度。风速通常指2米高处的风速。站点经纬度和海拔用于计算太阳常数、大气压等。为什么选它因为它物理基础扎实在全球多数地区表现稳定。但它的“娇贵”之处在于对数据质量要求高尤其是辐射数据。如果辐射数据缺失或不准计算结果误差会很大。2.2 哈格里夫斯公式数据匮乏时的“救星”在很多情况下尤其是历史数据或偏远地区我们可能只有温度数据。这时哈格里夫斯公式就派上大用场了。它只需要日均最高温、最低温以及站点纬度通过一个经验系数来估算太阳辐射进而计算潜在蒸散发。所需核心数据日均气温最高温、最低温。站点纬度用于估算地外辐射。它的优势与局限优势极其明显——数据需求极简。我在处理一些上世纪的气象数据时它几乎是唯一可行的选择。但它的精度通常低于彭曼-蒙蒂斯公式在非常潮湿或非常干燥的地区可能需要本地化校正系数。2.3 普里斯特利-泰勒公式湿润地区的简化方案这个公式基于能量平衡假设空气动力项可以用一个常数比例α通常取1.26与净辐射项关联。它适用于水分供应充足、大面积均匀的湿润表面比如茂密的森林或灌溉充分的农田。所需核心数据净辐射这是核心输入。气温用于计算饱和水汽压曲线斜率等热力学参数。适用场景当你关注的是能量限制为主的蒸发过程且风速、湿度数据不可靠时可以考虑它。但在干旱半干旱地区它会显著高估蒸散发。为了更直观地对比我整理了一个选型决策表公式名称核心数据需求计算复杂度适用场景主要局限彭曼-蒙蒂斯温度、湿度、风速、辐射、站点信息高数据齐全追求高精度全球多数地区数据要求高计算步骤多哈格里夫斯最高/最低温、纬度低只有温度数据历史数据或偏远地区精度相对较低需地区校准普里斯特利-泰勒净辐射、温度中湿润下垫面能量限制为主干旱区误差大依赖净辐射精度我的经验之谈在实际项目中我通常会做一个“数据审计”。先列出所有可用数据字段及其缺失率。如果辐射、风速数据质量尚可毫不犹豫用彭曼-蒙蒂斯。如果只有温度就用哈格里夫斯先跑出一个基准结果并在报告中明确说明其不确定性。有时我甚至会并行计算两种方法通过对比结果来交叉验证数据的合理性。3. 环境搭建与数据准备别在第一步就掉进坑里工欲善其事必先利其器。一个稳定、可复现的Python环境是后续所有工作的基础。很多人觉得安装包很简单但恰恰是这里隐藏着最多的版本冲突和依赖问题。3.1 创建独立的Python环境我强烈建议使用conda或venv为这个项目创建一个独立的虚拟环境。这能确保你的库版本不会干扰其他项目。# 使用 conda (假设你安装了Anaconda或Miniconda) conda create -n pet_calc python3.9 conda activate pet_calc # 或者使用 venv python -m venv pet_env # Windows pet_env\Scripts\activate # Linux/Mac source pet_env/bin/activate3.2 安装核心计算库在我们的环境中需要安装几个核心库pip install numpy pandas matplotlibnumpy: 数值计算的基石所有公式中的数组运算都靠它。pandas: 数据处理的瑞士军刀读取CSV/Excel、处理时间序列、数据清洗离不开它。matplotlib: 基础绘图库用于可视化结果和输入数据。一个关键的坑scipy的隐式依赖。在计算饱和水汽压、斜率等时我们可能会用到一些数学函数。虽然彭曼-蒙蒂斯公式可以手动实现但为了稳健和方便我们通常会用到scipy中的常量或优化函数。但请注意scipy在某些系统上安装可能因为编译依赖而失败。一个更轻量级的替代是使用metpy或pyet这类气象专用库它们封装了这些计算。这里我们先以纯手工计算为例确保通用性。如果需要可以后续安装pip install scipy3.3 气象数据的读取与清洗假设你有一份名为weather_data.csv的日尺度气象数据其格式可能如下DateTmax_CTmin_CRH_meanWindSpeed_2mSunshine_hoursLatitudeLongitudeElevation_m2023-07-0132.520.165.22.19.540.0116.5502023-07-0234.021.560.82.510.240.0116.550第一步用pandas加载数据import pandas as pd # 读取数据并指定日期列为索引 df pd.read_csv(weather_data.csv, parse_dates[Date], index_colDate) print(df.head()) print(df.info()) # 查看数据概要和缺失值第二步处理缺失值与异常值这是最耗时但也最重要的一步。气象数据常有缺失或明显错误如温度超过60°C。# 1. 简单查看缺失情况 print(df.isnull().sum()) # 2. 对于少量缺失可以用前后值插补需谨慎 # 例如用前一天的湿度填充当天的缺失值仅适用于连续缺失少的情况 df[RH_mean].fillna(methodffill, inplaceTrue) # 3. 对于明显异常值可以设定合理范围进行过滤或标记 # 例如假设研究地点在中国东部日最高温超过45°C或低于-20°C视为异常 df.loc[(df[Tmax_C] 45) | (df[Tmax_C] -20), Tmax_C] pd.NA df.loc[(df[Tmin_C] 35) | (df[Tmin_C] -30), Tmin_C] pd.NA # 4. 删除仍存在关键数据缺失的行例如温度数据缺失就无法计算 # 这里假设Tmax和Tmin是必须的 df.dropna(subset[Tmax_C, Tmin_C], inplaceTrue)注意插补方法需要根据数据特点选择。对于气象序列时间序列插值如线性插值、样条插值可能比简单的前向填充更合理。可以使用df.interpolate(methodtime)。但务必记录下你的处理步骤这在科研中至关重要。第三步单位检查与转换公式计算对单位非常敏感。确保你的数据单位与公式要求一致。常见转换包括温度公式通常使用摄氏度°C。如果你的数据是开尔文K需减去273.15。风速彭曼-蒙蒂斯公式通常需要2米高处的风速m/s。如果你的风速是10米高处或单位是km/h需要转换。辐射最易出错的地方。公式需要的是每日净辐射或太阳短波辐射单位是MJ/m²/day。如果你的日照时数单位是小时需要转换为辐射能量。4. 手把手实现三种主流公式的Python代码详解数据准备好了环境也搭好了现在让我们进入核心环节——编码实现。我将分别实现哈格里夫斯、普里斯特利-泰勒和彭曼-蒙蒂斯公式并解释每一步的物理意义和计算细节。4.1 哈格里夫斯公式实现简约而不简单哈格里夫斯公式的原始形式如下ET0 0.0023 * Ra * (Tmean 17.8) * sqrt(Tmax - Tmin)其中Ra是地外辐射MJ/m²/dayTmean是日均温Tmax和Tmin是日最高/最低温。这里的关键是计算Ra它取决于一年中的第几天DOY和站点纬度。import numpy as np import pandas as pd from math import pi, sin, cos, asin def calculate_ra(doy, lat_deg): 计算地外辐射 Ra (MJ/m²/day) Args: doy (int): 年积日1-365/366 lat_deg (float): 站点纬度度北纬为正南纬为负 Returns: float: 地外辐射 Ra # 1. 将纬度转换为弧度 lat_rad lat_deg * pi / 180.0 # 2. 计算日地相对距离倒数 dr 和太阳磁偏角 delta # 日角 theta 2 * pi * doy / 365.0 # 日地距离倒数 dr 1 0.033 * np.cos(theta) # 太阳磁偏角弧度 delta 0.409 * np.sin(theta - 1.39) # 3. 计算日落时角 omega_s (弧度) # 当 tan(lat)*tan(delta) -1 或 1 时表示极昼或极夜需要特殊处理 tan_term -np.tan(lat_rad) * np.tan(delta) # 限制值在[-1,1]之间避免数学错误 tan_term np.clip(tan_term, -1.0, 1.0) omega_s np.arccos(tan_term) # 4. 计算 Ra # 太阳常数 Gsc 0.0820 MJ/m²/min Gsc 0.0820 # 一天中的分钟数 day_minutes 24 * 60 Ra (day_minutes / pi) * Gsc * dr * ( omega_s * np.sin(lat_rad) * np.sin(delta) np.cos(lat_rad) * np.cos(delta) * np.sin(omega_s) ) return Ra def et0_hargreaves(tmax, tmin, lat, doy): 计算哈格里夫斯潜在蒸散发 ET0 (mm/day) Args: tmax (float): 日最高温 (°C) tmin (float): 日最低温 (°C) lat (float): 纬度 (度) doy (int): 年积日 Returns: float: ET0 (mm/day) tmean (tmax tmin) / 2.0 Ra calculate_ra(doy, lat) # 原始哈格里夫斯公式 et0 0.0023 * Ra * (tmean 17.8) * np.sqrt(tmax - tmin) return et0 # 应用到整个DataFrame # 首先我们需要为每一行计算年积日Day of Year df[DOY] df.index.dayofyear # 假设纬度存储在列Latitude中且所有行纬度相同 lat df[Latitude].iloc[0] # 使用apply函数逐行计算注意sqrt里温差可能为负需处理 df[ET0_Hargreaves] df.apply( lambda row: et0_hargreaves(row[Tmax_C], row[Tmin_C], lat, row[DOY]), axis1 )实操心得np.sqrt(tmax - tmin)这里有个隐患。在极少数情况下如某些海洋性气候tmax可能略低于tmin可能是数据错误或四舍五入导致导致对负数开方报错。一个稳健的做法是加上一个很小的数或取绝对值np.sqrt(np.abs(tmax - tmin) 1e-10)。另外FAO后来推荐了修正的哈格里夫斯公式系数不同如果需要更高精度可以查阅FAO-56手册。4.2 普里斯特利-泰勒公式实现抓住能量核心普里斯特利-泰勒公式ET0 alpha * (delta / (delta gamma)) * (Rn - G) / lambda其中alpha经验系数通常取1.26湿润地区。delta饱和水汽压曲线斜率kPa/°C。gamma干湿表常数kPa/°C。Rn地表净辐射MJ/m²/day。G土壤热通量MJ/m²/day对于日尺度通常近似为0。lambda水的汽化潜热约2.45 MJ/kg。def calculate_delta(tmean): 计算饱和水汽压曲线斜率 delta (kPa/°C) FAO-56 推荐公式 # 计算在温度Tmean下的饱和水汽压 e_s_t 0.6108 * np.exp((17.27 * tmean) / (tmean 237.3)) delta (4098 * e_s_t) / ((tmean 237.3) ** 2) return delta def calculate_gamma(pressure): 计算干湿表常数 gamma (kPa/°C) Args: pressure (float): 大气压 (kPa) Cp 1.013e-3 # 空气定压比热 (MJ/kg/°C) epsilon 0.622 # 水汽与干空气分子量之比 lambda_val 2.45 # 汽化潜热 (MJ/kg) gamma (Cp * pressure) / (epsilon * lambda_val) return gamma def et0_priestley_taylor(tmean, rn, pressure, alpha1.26, g0): 计算普里斯特利-泰勒潜在蒸散发 ET0 (mm/day) Args: tmean (float): 日均温 (°C) rn (float): 地表净辐射 (MJ/m²/day) pressure (float): 大气压 (kPa) alpha (float): 系数默认1.26 g (float): 土壤热通量 (MJ/m²/day)日尺度常取0 Returns: float: ET0 (mm/day) delta calculate_delta(tmean) gamma calculate_gamma(pressure) lambda_val 2.45 # 公式计算结果单位转换 (MJ/m²/day) / (MJ/kg) * 1000 mm/day et0 alpha * (delta / (delta gamma)) * (rn - g) / lambda_val return et0 # 应用到DataFrame # 假设我们有净辐射列 Rn_MJ 和大气压列 Pressure_kPa (可通过海拔估算) df[Tmean] (df[Tmax_C] df[Tmin_C]) / 2.0 # 估算大气压简化公式海拔单位米 df[Pressure_kPa_est] 101.3 * ((293 - 0.0065 * df[Elevation_m]) / 293) ** 5.26 df[ET0_PT] df.apply( lambda row: et0_priestley_taylor(row[Tmean], row[Rn_MJ], row[Pressure_kPa_est]), axis1 )注意净辐射Rn的计算本身就是一个复杂过程需要太阳辐射、反射率、长波辐射等数据。如果你的数据只有日照时数需要先通过安格斯-普雷斯科特等公式估算太阳辐射再计算净辐射。这往往是误差的主要来源。4.3 彭曼-蒙蒂斯公式实现挑战与细节这是最复杂的一个。我们将严格按照FAO-56手册的步骤来实现。公式如下ET0 (0.408 * delta * (Rn - G) gamma * (900/(T273)) * u2 * (es - ea)) / (delta gamma * (1 0.34 * u2))其中新增了u2: 2米高处的风速m/s。es: 饱和水汽压kPa。ea: 实际水汽压kPa。T: 日均温°C。def calculate_es(tmean): 计算饱和水汽压 es (kPa) es 0.6108 * np.exp((17.27 * tmean) / (tmean 237.3)) return es def calculate_ea(tmean, rh_mean): 根据平均相对湿度计算实际水汽压 ea (kPa) es calculate_es(tmean) ea es * (rh_mean / 100.0) return ea def calculate_rn_from_sunshine(tmax, tmin, sunshine_hours, lat, doy, elevation): 一个简化的净辐射估算示例基于日照时数。 实际应用应使用更可靠的辐射数据或完整公式。 这里仅作演示计算净短波辐射 Rns。 # 1. 计算地外辐射 Ra (复用之前的函数) Ra calculate_ra(doy, lat) # 2. 根据日照时数估算太阳辐射 Rs (FAO-56 公式) # 假设最大可能日照时数 N 可以粗略估算这里简化 # 实际应使用更精确的日照时角计算N N 2 * np.arccos(-np.tan(lat*pi/180) * np.tan(calculate_solar_declination(doy))) * 24 / (2*pi) N np.maximum(N, 0.1) # 避免除零 Rs (0.25 0.5 * (sunshine_hours / N)) * Ra # 简化公式 # 3. 计算净短波辐射 Rns (假设反照率0.23适用于参考作物) albedo 0.23 Rns (1 - albedo) * Rs # 4. 净长波辐射 Rnl 计算非常复杂依赖温度、水汽等此处大幅简化 # 仅作示意实际项目务必使用完整公式 Rnl 0.0 # 简化假设 # 5. 净辐射 Rn Rns - Rnl Rn Rns - Rnl return Rn def calculate_solar_declination(doy): 计算太阳磁偏角 delta (弧度) theta 2 * pi * doy / 365.0 delta 0.409 * np.sin(theta - 1.39) return delta def et0_fao56(tmax, tmin, rh_mean, wind_speed, sunshine_hours, lat, doy, elevation): FAO-56 彭曼-蒙蒂斯公式计算 ET0 (mm/day) 这是一个简化版本净辐射计算不完整仅用于演示流程。 tmean (tmax tmin) / 2.0 # 1. 计算饱和水汽压 es 和实际水汽压 ea es calculate_es(tmean) ea calculate_ea(tmean, rh_mean) # 2. 计算饱和水汽压曲线斜率 delta delta calculate_delta(tmean) # 3. 计算干湿表常数 gamma # 估算大气压 pressure 101.3 * ((293 - 0.0065 * elevation) / 293) ** 5.26 gamma calculate_gamma(pressure) # 4. 估算净辐射 Rn (这里调用简化函数实际应用需替换) Rn calculate_rn_from_sunshine(tmax, tmin, sunshine_hours, lat, doy, elevation) G 0 # 日尺度土壤热通量 # 5. 空气动力项计算 # 确保风速是2米高处如果不是需要转换 u2 wind_speed # 假设数据已是2米风速 numerator_wind gamma * (900 / (tmean 273)) * u2 * (es - ea) # 6. 能量项计算 numerator_energy 0.408 * delta * (Rn - G) # 7. 分母计算 denominator delta gamma * (1 0.34 * u2) # 8. 计算 ET0 et0 (numerator_energy numerator_wind) / denominator return et0 # 应用到DataFrame df[ET0_FAO56] df.apply( lambda row: et0_fao56( row[Tmax_C], row[Tmin_C], row[RH_mean], row[WindSpeed_2m], row[Sunshine_hours], row[Latitude], row[DOY], row[Elevation_m] ), axis1 )踩坑实录彭曼-蒙蒂斯公式的实现90%的问题出在净辐射Rn的计算上。上面的calculate_rn_from_sunshine函数是一个极度简化的示例绝对不能用于严肃的科研或业务计算。完整的净辐射计算需要分别计算入射短波辐射、反射短波辐射、入射长波辐射和射出长波辐射其中涉及云量、水汽压、 Stefan-Boltzmann常数等。在实际项目中我强烈建议直接使用可靠的辐射观测数据。如果必须估算使用FAO-56或ASCE标准中完整的净辐射计算模块。考虑使用成熟的第三方库如Python的pyet或refet库它们已经实现了经过严格测试的FAO-56公式。另一个常见错误是单位。确保风速是m/s温度是°C辐射是MJ/m²/day压力是kPa。一个单位错误会导致结果差一个数量级。5. 结果分析与可视化让数据自己说话计算完成后我们得到了三列或更多潜在蒸散发数据。如何验证它们的合理性并从中提取信息5.1 数据合理性检查首先进行简单的统计和逻辑检查# 查看基本统计信息 print(df[[ET0_Hargreaves, ET0_PT, ET0_FAO56]].describe()) # 检查是否存在负值或异常大的值ET0一般0-15 mm/day print((df[[ET0_Hargreaves, ET0_PT, ET0_FAO56]] 0).sum()) print((df[[ET0_Hargreaves, ET0_PT, ET0_FAO56]] 20).sum())负值通常不合理可能是计算错误如辐射为负且绝对值过大或输入数据异常如温差为负。过大值在极端炎热干燥多风的天气可能出现但超过20 mm/day需要谨慎核查输入数据特别是辐射和风速。季节性ET0应有明显的季节变化夏季高冬季低。如果曲线平坦可能有问题。5.2 时间序列可视化绘制全年的ET0变化曲线对比不同方法的结果。import matplotlib.pyplot as plt import matplotlib.dates as mdates plt.figure(figsize(14, 7)) plt.plot(df.index, df[ET0_Hargreaves], labelHargreaves, alpha0.7, linewidth1) plt.plot(df.index, df[ET0_PT], labelPriestley-Taylor, alpha0.7, linewidth1) plt.plot(df.index, df[ET0_FAO56], labelFAO-56 Penman-Monteith, alpha0.7, linewidth1) plt.xlabel(Date) plt.ylabel(Potential Evapotranspiration (mm/day)) plt.title(Daily Potential Evapotranspiration Calculated by Different Methods) plt.legend() plt.grid(True, alpha0.3) # 优化x轴日期显示 plt.gca().xaxis.set_major_formatter(mdates.DateFormatter(%Y-%m)) plt.gca().xaxis.set_major_locator(mdates.MonthLocator()) plt.gcf().autofmt_xdate() plt.tight_layout() plt.show()通过看图你可以直观地判断趋势是否一致三种方法应该表现出相似的年变化趋势。量级差异哈格里夫斯和彭曼-蒙蒂斯结果可能接近普里斯特利-泰勒在干旱季节可能偏低因为它忽略了空气动力项。异常波动某一天某个方法出现尖峰或低谷需要回去检查那天的原始数据如异常高温、大风或辐射数据。5.3 散点图与相关性分析定量比较不同方法之间的一致性。fig, axes plt.subplots(1, 2, figsize(12, 5)) # Hargreaves vs FAO-56 ax1 axes[0] ax1.scatter(df[ET0_FAO56], df[ET0_Hargreaves], alpha0.5, s10) # 添加1:1线 lims [0, max(df[ET0_FAO56].max(), df[ET0_Hargreaves].max())] ax1.plot(lims, lims, k--, alpha0.75, label1:1 Line) ax1.set_xlabel(FAO-56 PM ET0 (mm/day)) ax1.set_ylabel(Hargreaves ET0 (mm/day)) ax1.set_title(Hargreaves vs FAO-56 PM) ax1.legend() ax1.grid(True, alpha0.3) # Priestley-Taylor vs FAO-56 ax2 axes[1] ax2.scatter(df[ET0_FAO56], df[ET0_PT], alpha0.5, s10, colororange) ax2.plot(lims, lims, k--, alpha0.75, label1:1 Line) ax2.set_xlabel(FAO-56 PM ET0 (mm/day)) ax2.set_ylabel(Priestley-Taylor ET0 (mm/day)) ax2.set_title(Priestley-Taylor vs FAO-56 PM) ax2.legend() ax2.grid(True, alpha0.3) plt.tight_layout() plt.show() # 计算相关系数 corr_h_vs_pm df[ET0_Hargreaves].corr(df[ET0_FAO56]) corr_pt_vs_pm df[ET0_PT].corr(df[ET0_FAO56]) print(fCorrelation (Hargreaves vs PM): {corr_h_vs_pm:.3f}) print(fCorrelation (PT vs PM): {corr_pt_vs_pm:.3f})如果散点紧密分布在1:1线两侧说明两种方法一致性高。如果出现系统性的偏离如PT法点全部在1:1线下方则说明该方法在你的研究区存在系统偏差可能需要校准。6. 性能优化与工程化思考从脚本到工具当数据量很大比如全国站点数十年逐日数据时直接使用DataFrame.apply逐行计算可能会比较慢。我们可以利用numpy的向量化运算进行优化。6.1 向量化计算改造以哈格里夫斯公式为例我们可以重写函数使其直接接受数组输入def et0_hargreaves_vectorized(tmax_arr, tmin_arr, lat, doy_arr): 向量化版本的哈格里夫斯公式计算 Args: tmax_arr, tmin_arr, doy_arr: 一维numpy数组 lat: 标量纬度 tmean_arr (tmax_arr tmin_arr) / 2.0 # 计算Ra也需要向量化这里假设doy_arr是数组 # 我们需要一个向量化的calculate_ra def calculate_ra_vectorized(doy_arr, lat): lat_rad lat * np.pi / 180.0 theta 2 * np.pi * doy_arr / 365.0 dr 1 0.033 * np.cos(theta) delta 0.409 * np.sin(theta - 1.39) tan_term -np.tan(lat_rad) * np.tan(delta) tan_term np.clip(tan_term, -1.0, 1.0) omega_s np.arccos(tan_term) Gsc 0.0820 day_minutes 24 * 60 Ra (day_minutes / np.pi) * Gsc * dr * ( omega_s * np.sin(lat_rad) * np.sin(delta) np.cos(lat_rad) * np.cos(delta) * np.sin(omega_s) ) return Ra Ra_arr calculate_ra_vectorized(doy_arr, lat) # 处理温差可能为负的情况 delta_t tmax_arr - tmin_arr delta_t np.where(delta_t 0, 0, delta_t) # 将负温差设为0 et0_arr 0.0023 * Ra_arr * (tmean_arr 17.8) * np.sqrt(delta_t) return et0_arr # 使用向量化函数速度大幅提升 tmax_values df[Tmax_C].to_numpy() tmin_values df[Tmin_C].to_numpy() doy_values df[DOY].to_numpy() df[ET0_Hargreaves_vec] et0_hargreaves_vectorized(tmax_values, tmin_values, lat, doy_values)对于彭曼-蒙蒂斯等复杂公式向量化能带来数量级的性能提升。6.2 模块化与封装为了让代码更易用、易维护我们可以将相关函数组织成一个模块.py文件。例如创建一个名为pet_calculator.py的文件# pet_calculator.py import numpy as np import pandas as pd class PETCalculator: 潜在蒸散发计算器 def __init__(self, latitude, elevation): self.lat latitude self.elev elevation def calculate_ra(self, doy): # ... 实现计算Ra的代码 pass def calculate_delta(self, tmean): # ... 实现计算delta的代码 pass def et0_hargreaves(self, df): # 接收DataFrame返回计算好的Series pass def et0_fao56(self, df, use_column_names{tmax:Tmax_C, ...}): # 通过字典映射列名增加灵活性 pass # 使用时 from pet_calculator import PETCalculator calc PETCalculator(latitude40.0, elevation50) df[ET0] calc.et0_fao56(df)这样主程序会变得非常简洁并且计算逻辑可以复用。6.3 处理大规模数据与并行计算对于全国站点数据可以按站点分组利用multiprocessing或joblib库进行并行计算。from multiprocessing import Pool import pandas as pd def calculate_site_et0(site_data): 处理单个站点数据的函数 site_id, df_site site_data # 假设df_site包含该站点所有数据 calc PETCalculator(latitudedf_site[Latitude].iloc[0], elevationdf_site[Elevation_m].iloc[0]) df_site[ET0] calc.et0_fao56(df_site) return site_id, df_site # 假设all_data是一个字典键为站点ID值为该站点的DataFrame all_data {...} with Pool(processes4) as pool: # 使用4个进程 results pool.map(calculate_site_et0, all_data.items()) # 将结果合并7. 常见问题排查与经验分享即使按照步骤操作你也可能会遇到一些奇怪的问题。这里分享几个我踩过的坑和解决办法。7.1 结果全是NaN或无穷大可能原因1数据中存在NaN或inf。在计算过程中如果输入数据包含NaN大部分numpy运算结果也会是NaN。使用df.isnull().sum()和np.isinf(df).sum()检查数据。可能原因2数学域错误。例如计算sqrt(tmax - tmin)时遇到负数或者计算arccos(x)时x不在[-1,1]区间内。务必在函数中加入np.clip进行数值保护。可能原因3单位错误导致数值极端。例如风速单位是km/h但被当作m/s会导致空气动力项巨大。仔细检查所有输入数据的单位。7.2 计算结果与已知文献或软件结果对不上第一步核对公式版本。彭曼-蒙蒂斯公式有FAO-56、ASCE等多个版本系数略有不同。确保你实现的公式与对比来源一致。第二步逐步验证中间变量。不要只比较最终ET0。将你的中间变量如Ra,delta,es,ea,Rn与可靠来源如FAO-56手册附录的算例进行对比。这是定位问题最有效的方法。第三步检查常数取值。例如干湿表常数中的汽化潜热lambda是2.45 MJ/kg还是2.5太阳常数用0.0820还是0.0864 MJ/m²/min这些细微差别都会影响结果。第四步时区与日界。你的数据日期是当地时间还是世界时日平均值是从当地0点到24点吗辐射数据是日总量吗时间不一致会导致与基于不同时间基准的计算结果产生偏差。7.3 季节性曲线出现不合理的“锯齿”或突变检查原始数据质量直接绘制Tmax,Tmin,Sunshine_hours等原始数据的时序图。突变往往源于原始数据的错误或缺失值插补不当。辐射计算是重灾区如果使用日照时数估算辐射那么Sunshine_hours数据的质量直接决定了Rn和最终ET0的平滑度。日照数据本身可能波动很大。风速的影响在干燥季节风速对彭曼-蒙蒂斯公式结果影响显著。一个异常的风速峰值会导致ET0出现一个尖峰。7.4 如何选择最终结果如果计算了多种方法该信哪个以数据最全的方法为基准。如果你有完整的辐射、湿度、风速数据那么彭曼-蒙蒂斯FAO-56的结果通常最可靠。考虑研究区域。在干旱区普里斯特利-泰勒公式可能系统性低估。哈格里夫斯公式在有的地区需要本地化校准。进行交叉验证。如果有可能找到你研究区域内其他已发表的ET0数据或可靠的模型输出与你的结果进行对比。敏感性分析。可以稍微改变某个输入参数如温度±1°C看ET0的变化是否在合理范围内。这有助于理解结果的不确定性。最后记住文档和注释的重要性。在你的代码中清晰地注明每个公式的来源如FAO-56, ASCE Standardized Reference Evapotranspiration Equation等记录下所有的单位转换和假设。这不仅是良好的编程习惯更是科学研究可复现性的基本要求。当你半年后回头看这段代码或者其他人要使用它时详细的注释能节省大量时间。本文还有配套的精品资源点击获取