help with Ordinary Kriging with R (singular matrix error)

I tried to go ahead to fix my problem with R and Kriging ..

as the example i tried produce a bad varuiogram fitting (maybe my fault)
i tried the package automap using the following sintax :



d <- read.csv('/home/habcam/scl.csv')
e <- na.omit(d)
coordinates(e) <- ~ x+y
variogram = autofitVariogram(scallop_de~1,e)
plot(variogram)
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=300),
y=seq(from=y.range[1], to=y.range[2], by=300) )
gridded(grd) =~ x+y

kriging_result = autoKrige(scallop_de ~ 1, e, grd)
Warning in autoKrige(scallop_de ~ 1, e, grd) :
  Removed 2018 duplicate observation(s) in input_data:
              coordinates   cat scallop_de
2       (677919, 4546180)     2       2.05
20      (678847, 4545330)    20       2.02
59      (682175, 4541480)    59       0.81
...
...
...
21808   (684423, 4591090) 21808       1.35
21810   (684417, 4591090) 21810       1.71
21815   (684396, 4591090) 21815       1.53
21821   (684346, 4591090) 21821       1.06
21827   (684249, 4591090) 21827       2.66
21948   (699146, 4595160) 21948       1.42
[using ordinary kriging]




but R start to think ... and nver ends, after more than 1 hour i killed
the process :-( 

have you any suggestion on how can i "krige" my data [1] ?

 [1] files.me.com/epiesasha/ynl85n



i also tried a different approach using :

mydata <- read.csv("file.csv")
library(sp)
coordinates(mydata) <- c("x","y")
library(rgdal)
writeOGR(obj=mydata,driver="ESRI Shapefile",dsn=".",layer="mydata")


library(gstat)
library(rgdal)
mydata <- readOGR(dsn=".",layer="mydata")
plot(variogram(mydata at data$scallop_de~1,mydata, cloud=TRUE))

the last line produce an R crash ... after freezed the computer.


are 25000 points too mutch data for R ?



On Tue, 2010-12-14 at 15:34 -0500, Massimo Di Stefano wrote:

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