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

[Massimo Di Stefano]

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