Kriging with gstat in R: Resulting object is blank with no interpolated values

Question

I have kriged a dataset with the gstat package in R, but it has produced an empty variable as a result. I have 1,167 measurement points within a field, and I am trying to interpolate them across 3,464 interpolation points within the same field. How can I achieve a Kriged result with estimated values at each point?

Note: The input_data and interpolation points can be found within a text file at this link; the contents within the file just need to be copied and pasted into an R window and run to generate the data frames. In addition, if desired, the field boundary shape file referred to later in this post can be found here.

These are the points that will be interpolated:

#Import libraries
library(pacman)
p_load(raster, 
       sf, 
       dplyr, 
       ggplot2, 
       scales, 
       magrittr, 
       gstat, 
       gridExtra, 
       raster,
       sp,
       automap,
       mapview,
       leaflet,
       rgdal)

#Graph points for interpolation
input_data %>% 
  as.data.frame %>% 
  ggplot(aes(latitude, longitude)) +
  geom_point(aes(size = OM), color = 'red', alpha = 3/4) +
  ggtitle('Organic Matter Concentration') +
  coord_equal() +
  theme_bw()

Convert data frame into a spatial object

input_data_sf <- st_as_sf(input_data, coords = c('longitude', 'latitude'), crs = 4326)
crs(input_data_sf)

glimpse(input_data_sf)

Coordinate Reference System:
Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
WKT2 2019 representation:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]] 

The projection we are working with will be WGS 1984.

Making sure that the data was converted properly:

plot(st_geometry(input_data_sf))

That looks good.

Fitting the variogram using gstat and a Matern model:

lzn.vgm <- variogram(log(OM)~1, input_data_sf)

lzn.fit <- fit.variogram(lzn.vgm, vgm('Mat'), fit.kappa = TRUE)

The following warning message is produced after runnign the fit.variogram() command:

Warning messages:
1: In fit.variogram(o, m, fit.kappa = FALSE, fit.method = fit.method,  :
  No convergence after 200 iterations: try different initial values?
2: In fit.variogram(o, m, fit.kappa = FALSE, fit.method = fit.method,  :
  No convergence after 200 iterations: try different initial values?
3: In fit.variogram(o, m, fit.kappa = FALSE, fit.method = fit.method,  :
  No convergence after 200 iterations: try different initial values?

Plotting the variogram:

plot(lzn.vgm, lzn.fit)

After generating the semivariogram, we now must convert the input points (the points for interpolation) to a spatial object:

int_pnts <- st_as_sf(int_pnts, coords = c('POINT_X', 'POINT_Y'), crs = 4326)
crs(int_pnts)

glimpse(int_pnts)

Visualizing the measured points and interpolated points:

#Read field boundary shapefile
shp <- st_read("field_boundary.shp")

plot1 <- ggplot() +
  geom_sf(data = shp) +
  geom_sf(data = velv) +
  ggtitle("Sample points in field")

#Run this code if you don't download the field boundary
#plot1 <- ggplot() +
  #geom_sf(data = velv) +
  #ggtitle("Sample points in field")

plot1

plot2 <- ggplot() +
  geom_sf(data = shp) +
  geom_sf(data = int+pnts) +
  ggtitle("Interpolation points")

#Run this code if you don't download the field boundary:
#plot2 <- ggplot() +
  #geom_sf(data = int+pnts) +
  #ggtitle("Interpolation points")

plot2

Finally, perform the kriging:

lzn.kriged <- krige(log(OM) ~ 1, input_data_sf, int_pnts, model = lzn.fit)

The problem is here: When converting the lzn.kriged output to a dataframe to investigate how it turned out, both of the first two columns (var1 predict and actual) have NA values beside them. How can I fix this so I can graph my interpolation?

#Convert kriged output to data frame
krig <- lzn.kriged %>% 
  as.data.frame

Answer 1

When running the kriging I get some errors on Covariance matrix singular at location ... (see Kriging with gstat : "Covariance matrix singular at location" with predict). That means that your data input has duplicated locations. I followed this tidyverse approach to dedupe the dataset and it just worked fine ;):

source("input_data.R")

library(pacman)

p_load("sf", "dplyr", "gstat", "ggplot2")
nrow(input_data)
# [1] 1167

# Input data_sf dedupe
input_data_dedupe <- input_data %>%
  group_by(longitude, latitude) %>%
  # Just for illustatrion, you may sum or use another approach 
  summarise(OM = mean(OM))


nrow(input_data_dedupe)

# [1] 1165

# Ok, there were two duplicated coordinates!


# Now the rest of the analysis with deduped data
input_data_sf <- st_as_sf(input_data_dedupe,
  coords = c("longitude", "latitude"), crs = 4326
)
int_pnts_sf <- st_as_sf(int_pnts, coords = c("POINT_X", "POINT_Y"), crs = 4326)

# Variogram
lzn.vgm <- variogram(log(OM)~1, input_data_sf)


lzn.fit <- fit.variogram(lzn.vgm, vgm("Mat"), fit.kappa = TRUE)


# Krige
lzn.kriged <- krige(log(OM) ~ 1, input_data_sf, int_pnts_sf, model = lzn.fit)

lzn.kriged %>%
  as_tibble()


# # A tibble: 3,464 x 3
# var1.pred var1.var             geometry
# <dbl>    <dbl>          <POINT [°]>
#   1     0.646   0.0418  (-96.63478 39.1195)
# 2     0.417   0.0432 (-96.63515 39.11937)
# 3     0.812   0.0425  (-96.6346 39.12005)
# 4     0.654   0.0428 (-96.63523 39.12015)
# 5     0.654   0.0431 (-96.63568 39.12002)
# 6     0.686   0.0449  (-96.63485 39.1203)
# 7     0.905   0.0424  (-96.63453 39.1199)
# 8     0.633   0.0427  (-96.6352 39.12017)
# 9     0.657   0.0413 (-96.63458 39.11935)
# 10     0.853   0.0427 (-96.63548 39.11982)
# # ... with 3,454 more rows

ggplot(lzn.kriged) +
  geom_sf(aes(color=var1.pred))