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

[Mauricio Zambrano]

Dear Paul and Edzer,


2009/11/17 Tomislav Hengl <hengl at spatial-analyst.net>:
>
> 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).

I agree with Tom that it may be really frustrating not to get any
output, even most when this situation happens within a loop that
involves many time steps. I know that the user can easily get rid of
this manually for a single prediction, but it should be better if the
check could be integrated into the 'krige' or the 'autoKrige'
function, but of course this is your decision.


>
> 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.

In case you could allow us to run the predictions first, it shouldn't
be better to mask the prediction location with the constant value that
caused the error instead of NA ?. Or to have a flag that allow the
user to choose if to get an NA or the constant value ?


Kinds,

Mauricio

-- 
============================================
Ph.D. Candidate,
Dept. of Civil and Env. Engineering
University of Trento, Italy
=============================================
Linux user #454569 -- Ubuntu user #17469
=============================================
"Planning is bringing the future into the present
so that you can do something about it now"
(Alan Lakein)

>
> 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
>