数据清理后自动地图插值的问题

Question

我想对原始数据进行空间数据插值 使用 automap 库（使用 paar 库中的 depurate 函数清理数据后）网格（由 x 和 y 变量定义）。

安装后，我收到以下警告消息：

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

非常感谢有关此问题的任何帮助。

可重现的示例

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

# paar library instalation
# 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))

Answer 1

如果您在输入数据的确切位置进行克里金预测，则理论上该点的方差同样为零。因为数据告诉您，您知道该点表面的确切值。逐点测量误差为零。

但在计算上，有时这些零最终会变成非常非常小的负数，并且当求平方根以获得标准差时，在取负数的平方根时会出现错误。这就是您所得到的。

> subset(kriging_result$krige_output, is.na(var1.stdev))
             coordinates var1.pred      var1.var var1.stdev
11   (312432.8, 5800264)  3.671085 -1.387779e-17        NaN
16   (312412.8, 5800274)  3.694282 -2.775558e-17        NaN
17   (312422.8, 5800274)  3.669411 -4.163336e-17        NaN
28   (312442.8, 5800284)  3.627693 -4.163336e-17        NaN
36   (312412.8, 5800294)  3.679696 -4.163336e-17        NaN

因此在数据位置预测克里金法有点毫无意义。

如果您只想使用克里金法来填充清理过程中产生的孔，请将所有剩余的数据点作为输入，并将孔位置作为预测位置，然后将两者合并。数据位置的不确定性应为零，孔处的不确定性应为正方差和正标准差。

Answer 2

这是解决方案的当前版本。

任何如何改进或纠正代码的想法都非常受欢迎。

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

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)
}

# cleaning the data ---
tg_cleaned <- clean_variable(dt_sf, "Tg")

# identifying holes (missing locations after cleaning) ---
holes <- dt_sf[!dt_sf$Tg %in% tg_cleaned$Tg, ]

# kriging function - prediction only at "hole" locations
interpolate_holes <- function(cleaned_data, holes) {
  cleaned_data_sp <- as(cleaned_data, "Spatial")
  holes_sp <- as(holes, "Spatial")
  
  vgm_model <- autofitVariogram(Tg ~ 1, cleaned_data_sp)$var_model
  kriging_result <- krige(Tg ~ 1, cleaned_data_sp, holes_sp, model = vgm_model)
  
  holes_sp@data$Tg_pred <- kriging_result@data$var1.pred
  holes_sp@data$Tg_var <- kriging_result@data$var1.var
  
  return(st_as_sf(holes_sp))
}

# kriging to the holes ---
holes_with_predictions <- interpolate_holes(tg_cleaned, holes)

tg_cleaned_sf <- st_as_sf(tg_cleaned)
tg_cleaned_sf$Tg_pred <- tg_cleaned_sf$Tg
tg_cleaned_sf$Tg_var <- NA
holes_with_predictions$Tg <- NA
holes_with_predictions <- holes_with_predictions[, colnames(tg_cleaned_sf)]
final_data_sf <- rbind(tg_cleaned_sf, holes_with_predictions)

# spatial map ---
ggplot(final_data_sf) +
  geom_sf(aes(color = Tg_pred)) +
  scale_color_viridis_c(option = "turbo") +
  theme(plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(vjust = 1, hjust = 1, angle = 25))