[R-sig-Geo] Problem in regression kriggin

[halim10-fes]

Dear Group,

Happy New Year to Everybody.

I am a novice user of gstat and geostatistics. I have a data set (please see 
attached: data.csv) containing coordinates (East_x, North_y),Elev(Elevation), 
and Richness (Plant Species) collected from 125 plots of 10m x 10m plots. I 
want to interpolate plant species richness from the relationship between 
Richness and Elevation of the study area using Regression Krigging. I have 
written the following codes:

#Variogram of residuals
#===========================================================================

data<-read.csv("E://data.csv",header=TRUE)
attach(data)
names(data)
m2=glm(Richness~Elev)
summary(m2)
residuals<-m2$residuals
data$Residuals<-residuals

library(maptools)

coordinates(data) = ~ East_x+North_y

library(gstat)

var <- variogram(Residuals~1, data=data,alpha=c(0,22.5,45,67.5,90,112.5,135),
                 tol.hor=30,cutoff=1000,width=67)
plot(var)

#Variogram shows that major axis is 45 & minor 112.5; Fitting variogram in 
#major axis

var.45<-var$dir.hor==45
plot(var[var.45,])

eye.sph<-vgm(psill=10.5,model="Sph",range=705,nugget=2)
plot(var[var.45,],eye.sph)


#Fit anisotropic variogram

fit.sph<-fit.variogram(var,model=vgm(psill=10.5,model="Sph",range=705,
                                     nugget=2,anis=c(45,0.72)))
plot(var, model=fit.sph, as.table=TRUE)

fit.exp<-fit.variogram(var, model=vgm(psill=10.5,model='Exp',range=705/3,
                                      nugget=2, anis=c(45, 0.72)))
plot(var, model=fit.exp, as.table=TRUE)

fit.gau<-fit.variogram(var,model=vgm(psill=10.5,model='Exp',range=705/sqrt(3),
                                     nugget=2, anis=c(45, 0.72)))
plot(var, model=fit.gau, as.table=TRUE)

#Cross Variogram

cv.sph<-krige.cv(Residuals~1,data,fit.sph)
View(cv.sph)
mean(cv.sph$zscore) #-0.0001968425
var(cv.sph$zscore)  #1.114187

cv.exp<-krige.cv(Residuals~1,data,fit.exp)
View(cv.exp)
mean(cv.exp$zscore)  #-1.243423e-05
var(cv.exp$zscore)   #1.116687

cv.gau<-krige.cv(Residuals~1,data,fit.gau)
View(cv.gau)
mean(cv.gau$zscore)  #-1.235693e-05
var(cv.gau$zscore)   #1.116684


#Cross variogram shows that spherical model shows the best fit. The sill 
varies in different direction of the varigram. So, we will have to fit a zonal 
#anisotropic variogram with spherical model (fit.sph). Variogram shows that
#22.5 deg has the highest sill and 112.5 deg should have the lowest sill. We 
#will fit the zonal variogram to the 112.5 deg direction with an anisotropic 
#ratio of 293/(1600*10000)

#Fitting zonal anisotropic variogram

fit.sph1<-vgm(psill=2.5,"Sph",range=1600*10000,anis=c(122.5,1.825e-05),
              add.to=fit.sph)
plot(var,fit.sph1,as.table=TRUE)

#Cross Variogram of the final model

cv.sph1<-krige.cv(Residuals~1,data,fit.sph1)
View(cv.sph1)
mean(cv.sph1$zscore)  #-0.0009952294
var(cv.sph1$zscore)   #1.071403

#Creating a grid file
#============================================================================

grd <- expand.grid(x=seq(from=341960,to=343195,by=10), 
                   y=seq(from=2667253,to=2668873,by=10))

library(maptools)

coordinates(grd) <- ~ x+y
gridded(grd) <- TRUE

So far,I think, things were going well but I found a error while running the 
following code

#Regression kriggin
#===========================================================================

prk.rich<-krige(Richness~Elev, data=data, newdata=grd,model=fit.sph1)

Error in function (classes, fdef, mtable)  : 
  unable to find an inherited method for function ‘krige’ for signature 
‘"formula", "missing"’

Can anyone help me out of this please? I shall be grateful.

Best regards,

---------------
Md. Abdul Halim
Assistant Professor
Department of Forestry and Environmental Science
Shahjalal University of Science and Technology,Sylhet-3114,
Bangladesh.
Cell: +8801714078386.
alt. e-mail: xou03 at yahoo.com

-- 
This message has been scanned for viruses and
dangerous content by MailScanner, and is
believed to be clean.

-------------- next part --------------
A non-text attachment was scrubbed...
Name: data.csv
Type: application/vnd.ms-excel
Size: 4212 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20130108/36efe1ae/attachment.xlb>