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

[Tomislav Hengl]

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