企业官网建设流程全解析

简介：本资源是一套面向水文、气象及农业科研工作者的潜在蒸散发（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 python=3.9 conda activate pet_calc # 或者使用 venv python -m venv pet_env # Windows pet_env\Scripts\activate # Linux/Mac source pet_env/bin/activate

3.2 安装核心计算库

在我们的环境中，需要安装几个核心库：

pip install numpy pandas matplotlib

• numpy: 数值计算的基石，所有公式中的数组运算都靠它。
• pandas: 数据处理的瑞士军刀，读取CSV/Excel、处理时间序列、数据清洗离不开它。
• matplotlib: 基础绘图库，用于可视化结果和输入数据。

一个关键的坑：scipy的隐式依赖。在计算饱和水汽压、斜率等时，我们可能会用到一些数学函数。虽然彭曼-蒙蒂斯公式可以手动实现，但为了稳健和方便，我们通常会用到scipy中的常量或优化函数。但请注意，scipy在某些系统上安装可能因为编译依赖而失败。一个更轻量级的替代是使用metpy或pyet这类气象专用库，它们封装了这些计算。这里我们先以纯手工计算为例，确保通用性。如果需要，可以后续安装：

pip install scipy

3.3 气象数据的读取与清洗

假设你有一份名为weather_data.csv的日尺度气象数据，其格式可能如下：

第一步：用pandas加载数据

import pandas as pd # 读取数据，并指定日期列为索引 df = pd.read_csv('weather_data.csv', parse_dates=['Date'], index_col='Date') print(df.head()) print(df.info()) # 查看数据概要和缺失值

第二步：处理缺失值与异常值这是最耗时但也最重要的一步。气象数据常有缺失或明显错误（如温度超过60°C）。

# 1. 简单查看缺失情况 print(df.isnull().sum()) # 2. 对于少量缺失，可以用前后值插补（需谨慎） # 例如，用前一天的湿度填充当天的缺失值（仅适用于连续缺失少的情况） df['RH_mean'].fillna(method='ffill', inplace=True) # 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'], inplace=True)

注意：插补方法需要根据数据特点选择。对于气象序列，时间序列插值（如线性插值、样条插值）可能比简单的前向填充更合理。可以使用df.interpolate(method='time')。但务必记录下你的处理步骤，这在科研中至关重要。

第三步：单位检查与转换公式计算对单位非常敏感。确保你的数据单位与公式要求一致。常见转换包括：

• 温度：公式通常使用摄氏度（°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²/day），Tmean是日均温，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']), axis=1 )

实操心得：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, alpha=1.26, g=0): """ 计算普里斯特利-泰勒潜在蒸散发 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']), axis=1 )

注意：净辐射Rn的计算本身就是一个复杂过程，需要太阳辐射、反射率、长波辐射等数据。如果你的数据只有日照时数，需要先通过安格斯-普雷斯科特等公式估算太阳辐射，再计算净辐射。这往往是误差的主要来源。

4.3 彭曼-蒙蒂斯公式实现：挑战与细节

这是最复杂的一个。我们将严格按照FAO-56手册的步骤来实现。公式如下：ET0 = (0.408 * delta * (Rn - G) + gamma * (900/(T+273)) * 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'] ), axis=1 )

踩坑实录：彭曼-蒙蒂斯公式的实现，90%的问题出在净辐射Rn的计算上。上面的calculate_rn_from_sunshine函数是一个极度简化的示例，绝对不能用于严肃的科研或业务计算。完整的净辐射计算需要分别计算入射短波辐射、反射短波辐射、入射长波辐射和射出长波辐射，其中涉及云量、水汽压、 Stefan-Boltzmann常数等。在实际项目中，我强烈建议：

1. 直接使用可靠的辐射观测数据。
2. 如果必须估算，使用FAO-56或ASCE标准中完整的净辐射计算模块。
3. 考虑使用成熟的第三方库，如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'], label='Hargreaves', alpha=0.7, linewidth=1) plt.plot(df.index, df['ET0_PT'], label='Priestley-Taylor', alpha=0.7, linewidth=1) plt.plot(df.index, df['ET0_FAO56'], label='FAO-56 Penman-Monteith', alpha=0.7, linewidth=1) plt.xlabel('Date') plt.ylabel('Potential Evapotranspiration (mm/day)') plt.title('Daily Potential Evapotranspiration Calculated by Different Methods') plt.legend() plt.grid(True, alpha=0.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()

通过看图，你可以直观地判断：

1. 趋势是否一致：三种方法应该表现出相似的年变化趋势。
2. 量级差异：哈格里夫斯和彭曼-蒙蒂斯结果可能接近，普里斯特利-泰勒在干旱季节可能偏低（因为它忽略了空气动力项）。
3. 异常波动：某一天某个方法出现尖峰或低谷，需要回去检查那天的原始数据（如异常高温、大风或辐射数据）。

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'], alpha=0.5, s=10) # 添加1:1线 lims = [0, max(df['ET0_FAO56'].max(), df['ET0_Hargreaves'].max())] ax1.plot(lims, lims, 'k--', alpha=0.75, label='1: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, alpha=0.3) # Priestley-Taylor vs FAO-56 ax2 = axes[1] ax2.scatter(df['ET0_FAO56'], df['ET0_PT'], alpha=0.5, s=10, color='orange') ax2.plot(lims, lims, 'k--', alpha=0.75, label='1: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, alpha=0.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(f"Correlation (Hargreaves vs PM): {corr_h_vs_pm:.3f}") print(f"Correlation (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(latitude=40.0, elevation=50) 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(latitude=df_site['Latitude'].iloc[0], elevation=df_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(processes=4) 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 如何选择最终结果？

如果计算了多种方法，该信哪个？

1. 以数据最全的方法为基准。如果你有完整的辐射、湿度、风速数据，那么彭曼-蒙蒂斯（FAO-56）的结果通常最可靠。
2. 考虑研究区域。在干旱区，普里斯特利-泰勒公式可能系统性低估。哈格里夫斯公式在有的地区需要本地化校准。
3. 进行交叉验证。如果有可能，找到你研究区域内其他已发表的ET0数据或可靠的模型输出，与你的结果进行对比。
4. 敏感性分析。可以稍微改变某个输入参数（如温度±1°C），看ET0的变化是否在合理范围内。这有助于理解结果的不确定性。

最后，记住文档和注释的重要性。在你的代码中，清晰地注明每个公式的来源（如FAO-56, ASCE Standardized Reference Evapotranspiration Equation等），记录下所有的单位转换和假设。这不仅是良好的编程习惯，更是科学研究可复现性的基本要求。当你半年后回头看这段代码，或者其他人要使用它时，详细的注释能节省大量时间。

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