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

Question

对于我研究站点中的几个位置，我有水深 (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')

Answer 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')

输出