Kriging with gstat package in R produces empty interpolation points

gstatkrigingrsf

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

enter image description here

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

enter image description here

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)

enter image description here

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

Best Answer

I don't really recommend doing kriging in lat-long degree coordinates. Project your data to a cartesian coordinate system and try again. For example, if I convert to 3857 (web mercator) then it works.

> int_pnts=st_transform(int_pnts, 3857)
> input_data_sf = st_transform(input_data_sf, 3857)
> lzn.vgm <- variogram(log(OM)~1, input_data_sf)
> lzn.fit <- fit.variogram(lzn.vgm, vgm('Mat'), fit.kappa = TRUE)

Some warnings we'll ignore, and then kriging produces meaningful results:

> lzn.kriged <- krige(log(OM) ~ 1, input_data_sf, int_pnts, model = lzn.fit)
[using ordinary kriging]
> head(lzn.kriged)
Simple feature collection with 6 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -10757430 ymin: 4738785 xmax: -10757310 ymax: 4738918
Projected CRS: WGS 84 / Pseudo-Mercator
  var1.pred   var1.var                  geometry
1 0.6545492 0.04198732 POINT (-10757334 4738803)
2 0.3821123 0.04326622 POINT (-10757376 4738785)
3 0.8089032 0.04263483 POINT (-10757315 4738882)
4 0.6558402 0.04288256 POINT (-10757384 4738896)

CRS 3857 is not the best choice, its just something that is likely to "work" in terms of "producing a result" in most places on the globe. You should maybe find a local coordinate system for your data, or use a UTM zone appropriate for your longitude.

Related Question