[R-sig-Geo] cross variogram

Tom Gottfried tom.gottfried at tum.de
Sun Jan 1 15:11:12 CET 2012


Hi,

Am 30.12.2011 18:45, schrieb Saman Monfared:
> Hi.
> I want fit a cross variogram with two variables. but it is non-possitive
> definite.
> even if I use fit.range=T, prediction is performed by inverse distance not
> cokriging and all of  variance prediction are NA.
> my program is like below.
>
> rm(list=ls())
> library(gstat)
> dd<-read.table("bees.txt",header=TRUE)
> x.range<- as.integer(range(dd at coords[,1]))
> y.range<- as.integer(range(dd at coords[,2]))
> grd<- expand.grid(x=seq(from=x.range[1]-20000, to=x.range[2]+90000,
> by=500), y=seq(from=y.range[1]-150000, to=y.range[2]+90000, by=500) )
> coordinates(grd)<- ~ x+y
> gridded(grd)<- TRUE
> plot(grd, cex=0.5)
> points(dd, pch=3, col='blue', cex=0.7)
> coordinates(dd)<- ~ x+y
> m<-vgm(1954,"Gau",80017,225)
> g<- gstat(NULL, id = " B", form =  B ~ 1,data=dd)
> g<- gstat(g, id = "E", form = E~ 1,data=dd)
> g<- gstat(g, id = "T", form =  T~ 1,data=dd)
> vm<- variogram(g)
> plot(vm)
> vm.fit<- fit.lmc(vm, g, model=m,fit.ranges =T)

?fit.lmc says "returns an object of class gstat, with fitted 
variograms". Thus vm.fit is your gstat-object with the fitted models. In 
contrast, the object g still is without any model.

> plot(vm, vm.fit,col="black",type="b",lwd=3,main="Fig1")
>   cv.c2<- gstat.cv(g)

use `gstat.cv(vm.fit)' to get gstat using cokriging. But beware, the 
default of the argument fit.lmc in fit.lmc! With fit.ranges=T you will 
probably get non-positive definite coefficient matrices.

HTH,
Tom

> summary(cv.c2$residual)
> sqrt(mean(cv.c2$residual^2))
> mean(cv.c2$residual)
> k.c2<- predict.gstat(g, dd)
> image(k.c2)
> summary(k.c2)
> contour(k.c2,add=T)
>
> Thanks for your help.
>

-- 
Technische Universität München
Department für Pflanzenwissenschaften
Lehrstuhl für Grünlandlehre
Alte Akademie 12
85350 Freising / Germany
Phone: ++49 (0)8161 715324
Fax:   ++49 (0)8161 713243
email: tom.gottfried at wzw.tum.de
http://www.wzw.tum.de/gruenland



More information about the R-sig-Geo mailing list