[R-sig-Geo] help with Ordinary Kriging with R (singular matrix error)
Edzer Pebesma
edzer.pebesma at uni-muenster.de
Wed Dec 15 21:32:29 CET 2010
The default is global kriging, meaning it wants to solve a system of
25000 x 25000. To avoid that, specify nmax to autokrige:
autokrige(..., nmax = 100) # or so
On 12/15/2010 08:58 PM, Massimo Di Stefano wrote:
> 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
>>
>>
>>
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of Münster
Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de
http://www.52north.org/geostatistics e.pebesma at wwu.de
More information about the R-sig-Geo
mailing list