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

[Paul Hiemstra]

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
>   
Try to check before running the interpolation if all the observations 
are zero. If this is the case, then instead of kriging return a grid 
where all the values are constant. Most probably this only occurs with 
all zeros and with e.g. all 4.67. To check if all data is the same use a 
cod structure like:

if(length(unique(meuse$zinc)) != 1) {
    autoKrige
} else {
    return grid with only zeros (or another value is that is constant
}

cheers,
Paul
> # 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