'LDLfactor' error in 'krige' function

Tom, you can already do this:

library(gstat)
data(meuse)
coordinates(meuse)=~x+y
data(meuse.grid)
gridded(meuse.grid)=~x+y
meuse=meuse[c(1,1:155),] # replicate first observation
pr1 = predict(zinc~1,meuse,meuse.grid,vgm(1, "Exp", 300)) # will break
# the following will generate NA's for cells where the estimated
condition number of the covariance matrix exceeds 1e10:
pr1 = krige(zinc~1,meuse,meuse.grid,vgm(1, "Exp",
300),nmax=30,                 set=list(cn_max=1e10))
# generate a full interpolation as comparison:
pr1$idw = idw(zinc~1,meuse,meuse.grid)[[1]]
# show side by side:
spplot(pr1, c("var1.pred", "idw"), col.regions=bpy.colors())

... missing values are generated for those cells that result in a
singular covariance matrix, given their local neighbourhood of 30
nearest observations.

cn_max referers to the maximum allowed condition number, see
http://en.wikipedia.org/wiki/Condition_number

Two issues are (i) what singularity means when we use floating point
representations for real numbers, and (ii) that condition numbers of
matrices are expensive to evaluate, and an estimate based on the LU
decomposition is used. I threshold to 1e10 here, but that is purely for
illustrational purposes.

Note, for you of interest, that IIRC this thresholding is also done for
singularity of the X matrix, holding the predictors.

All this information didn't make it to the help pages of the krige
function. This help page would cover tens of pages otherwise. For full
documentation, still the "old", gstat stand-alone manual on gstat.org is
needed. I agree that this is not optimal.

With best wishes,
--
Edzer

'LDLfactor' error in 'krige' function

Thread (11 messages)