[R-sig-Geo] 'LDLfactor' error in 'krige' function

Tue Nov 17 15:01:17 CET 2009

Hi Edzer,

Thanks for the info. I was not aware that I can simply set cn_max and the predictions will not break
(and I still do read the gstat manual).

I guess that this is then the solution to our problem. 

For me it enough to generate a map with NA's, then zoom into map to see where the problematic
areas/points are, then either filter the maps, consider using different search radius or threshold.

T. Hengl
http://home.medewerker.uva.nl/t.hengl/

> -----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
> happen,
> > 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
> let
> > the users judge if one should increase the search radius, remove some too-uniform predictors etc
> > etc.
> >
> > I hope you agree with me,
> >
> > T. Hengl
> > http://home.medewerker.uva.nl/t.hengl/
> >
> >
> >
> >> -----Original Message-----
> >> From: r-sig-geo-bounces at stat.math.ethz.ch [mailto:r-sig-geo-bounces at stat.math.ethz.ch] On
> Behalf
> >> 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
> >>
> >> _______________________________________________
> >> R-sig-Geo mailing list
> >> R-sig-Geo at stat.math.ethz.ch
> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >>
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > R-sig-Geo at stat.math.ethz.ch
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
> 
> --
> Edzer Pebesma
> Institute for Geoinformatics (ifgi), University of Münster
> Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
> 8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de/
> http://www.springer.com/978-0-387-78170-9 e.pebesma at wwu.de