[R-sig-Geo] help with Ordinary Kriging with R (singular matrix error)

[Massimo Di Stefano]

Hello All,

I'm tring to perfotm an ordinary Kriging interpolation on aset of data.
It's my 'first real attempt' with R and Kriging, so please apologize my
mistake ... i'm slowly tring to learn it.

The dataset i'm working on is a shp file, for convenience (to better
follow the notes i find on [1] ) i converted it to sc.csv (attached
there [2]).


I tried using the following code in R :




"""
library(gstat)
d <- read.csv('sc.csv')
e <- na.omit(d)
coordinates(e) <- ~ x+y
bubble(e, zcol='scallop_co', fill=FALSE, do.sqrt=FALSE, maxsize=2)
x.range <- as.integer(range(e at coords[,1]))
y.range <- as.integer(range(e at coords[,2]))
grd <- expand.grid(x=seq(from=x.range[1], to=x.range[2], by=500),
y=seq(from=y.range[1], to=y.range[2], by=500) )
coordinates(grd) <- ~ x+y
gridded(grd) <- TRUE
plot(grd, cex=0.5)
points(e, pch=1, col='red', cex=0.7)
title("Interpolation Grid and Sample Points")

g <- gstat(id="elev", formula=scallop_co ~ 1, data=e)
plot(variogram(g, map=TRUE, cutoff=4000, width=200), threshold=10)
v <- variogram(g, alpha=c(0,45,90,135))
v.fit <- fit.variogram(v, model=vgm(model='Lin' , anis=c(0, 0.5)))
plot(v, model=v.fit, as.table=TRUE)
g <- gstat(g, id="scallop_co", model=v.fit )

## perform ordinary kriging prediction:
p <- predict(g, model=v.fit, newdata=grd)
"""

It gived me the error :

> p <- predict(g, model=v.fit, newdata=grd)
[using ordinary kriging]

"chfactor.c", line 130: singular matrix in function LDLfactor()
Error in predict.gstat(g, model = v.fit, newdata = grd) : LDLfactor


it is probably that i have to change some parameters in the previouse
code, i'm tring to learn but i've not yet any clue on how to fix it.


i tried to add more value in when i create the directional variogram :

"""
v <- variogram(g, alpha=c(0,10,20,30,40,45,50,60,70,80,90,100,110,135))

"""
the results is attached at [3]
 
my sessionInfo is :

"""
> sessionInfo()
R version 2.12.0 (2010-10-15)
Platform: x86_64-pc-linux-gnu (64-bit)

locale:
 [1] LC_CTYPE=en_US.utf8       LC_NUMERIC=C             
 [3] LC_TIME=en_US.utf8        LC_COLLATE=en_US.utf8    
 [5] LC_MONETARY=C             LC_MESSAGES=en_US.utf8   
 [7] LC_PAPER=en_US.utf8       LC_NAME=C                
 [9] LC_ADDRESS=C              LC_TELEPHONE=C           
[11] LC_MEASUREMENT=en_US.utf8 LC_IDENTIFICATION=C      

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods
base     

other attached packages:
[1] gstat_0.9-74 sp_0.9-74   

loaded via a namespace (and not attached):
[1] grid_2.12.0     lattice_0.19-13
> 
"""

thankyou a lot for any help!

thanks,

Massimo.



[1] http://casoilresource.lawr.ucdavis.edu/drupal/node/442
[2] http://www.geofemengineering.it/data/sc.csv
[3] http://www.geofemengineering.it/data/sc_variogram.png