使用R中ET.PenmanMonteith函数计算FAO56参考作物ET时全为NA的问题排查

问题概述 #

使用R语言Evapotranspiration包的ET.PenmanMonteith函数计算FAO56参考作物蒸散量时，所有ET结果均为NA，统计值异常（均值NaN、最大值-Inf、最小值Inf），同时触发max/min函数警告（无非缺失参数）。用户已确认原始数据与格式化数据无缺失值，日期列格式正确。

运行输出 #

Penman-Monteith FAO56 Reference Crop ET
Evaporative surface: FAO-56 hypothetical short grass, albedo = 0.23 ; surface resistance = 70 sm^-1; crop height = 0.12  m; roughness height = 0.02 m
Sunshine hour data have been used for calculating incoming solar radiation
Wind data have been used for calculating the reference crop evapotranspiration
Timestep: daily
Units: mm
Time duration: 2022-01-01 to 2023-05-17
501 ET estimates obtained; 501 NA output entries due to missing data
Basic stats (NA excluded)
Mean: NaN
Max: -Inf
Min: Inf
Warning messages:
1: In max(res_ts, na.rm = T) :
  no non-missing arguments to max; returning -Inf
2: In min(res_ts, na.rm = T) :
  no non-missing arguments to min; returning Inf

核心代码 #

# 读取数据
data2 <- read.csv("climate_data_2.csv", head = TRUE, sep=",")

# 转换日期格式
data2$Date.daily <- as.Date(data2$Date.daily, format = "%Y-%m-%d")

# 创建数据列表
data <- list(Date.daily = data2$Date.daily,
             Date.monthly = format(data2$Date.daily, "%b %Y"))

# 向列表添加变量（转为zoo对象）
variable_names <- c("Tmax", "Tmin", "RHmax", "RHmin", "n", "uz")
for (var_name in variable_names) {
  data[[var_name]] <- zoo::zoo(data2[[var_name]], order.by = data2$Date.daily)
}

# 调用ET.PenmanMonteith函数
results <- ET.PenmanMonteith(data, constants, ts = "daily", solar = "sunshine hours", wind = "yes",
                             crop = "short", message = "yes", AdditionalStats = "yes", save.csv = "no")

原因分析 #

1. constants参数缺失关键地理信息：FAO56公式依赖站点纬度、海拔等参数计算辐射、气压项，若constants未定义lat（纬度）、elev（海拔）、z（风速测量高度）等核心元素，函数无法完成核心计算，直接返回全NA。
2. 输入变量数值逻辑异常：即使无缺失值，若存在Tmax<Tmin、RHmax<RHmin、日照时数n为负或超过当地理论最大日照时数等情况，函数内部校验会返回NA。
3. zoo对象时间索引不匹配：手动转为zoo对象时，若各变量的时间索引与Date.daily存在隐性差异（如时区、格式细微问题），函数无法正确匹配每日数据计算。
4. 参数配置不兼容：solar = "sunshine hours"模式下，函数需要计算理论最大日照时数，若缺少纬度参数或日期处理错误，辐射项计算失败导致全NA。

解决建议 #

1. 完善constants参数 #

确保constants包含FAO56公式必需的参数，示例如下：

constants <- list(
  lat = 30.5,         # 站点纬度（北纬正，南纬负，单位：度）
  elev = 150.0,       # 站点海拔（单位：米）
  z = 2.0,            # 风速测量高度（需与uz的测量高度一致，默认2m）
  lambda = 2.45,      # 汽化潜热（默认2.45 MJ/kg）
  Gsc = 0.0820        # 太阳常数（默认0.0820 MJ/m²/min）
)

• 纬度和海拔是必须项，缺失会直接导致辐射、气压计算失败；
• 风速测量高度z需与输入的uz（风速）的实际测量高度匹配，否则无法转换为FAO56要求的2m标准高度风速。

2. 校验输入变量的合理性 #

排查数值逻辑错误：

# 检查Tmax是否小于Tmin
sum(data2$Tmax < data2$Tmin, na.rm = TRUE)
# 检查RHmax是否小于RHmin
sum(data2$RHmax < data2$RHmin, na.rm = TRUE)
# 检查日照时数是否为负
sum(data2$n < 0, na.rm = TRUE)

若存在异常值，需修正原始数据（替换错误值或剔除异常记录）。

3. 简化数据输入格式 #

放弃手动转为zoo对象，直接使用向量输入（函数支持向量格式，兼容性更好）：

# 重构data列表，使用原始向量
data <- list(
  Date.daily = data2$Date.daily,
  Date.monthly = format(data2$Date.daily, "%b %Y"),
  Tmax = data2$Tmax,
  Tmin = data2$Tmin,
  RHmax = data2$RHmax,
  RHmin = data2$RHmin,
  n = data2$n,
  uz = data2$uz
)

4. 分步调试定位问题 #

若上述方法无效，手动计算FAO56公式的中间变量，确认哪个环节出错：

• 计算理论最大日照时数N；
• 计算净辐射Rn；
• 计算饱和水汽压亏缺es - ea；
• 计算空气动力学阻力ra。

通过分步验证，可精准定位导致全NA的计算环节。

内容的提问来源于stack exchange，提问作者Sondergaard