R, Gstat predict() error

Question

I'm trying to use gstat to predict spatially autocorrelated random variables. I'm using these sea surface temperature data. I first fit a variogram model and use this in predict. The approach works fine if I use the global dataset, but fails with some regional subsets.

My code is below:

env.dat.basin <- env.dat
# remove duplicate locations
env.dat.basin <- env.dat.basin[!duplicated(env.dat.basin[,-c(3,4)]),]

# convert to a spatial points df
coordinates(env.dat.basin) <- ~x + y
proj4string(env.dat.basin) <- CRS("+proj=longlat +datum=WGS84")

# fit variogram
v <- variogram(SST~1, env.dat.basin)
fit.v <- fit.variogram(v, model = vgm(psill = 30, model = 'Gau', range = 5000, nugget = 0), fit.sills = T, fit.ranges = T)
plot(v, fit.v)

# define a gstat object (spatial model) based on variogram fit
g.dummy <- gstat(formula=SST~1, data = env.dat.basin, beta = 1, dummy=T, model=fit.v, nmax = 30)

# make nsim simulations based on the stat object
yy <- predict(g.dummy, newdata=env.dat.basin, nsim=1)

this works, and if I subset the global dataset

env.dat.basin <- env.dat[which(env.dat$ocean == 7 | env.dat$ocean == 513), ]

it works fine too, but with a different subset

env.dat.basin <- env.dat[c(which(env.dat$ocean == 11), which(env.dat$ocean == 7 & env.dat$y <10)), ]

which looks like this:

the prediction returns many warnings that cause R to crash if nsim >1

Warning messages:
1: In predict.gstat(g.dummy, newdata = env.dat.basin, nsim = 1) :
  singular simulation covariance matrix
2: In predict.gstat(g.dummy, newdata = env.dat.basin, nsim = 1) :
  Covariance matrix singular at location [6.3,-21.4333,0]: skipping...
3: In predict.gstat(g.dummy, newdata = env.dat.basin, nsim = 1) :
  singular simulation covariance matrix
4: In predict.gstat(g.dummy, newdata = env.dat.basin, nsim = 1) :
  Covariance matrix singular at location [-10.6,-34.45,0]: skipping...
5: In predict.gstat(g.dummy, newdata = env.dat.basin, nsim = 1) :
  Covariance matrix singular at location [-14.885,-15.5333,0]: skipping...

I read elsewhere that this could be due to the presence of duplicate observations, but I think I removed those and the code works for the global dataset.

Answer 1

It is a consequence of using a Gaussian variogram model with zero nugget; this easily leads to near-singular covariance matrices even if points don't (exactly) overlap. Solutions: choose a different variogram model (Matern, with strong smoothing?) or add a small nugget effect (very small may already help).