Search code examples
rkriginggstat

Kriging with gstat : "Covariance matrix singular at location" with predict

I am trying to do an estimation by kriging with gstat, but can never achieve it because of an issue with the covariance matrix. I never have estimates on the locations I want, because they are all skipped. I have the following warning message, for each location :

1: In predict.gstat(g, newdata = newdata, block = block, nsim = nsim,  : 
Covariance matrix singular at location [-8.07794,48.0158,0]: skipping...

And all estimates are NA.

So far I've browsed many related StackOverflow threads, but none resolved my problems (https://gis.stackexchange.com/questions/222192/r-gstat-krige-covariance-matrix-singular-at-location-5-88-47-4-0-skipping ; https://gis.stackexchange.com/questions/200722/gstat-krige-error-covariance-matrix-singular-at-location-917300-3-6109e06-0 ; https://gis.stackexchange.com/questions/262993/r-gstat-predict-error?rq=1)

I checked that :

  • there is actually a spatial structure in my dataset (see bubble plot with code below)

  • there are no duplicate locations

  • the variogram model is not singular and has a good fit to the experimental variogram (see plot with code below)

  • I also tried several values of range, sill, nugget and all the models in the gstat library

  • The covariance matrix is positive definite and has positive eigen values. It is singular according to gstat, but not to is.singular.matrix function

  • There were enough pair of points to do the experimental variogram

How to overcome this problem? What tips to avoid a singular covariance matrix? I also welcome any "best practice" for kriging.

Code (requires forSO.Rdata : https://www.dropbox.com/s/5vfj2gw9rkt365r/forSO.Rdata?dl=0 ) : 

library(ggplot2)
library(gstat)

#Attached Rdata
load("forSO.Rdata")

#The observations
str(abun)

#Spatial structure
abun %>% as.data.frame %>% 
  ggplot(aes(lon, lat)) +
  geom_point(aes(colour=prop_species_cells), alpha=3/4) + 
  coord_equal() + theme_bw()

#Number of pair of points
cvgm <- variogram(prop_species_cells ~1, data=abun, width=3,  cutoff=300)
plot(cvgm$dist,cvgm$np)

#Fit a model covariogram
efitted = fit.variogram(cvgm, vgm(model="Mat", range=100, nugget=1), fit.method=7, fit.sills=TRUE, fit.ranges=TRUE)
plot(cvgm,efitted)

#No warning, and the model is non singular
attr(efitted, "singular")

#Covariance matrix (only on a small set of points, I have more than 25000 points) : positive-definite, postiive eigen values and not singular
hex_pointsDegTiny=hex_pointsDeg
hex_pointsDegTiny@coords=hex_pointsDegTiny@coords[1:10,]
dists <- spDists(hex_pointsDegTiny)
covarianceMatrix=variogramLine(efitted, maxdist = max(cvgm$dist), n = 10*max(cvgm$dist), dir = c(1,0,0), dist_vector = dists, covariance = TRUE)
eigen(covarianceMatrix)$values
is.positive.definite(covarianceMatrix)
is.singular.matrix(covarianceMatrix)

# No duplicate locations
zerodist(hex_pointsDegTiny)

# Impossible to krig
OK_fit <- gstat(id = "OK_fit", formula = prop_species_cells ~ 1, data = abun, model = efitted)
dist <- predict(OK_fit, newdata = hex_pointsDegTiny)
dist@data
Solution

  • Actually, there were duplicate locations in abun dataset (zerodist(abun)), they were not to be seeked into the grid on which I wanted to krig estimates. After getting rid of the duplicates, kriging worked fine.