Error using gaussian model with kriging

[jwkelley]

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)

And outputs:

[edzer]

That is a well-known property of the Gaussian variogram - adding a small nugget often helps.