1

I am trying to produce kriging models from elevations points using the gstat package. I can fit models to the empirical variogram using exponential, spherical, and Gaussian curves. However, despite the Gaussian curve arguably fitting the variogram best, it produces terrible outputs compared to the other two models.

CSV of example data can be found here: https://drive.google.com/file/d/1Lwmvlwv0h1kyXbN60kOobn4a2_PkpBh0/view?usp=sharing

Below is the code I am using and the outputs.

library("terra")
library("stars")
library("tmap")
library("sf")
library("spatstat")
library("sp")
library("dplyr")
library("gstat")

data <- read.csv("./TestPoints.csv")

pnts <- st_as_sf(data, coords = c("X","Y"), crs = 26910)

# Create an empty grid
grd <- as.data.frame(spsample(as_Spatial(pnts), "regular", n=100000))
names(grd)       <- c("X", "Y")
coordinates(grd) <- c("X", "Y")
gridded(grd)     <- TRUE  
fullgrid(grd)    <- TRUE  
crs(grd) <- crs(pnts)

#Spatial Interpolation with Kriging
f.0 <- as.formula(Z ~ 1) 

#Create variogram and models
var.smpl <- variogram(f.0, pnts, cloud = FALSE, cutoff = 300) 

exp.fit  <- fit.variogram(var.smpl, fit.ranges = FALSE, fit.sills = FALSE,
                          vgm(model="Exp"))

sph.fit  <- fit.variogram(var.smpl, fit.ranges = FALSE, fit.sills = FALSE,
                          vgm(model="Sph"))

gau.fit  <- fit.variogram(var.smpl, fit.ranges = FALSE, fit.sills = FALSE,
                          vgm(model="Gau"))

#Plot variogram and models
plot(gamma ~ dist, var.smpl, ylim = c(0, 1.05*max(var.smpl$gamma)), col='blue', ylab = 'semivariance', xlab  = 'distance')
lines(variogramLine(exp.fit, 300), lty =1, lwd=1)
lines(variogramLine(sph.fit, 300), lty=2, lwd =1)
lines(variogramLine(gau.fit, 300), lwd=2, lty=2)
legend(5, 140, c("Exponential model", "Spherical model", "Gaussian model"), lty = c(1,2,2), lwd = c(1,1,2))


# Perform the EXP interpolation 
exp.krg <- krige(f.0, as_Spatial(pnts), grd, exp.fit)

# Plot the map
tm_shape(exp.krg) + 
  tm_raster(col = "var1.pred", n=10,palette = "-RdBu",
            title="TITLE") + 
  tm_shape(pnts) + tm_dots(size=0.01, alpha = 0.9) +
  tm_legend(legend.outside=TRUE)


# Perform the SPH interpolation 
sph.krg <- krige(f.0, as_Spatial(pnts), grd, sph.fit)

# Plot the map
tm_shape(sph.krg) + 
  tm_raster(col = "var1.pred", n=10,palette = "-RdBu",
            title="TITLE") + 
  tm_shape(pnts) + tm_dots(size=0.01, alpha = 0.9) +
  tm_legend(legend.outside=TRUE)


# Perform the GAU interpolation 
gau.krg <- krige(f.0, as_Spatial(pnts), grd, gau.fit)

# Plot the map
tm_shape(gau.krg) + 
  tm_raster(col = "var1.pred", n=10,palette = "-RdBu",
            title="TITLE") + 
  tm_shape(pnts) + tm_dots(size=0.01, alpha = 0.9) +
  tm_legend(legend.outside=TRUE)

ModelFit

EXP Model Map

SPH Model Map

GAU Model Map

1 Answer 1

0

I think what's happening is that the big blue blob (and the little red bits) is forcing the colours to rescale so you can't actually see the detail. Those regions are probably away from your data (with respect to your correlation scale) and I suspect have even larger standard errors, so that the big blue blob is not significantly different from the mean of your data. If you consider the uncertainty you'll probably find that this isn't as terrible as the other covariance fits.

If you can do the plot with the same colour scale as the other maps then I think you'll see the same sort of detail as the other maps in the locations where your data are.

One thing you could do is crop the output only to regions that are within a certain distance of your data, or within a concave hull constructed on the data points.

Kriging can go a bit wild once it gets away from your data and before its had a chance to settle on a value of the global mean at very large distances, which is what it should do. I suspect there's a little bit of downward curvature at the edge of the data and the specifics of the Gaussian correlation structure force it into that huge blue blob dip before it can come back up to hit the mean.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.