-----Original Message-----
From: Edzer Pebesma [mailto:edzer.pebesma at uni-muenster.de]
Sent: Tuesday, November 17, 2009 1:01 PM
To: Tomislav Hengl
Cc: 'r-sig-geo'
Subject: Re: [R-sig-Geo] '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
Tomislav Hengl wrote:
Hi Paul, Edzer,
I understand why the singular matrix problem happens and I know that there is not really a
mathematical solution around it:
x <- matrix(runif(100), nrow=10)
x.i <- solve(x)
str(x.i)
num [1:10, 1:10] 0.8191 -1.0293 0.0826 1.068 -0.2106 ...
# add a 'singular' column
x[,1:10] <- rep(1, 10)
x.i <- solve(x)
Error in solve.default(x) :
Lapack routine dgesv: system is exactly singular
However, I would very much support if you would integrate an "if" loop to check if it will
and then mask the prediction location using "NA".
I mean, at the moment every time we want to run predictions, even if only at a single prediction
location the model fails, we are not able to generate any output (this is sometimes really
frustrating).
Hence, I would support that you allow us to run predictions for all locations first, and then
-----Original Message-----
From: r-sig-geo-bounces at stat.math.ethz.ch [mailto:r-sig-geo-bounces at stat.math.ethz.ch] On
Of Paul Hiemstra
Sent: Tuesday, November 17, 2009 9:53 AM
To: Edzer Pebesma
Cc: r-sig-geo
Subject: Re: [R-sig-Geo] 'LDLfactor' error in 'krige' function
Edzer Pebesma wrote:
My guess is that from constant data, the (co)variance is constant and
zero, so the covariance matrix cannot be decomposed (hence: LDLfactor
errors).
Is this a case that autokrige should catch?
I added a check in autoKrige. The output for the example below is now:
kriging_result = autoKrige(zinc~1, meuse)
Error in autoKrige(zinc ~ 1, meuse) :
All data in attribute 'zinc' is identical and equal to 0
Can not interpolate this data
A new version of automap with this feature has been uploaded to CRAN.
cheers,
Paul
--
Edzer
Mauricio Zambrano wrote:
Dear List,
During some OK interpolations of daily precipitation, with the
'automap' library, I got the following error:
[using ordinary kriging]
"chfactor.c", line 130: singular matrix in function LDLfactor()
Error en predict.gstat(g, newdata = newdata, block = block, nsim = nsim, :
LDLfactor
The code I'm using works fine for other days, so I assume that the
distance between the measurement points is no the problem. When
looking at the data that rose the error, I realized that all the
measured values were equal to zero (I can not skip those days in which
all the measured points have the same value in advance, because the
measured value in those points change with time).
According to a traceback that is given below, it seems that the cause
is in the 'predict.gstat' function of the 'gstat' package.
The same error can be risen with:
# Data preparation
data(meuse)
coordinates(meuse) =~ x+y
data(meuse.grid)
gridded(meuse.grid) =~ x+y
meuse$zinc <- meuse$zinc*0
meuse$zinc
# Ordinary kriging, no new_data object
kriging_result = autoKrige(zinc~1, meuse)
traceback()
6: .Call("gstat_predict", as.integer(nrow(as.matrix(new.X))),
as.double(as.vector(raw$locations)),
as.vector(new.X), as.integer(block.cols), as.vector(block),
as.vector(bl_weights), as.integer(nsim), as.integer(BLUE))
5: predict.gstat(g, newdata = newdata, block = block, nsim = nsim,
indicators = indicators, na.action = na.action, debug.level =
debug.level)
4: .local(formula, locations, ...)
3: krige(formula, input_data, new_data, variogram_object$var_model,
block = block, ...)
2: krige(formula, input_data, new_data, variogram_object$var_model,
block = block, ...)
1: autoKrige(zinc ~ 1, meuse)
I would really appreciate any hint about the possible reason of this
error and if it could be overcome in some way.
Thanks in advance for any help.
Mauricio
--
Drs. Paul Hiemstra
Department of Physical Geography
Faculty of Geosciences
University of Utrecht
Heidelberglaan 2
P.O. Box 80.115
3508 TC Utrecht
Phone: +3130 274 3113 Mon-Tue
Phone: +3130 253 5773 Wed-Fri
http://intamap.geo.uu.nl/~paul