数据清洗后使用automap空间插值出现NaN警告问题求助

问题背景 #

我希望使用automap库，基于原始网格（由x、y变量定义）对经paar库depurate函数清洗后的空间数据进行插值，拟合后出现如下警告：

Warning message:
In sqrt(krige_result$var1.var) : NaNs produced

可复现代码 #

library(sf)
library(paar)
library(automap)
library(ggplot2)
library(viridis)

# paar 库安装
# install.packages("devtools")
# devtools::install_github("PPaccioretti/paar")

data("wheat", package = "paar")
dt_sf <- st_as_sf(wheat, coords = c("x", "y"), crs = 32720)

clean_variable <- function(data_sf, var_name) {
  cleaned_data <- depurate(data_sf,
                           y = var_name,
                           toremove = c("edges", "outlier", "inlier"),
                           buffer = -10)
  return(cleaned_data$depurated_data)
}

tg_cleaned <- clean_variable(dt_sf, "Tg")
original_grid <- dt_sf
original_grid_sp <- as(original_grid, "Spatial")

interpolate_tg <- function(cleaned_data, original_grid_sp) {
  cleaned_data_sp <- as(cleaned_data, "Spatial")
  kriging_result <- autoKrige(Tg ~ 1, cleaned_data_sp, new_data = original_grid_sp)
  original_grid_sp@data[["tg_interpol"]] <- kriging_result$krige_output@data$var1.pred
  original_grid_sp@data[["tg_var"]] <- kriging_result$krige_output@data$var1.var
  return(original_grid_sp)
}

final_data <- interpolate_tg(tg_cleaned, original_grid_sp)
final_data_sf <- st_as_sf(final_data)

ggplot(final_data_sf) +
  geom_sf(aes(color = tg_interpol)) +
  scale_color_viridis_c(option = "turbo") +
  ggtitle("spatial distribution of var1.pred") +
  theme_void() +
  theme(plot.title = element_text(hjust = 0.5))

ggplot(final_data_sf) +
  geom_sf(aes(color = tg_var)) +
  scale_color_viridis_c(option = "turbo") +
  ggtitle("spatial distribution of var1.var") +
  theme_void() +
  theme(plot.title = element_text(hjust = 0.5))

问题原因 #

这个警告的核心是krige_result$var1.var（克里金预测方差）中出现了负值，对负数开平方会直接产生NaN。常见诱因包括：

• 数据清洗后有效样本点过少，导致变异函数拟合不稳定，预测方差计算出现异常
• 部分预测点位于清洗后数据的空间范围之外，克里金外插时容易出现方差负值
• 自动拟合的变异函数参数不合理（如块金值、基台值、范围设置偏差）

解决步骤 #

1. 调整数据清洗策略，保留足够样本 #

先检查清洗后的数据量：

nrow(tg_cleaned)

如果样本量过少，修改depurate函数的参数，减少不必要的数据剔除：

clean_variable <- function(data_sf, var_name) {
  cleaned_data <- depurate(data_sf,
                           y = var_name,
                           toremove = c("outlier"),  # 仅保留异常值剔除
                           buffer = -5)  # 缩小边缘点剔除的范围
  return(cleaned_data$depurated_data)
}

2. 限制预测范围，避免外插 #

只在清洗后数据的空间边界内进行插值，过滤原始网格中的外插点：

# 获取清洗后数据的空间边界
clean_bbox <- st_bbox(tg_cleaned)
# 过滤原始网格，仅保留边界内的点
filtered_grid <- original_grid[st_within(original_grid, st_as_sfc(clean_bbox), sparse = FALSE), ]
original_grid_sp <- as(filtered_grid, "Spatial")

3. 手动指定变异函数模型 #

如果自动拟合的变异函数不稳定，手动定义更合理的模型参数：

# 先查看自动拟合的变异函数
variogram_model <- autofitVariogram(Tg ~ 1, as(tg_cleaned, "Spatial"))
print(variogram_model$var_model)

# 手动指定模型（示例为球状模型，参数需根据实际数据调整）
manual_vgm <- vgm(psill = 8, model = "Sph", range = 150, nugget = 1.2)

# 使用手动模型进行插值
kriging_result <- autoKrige(Tg ~ 1, as(tg_cleaned, "Spatial"), new_data = original_grid_sp, model = manual_vgm)

4. 处理预测方差的负值 #

如果以上方法仍有负值，直接替换为0或极小值，避免NaN：

# 在插值函数中添加负值处理逻辑
interpolate_tg <- function(cleaned_data, original_grid_sp) {
  cleaned_data_sp <- as(cleaned_data, "Spatial")
  kriging_result <- autoKrige(Tg ~ 1, cleaned_data_sp, new_data = original_grid_sp)
  
  # 替换预测方差中的负值为0
  pred_var <- kriging_result$krige_output@data$var1.var
  pred_var[pred_var < 0] <- 0
  
  original_grid_sp@data[["tg_interpol"]] <- kriging_result$krige_output@data$var1.pred
  original_grid_sp@data[["tg_var"]] <- pred_var
  return(original_grid_sp)
}

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