使用 gstat 和 automap 包的克里金法 - 复制教程时的问题

Posted 2023-03-12

【中文标题】使用 gstat 和 automap 包的克里金法 - 复制教程时的问题

【英文标题】：Kriging using the gstat and automap packages - problems when copying a tutorial

【发布时间】：2022-01-03 04:06:48

【问题描述】：

对于我研究站点中的几个位置，我有水深 (m) 数据，并且我正在尝试使用克里金法将深度插入到我没有数据的位置。我找到了 Wilke 博士 here 写的一篇有用的博客文章，并尝试将他的代码应用于我的数据集。但是，当我运行所有代码并绘制结果时，所有图都返回为空（灰色框；深度样本以灰点返回，而不是在连续尺度上反映深度的颜色）。谁能告诉我我做错了什么？下面是我尝试过的代码。我怀疑问题出在第 2 步，因为这里的代码输出开始偏离博客代码输出。

设置

depth <- structure(list(X = c(668875.407301466, 675555.822548817, 677165.984257575, 
671427.139500822, 670565.010067336, 676196.888461476, 669435.428473297, 
677288.114994538, 679355.228144992, 677003.998083734, 677585.341252976, 
679549.867927876, 677659.023812075, 679730.677304853, 677749.442051316, 
669617.624346464, 671569.53262659, 673306.072271015, 670280.828937677, 
677894.966356866, 677192.79791361, 672207.85937462, 674600.067203517, 
674244.297764696, 673763.964359285, 669723.003482714, 677106.449499354, 
671172.528780568, 673040.137152061, 672415.374390262, 675286.888306849, 
674272.583595863, 672978.2232579, 672207.731105759, 673057.984061267, 
673098.008775404, 670469.846198017, 670461.213900251, 671449.443237703, 
671476.810705895, 672429.201320869, 672426.444205806, 673397.966944431, 
673397.001177047, 676506.198002199, 679013.454950302, 668360.513631175, 
669309.125013407, 675540.056558172, 673946.851010242, 668736.385947233, 
669581.113779484, 669531.260804133, 668572.915079444, 668562.271762613
), Y = c(2843306.80511914, 2846996.53521373, 2847483.31959936, 
2843061.483504, 2842817.70752641, 2853858.66887037, 2843862.35002759, 
2844357.28754673, 2844326.38872235, 2843760.75211182, 2845096.92928085, 
2844999.32573033, 2845739.40590226, 2845654.35288888, 2846405.37514488, 
2846251.04899388, 2856292.70424562, 2855857.84577301, 2848289.20711897, 
2848868.04137537, 2851322.64880942, 2853330.72832414, 2854696.08607969, 
2850951.06030575, 2849409.17353187, 2840685.32365573, 2837738.38794808, 
2850715.9222808, 2846526.89253892, 2842945.78513155, 2844114.61778685, 
2844613.08986426, 2844813.27348601, 2845405.92798845, 2846082.85831243, 
2848010.02032868, 2841631.05374873, 2840583.99941216, 2841647.00399675, 
2840617.02709824, 2841650.84197688, 2840627.11766348, 2841420.84288763, 
2840419.29474805, 2843303.19107449, 2843694.67766026, 2838715.91143071, 
2838736.83548694, 2839850.31020878, 2842589.10652076, 2840663.87712246, 
2842046.19194071, 2839662.52414273, 2842072.11298559, 2839670.12536644
), Z = c(27, 2, 1, 1.5, 1.5, 4, 6, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 
1.5, 1.5, 14, 26, 4, 12, 3, 3, 7, 4, 1, 1, 14.4, 2, 6, 1, 1, 
1, 1, 1, 1, 1, 1, 5, 4, 3, 3, 2, 4, 2, 3, 2, 2, 27.8, 10, 2, 
1.5, 28, 8, 16.6, 28, 28.2)), row.names = c(NA, -55L), class = c("tbl_df", 
"tbl", "data.frame"))

X 和 Y 是包含空间坐标的列，Z 包含这些位置的深度值。

# Convert to sf because that is the best way to store spatial points

  depth_sf <- st_as_sf(depth, coords = c("X", "Y"), crs = 32617) %>% 
      cbind(st_coordinates(.))

创建变异函数

# We will discuss later, what Z~1 does actually mean in this context
    
  v_emp_OK <- gstat::variogram(
     Z~1,
      as(depth_sf, "Spatial") # switch from sf to sp
  )

# automap's autofitVariogram actually produces more info than we need. I will only keep the var_model part.
    
  v_mod_OK <- automap::autofitVariogram(Z~1, as(depth_sf, "Spatial"))$var_model

# To inspect the automatic fit that was chosen for us we can use automap's excellent build in methods for base::plot

  plot(automap::autofitVariogram(Z~1, as(depth_sf, "Spatial")))

定义目标网格

# technically we already have a grid from the initial dataset, but as we are still working under the pretense that our only available data are the simulated observation wells, we will construct our grid from that object.'

# Step 1: define a grid based on the bounding box of our observations
  grd_100_sf <- depth_sf %>% 
    st_bbox() %>% 
    st_as_sfc() %>% 
    st_make_grid(
     cellsize = c(100, 100), # 100m pixel size
      what = "centers"
    ) %>%
    st_as_sf() %>%
    cbind(., st_coordinates(.))

#Step 2: making our grid work for gstat

 grd_100_sp <- as(grd_100_sf, "Spatial") # converting to sp format
 gridded(grd_100_sp) <- TRUE             # informing the object that it is a grid
 grd_100_sp <- as(grd_100_sp, "SpatialPixels") # specifying what kind of grid

