rkriginggstatautomap

Kriging using the gstat and automap packages - problems when copying a tutorial

For several locations in my study site i have water depth (m) data, and i'm trying to use kriging to interpolate depth to locations for which i do not have data. I found a useful blog post written by Dr. Wilke here, and tried to apply his code to my data set. However, when i run all the code and plot the results, the plots all come back empty (grey boxes; the depths samples comes back with grey dots instead of colours that reflect depth on continuous scale). Could anyone please tell me what i'm doing wrong? Below is the code that i tried. I suspect the issue lies with step 2, as the code output starts to deviate here from the blog code output.

Setup
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 and Y are columns that contain the spatial coordinates, Z contains the values for depth at those locations.

# 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(.))
Creating a variogram
# 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")))
Defining a target grid
# 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
Kriging
# 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
)
Inspect the results
# 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')
Solution

  • It looks like the culprit is in the ggplot sections of the code. When using scale_color_viridis(), you need to change the limits to fit your data. Since the depths in your data are 0 to 28, then you should change the limits from c(50, 100) to something like 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')

    Output

    enter image description here