使用HeatStress包Liljegren法计算WBGT时闰年报错问题

问题场景 #

使用HeatStress包的wbgt.Liljegren方法计算湿球黑球温度（WBGT），所有NC格式输入文件（温度、露点、风速、辐射）读取及变量提取均正常，但执行计算约1分钟后触发如下错误：

Error in if (is.leapyear(year)) dpy = 366 else dpy = 365: Missing value where TRUE/FALSE is required

已确认时间数据符合手册要求，尝试过多个自定义闰年函数（如下），但错误仍未解决。需明确解决方向、是否需额外闰年函数，以及如何规避布尔值缺失错误。

原始代码 #

library(HeatStress)
library(ncdf4)
library(lubridate)

ncfile_t2m <- nc_open("t2m.nc") 
ncfile_d2m <- nc_open("d2m.nc")  
ncfile_wind <- nc_open("wind.nc")
ncfile_ssrd <- nc_open("ssrd.nc")  

tas <- ncvar_get(ncfile_t2m, "t2m")  
dewp <- ncvar_get(ncfile_d2m, "d2m")  
wind <- ncvar_get(ncfile_wind, "wind_speed")  
radiation <- ncvar_get(ncfile_ssrd, "ssrd") 

lon <- ncvar_get(ncfile_t2m, "longitude") 
lat <- ncvar_get(ncfile_t2m, "latitude") 
time <- ncvar_get(ncfile_t2m, "time")
origin <- ymd_hms("1900-01-01 00:00:00")  # 时间原点
time_posix <- origin + hours(time)  # 将小时数转换为POSIX时间


nc_close(ncfile_t2m)
nc_close(ncfile_d2m)
nc_close(ncfile_wind)
nc_close(ncfile_ssrd)


wbgt_result <- wbgt.Liljegren(tas = tas, dewp = dewp, wind = wind, radiation = radiation, 
dates = time_posix, lon = lon[1], lat = lat[1])

print(wbgt_result)

尝试过的自定义闰年函数 #

函数1 #

days_per_year <- function(year) {
  if ((year %% 4 == 0 & year %% 100 != 0) | year %% 400 == 0) {
    return(366)
  else {
    return(365)
  }
}

函数2 #

is.leapyear <- function(year) {
  leap_year <- function(y) {
    as.logical((y %% 4 == 0 & year %% 100 != 0) | year %% 400 == 0)
  }
  if (anyNA(year)) {
    return(rep(FALSE, length(year)))  # 显式处理NA值
  }
  sapply(year, leap_year)
}

函数3 #

is.leapyear <- function(year) {
  if (is.na(year)) {
    return(FALSE)  # 显式处理NA值
  }
  as.logical((year %% 4 == 0 & year %% 100 != 0) | year %% 400 == 0)
}

解决方案 #

问题根源 #

wbgt.Liljegren内部调用的默认is.leapyear函数未处理NA值，若time_posix中存在NA年份，会导致条件判断时因缺失布尔值触发错误。此外，自定义函数未正确覆盖包内原函数，或未处理向量化输入的NA场景。

具体解决步骤 #

1. 检查并清理时间数据
   先确认time_posix是否存在NA值：

   sum(is.na(time_posix))
   

   若存在NA，回溯NC文件的time变量，确认时间序列完整性，必要时过滤或补全缺失时间点。

2. 覆盖包内的is.leapyear函数
   在加载HeatStress包后，立即定义一个支持向量化且处理NA的函数，覆盖原函数：

   library(HeatStress)
   # 覆盖内置的is.leapyear函数
   is.leapyear <- function(year) {
     ifelse(is.na(year), FALSE, (year %% 4 == 0 & year %% 100 != 0) | year %% 400 == 0)
   }
   

   该函数能批量处理年份向量，自动将NA转换为FALSE，避免条件判断出错。

3. 验证输入维度匹配
   确认tas、dewp、wind、radiation的时间维度与time_posix长度一致：

   length(time_posix) == dim(tas)[3]  # 假设时间是第三维度
   

   维度不匹配可能导致隐性数据缺失，进而触发NA相关错误。

4. 验证年份提取结果
   提取时间中的年份并检查NA：

   year_vec <- year(time_posix)
   sum(is.na(year_vec))
   

   若仍有NA，需检查origin与NC文件的时间基准是否一致（部分NC文件时间基准可能为1970-01-01）。

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