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