克里金法

# Ordinary Kriging
OK <- krige(
  Z~1,                       # Z is our variable and "~1" means "depends on mean"
  as(depth_sf, "Spatial"), # input data in sp format
  grd_100_sp,                # locations to interpolate at
  model = v_mod_OK           # the variogram model fitted above
)

## [using ordinary kriging]

# Simple Kriging
SK <- krige(
  Z~1,                       # Z still depends on mean
  beta = mean(depth$Z), # but this time we know the mean's value
  as(depth_sf, "Spatial"), # input data in sp format
  grd_100_sp,                # locations to interpolate at
  model = v_mod_OK           # the variogram model fitted above
)

## [using simple kriging]

# Universal Kriging
# Implementing this method is somewhat different.
# we no longer assume that Z is essentially depending on a single mean but
# rather on the position of the interpolation location within our target grid
UK <- krige(
  Z~coords.x1+coords.x2, # Think "Z~X+Y" but sp conversion alters variable naming
  as(depth_sf, "Spatial"), # input data in sp format (`X` --> `coords.x1`)
  grd_100_sp,                # locations to interpolate at
  model = autofitVariogram(  # we need an appropriate variogram fit
    Z~X+Y,                   # here we can keep "X+Y" - it's just how it is
    as(depth_sf, "Spatial")
  )$var_model
)

## [using universal kriging]

# I'll also add an inverse distance weighted model to provide a baseline
# for model evaluation
# Note how the only difference to Ordinary Kriging is the absence of a
# fitted variogram model
idwres <- idw(
  Z~1,                       # idw also depends on mean
  as(depth_sf, "Spatial"), # input data in sp format
  grd_100_sp,                # locations to interpolate at
) 

检查结果

# A function to plot rasters
plot_my_gstat_output <- function(raster_object, object_name)
  
  df <- rasterToPoints(raster_object) %>% as_tibble()
  colnames(df) <- c("X", "Y", "Z")
  
  ggplot(df, aes(x = X, y = Y, fill = Z)) +
    geom_raster() +
    ggtitle(label = object_name) +
    scale_fill_viridis(option = "B", limits = c(50, 100)) +
    theme_void() +
    theme(
      plot.title = element_text(hjust = 0.5)
    )


p_idw <- plot_my_gstat_output(raster(idwres), "IDW")
p_SK <- plot_my_gstat_output(raster(SK), "Simple Kriging")
p_OK <- plot_my_gstat_output(raster(OK), "Ordinary Kriging")
p_UK <- plot_my_gstat_output(raster(UK), "Universal Kriging")

# I also want to display sampling locations
p_depth <- ggplot(
  data = depth,
  mapping = aes(x = X, y = Y, color = Z)
) +
  geom_point(size = 3) + 
  scale_color_viridis(option = "B",  limits = c(50, 100)) +
  ggtitle(label = "Depths Sampled") +
  theme_void() +
  theme(
    plot.title = element_text(hjust = 0.5)
  )

# This works because of library(patchwork)
(p_depth + p_idw) / 
  (p_SK + p_OK + p_UK) + 
  plot_layout(guides = 'collect')

【问题讨论】：

除了在 ggplot 中调用的深度范围值之外，您的代码和博客代码之间的一个区别是您的对象 v_mod_OK 缺少变量 kappa。不确定代码中其他地方的调用位置，但有一点需要注意。

@JessicaBurnett 感谢您的指出。我不确定没有 kappa 值的后果是什么，但我认为在这种情况下，它是一个与 Matern 模型有关的参数？因此不确定这是否真的那么重要，因为变异函数确定了球形模型以提供最佳拟合？

【参考方案1】：

看起来罪魁祸首在代码的ggplot 部分。使用scale_color_viridis() 时，您需要更改limits 以适合您的数据。由于数据中的深度为 0 到 28，因此您应该将限制从 c(50, 100) 更改为 c(0, 30)。

# A function to plot rasters
plot_my_gstat_output <- function(raster_object, object_name)
  
  df <- rasterToPoints(raster_object) %>% as_tibble()
  colnames(df) <- c("X", "Y", "Z")
  
  ggplot(df, aes(x = X, y = Y, fill = Z)) +
    geom_raster() +
    ggtitle(label = object_name) +
    scale_fill_viridis(option = "B", limits = c(0, 30)) +
    theme_void() +
    theme(
      plot.title = element_text(hjust = 0.5)
    )


p_idw <- plot_my_gstat_output(raster(idwres), "IDW")
p_SK <- plot_my_gstat_output(raster(SK), "Simple Kriging")
p_OK <- plot_my_gstat_output(raster(OK), "Ordinary Kriging")
p_UK <- plot_my_gstat_output(raster(UK), "Universal Kriging")

# I also want to display sampling locations
p_depth <- ggplot(
  data = depth,
  mapping = aes(x = X, y = Y, color = Z)
) +
  geom_point(size = 3) + 
  scale_color_viridis(option = "B",  limits = c(0, 30)) +
  ggtitle(label = "Depths Sampled") +
  theme_void() +
  theme(
    plot.title = element_text(hjust = 0.5)
  )

# This works because of library(patchwork)
(p_depth + p_idw) / 
  (p_SK + p_OK + p_UK) + 
  plot_layout(guides = 'collect')

输出

【讨论】：

感谢您的发现！在逐步进行之前，我没有仔细检查所有代码。如果我这样做，我可能会发现问题。

@FlyingDutch 没问题！我实际上已经参考并使用过这篇确切的博文，所以已经很熟悉了。