[R-sig-Geo] Problem with universal kriging using gstat

[Tomislav Hengl]

Auxiliary variables that are used to explain the trend-part of variation need to be available also
at all new prediction locations.

see ?krige:

"newdata - data frame or Spatial object with prediction/simulation locations; should contain
attribute columns with the independent variables (if present) and (if locations is a formula) the
coordinates with names as defined in locations"

Gstat obviously can not find X and Y coordinates of your new locations (data.grid).

Try instead:

> data.grid <- expand.grid(x=seq(xLL,xUL,xN), y=seq(xLL,xUL,xN))
> names(data.grid) = c("X","Y")
> gridded(data.grid) <- ~X+Y

cheers,

Tom Hengl
http://spatial-analyst.net/ 


-----Original Message-----
From: r-sig-geo-bounces at stat.math.ethz.ch [mailto:r-sig-geo-bounces at stat.math.ethz.ch] On Behalf Of
Yong Li
Sent: maandag 14 april 2008 6:45
To: r-sig-geo at stat.math.ethz.ch
Subject: [R-sig-Geo] Problem with universal kriging using gstat

Hi ALL,
 
Can any expert to see my mistake in the following R script when I try the gstat universal kriging?
The error occurs at the last step. It sounds I missed the setup for the names of coordinates of the
SpatialGrid "data.grid".
 
Regards
 
Yong
 
> memory.size()
[1] 47930616
> memory.size(TRUE)
[1] 85377024
> memory.limit(size=2048)
NULL
> round(memory.limit()/1048576.0, 2)
[1] 2048
> options(scipen=3)
> 
> options(digits=15)
> rm(list=ls())
> require(sp)
[1] TRUE
> require(gstat)
[1] TRUE
> require(maptools)
[1] TRUE
> require(foreign)
[1] TRUE
> 
> cellsize <- 10
> 
> #                 1           2               3             4
> VarName=c("logit(SOM)","logit(TN)","logit(OLSENP)","logit(EXTK)")
> out=paste("Variable","No","Model","C0","C","Spatial","Range","R2","p","SSE",sep=",")
> setwd("E:\\Rcode\\EJSS\\Hongtong\\village")
> data <- read.dbf("libaocun_gs.dbf")
> names(data)
 [1] "ID"       "XCOOR"    "YCOOR"    "NO"       "CDBH"     "NO3N"     "NH4N"     "MINERALN" "SOM"
"TOTALN"   "OLSENP"   "EXTK"    
[13] "TN"       "X"        "Y"       
> class(data)
[1] "data.frame"
> names(data)[14] <- "X"
> names(data)[15] <- "Y"
> 
> data$X <- data$X  + 1750
> data$Y <- data$Y +  650
>  
> coordinates(data) <- ~X+Y
> proj4string(data) <- CRS("+proj=tmerc +lat_0=0 +lon_0=111 +k=1.0 +x_0=500000 +y_0=0 +ellps=krass
+units=m +no_defs +towgs84=22,-118,30.5,0,0,0,0")
> 
> #read in a shape file of the boundary
> aSHAPE <- readShapePoly("libao_bd.shp", IDvar="BSM", proj4string=CRS(proj4string(data))) 
> 
> #generate the grid coordinates (in meters)
> xLL <- round(slot(aSHAPE,"bbox")[1]) - 100
> xUR <- round(slot(aSHAPE,"bbox")[3]) + 100
> yLL <- round(slot(aSHAPE,"bbox")[2]) - 100
> yUR <- round(slot(aSHAPE,"bbox")[4]) + 100
> xN <- round((xUR-xLL)/cellsize)
> yN <- round((yUR-yLL)/cellsize)
> data.grid = SpatialGrid(GridTopology(c(xLL,yLL),c(cellsize,cellsize),c(xN,yN)))
> proj4string(data.grid) <- CRS(proj4string(aSHAPE))
> 
> 
> #data logit transformation
> SOM_MAX      <- max(data$SOM)*1.1
> SOM_MIN      <- min(data$SOM)*0.9
> SOM_DIFF     <- SOM_MAX - SOM_MIN
> data$SOMt    <- log(((data$SOM-SOM_MIN)/SOM_DIFF)/(1-(data$SOM-SOM_MIN)/SOM_DIFF))
> 
> TN_MAX      <- max(data$TN)*1.1
> TN_MIN      <- min(data$TN)*0.9
> TN_DIFF     <- TN_MAX - TN_MIN
> data$TNt    <- log(((data$TN-TN_MIN)/TN_DIFF)/(1-(data$TN-TN_MIN)/TN_DIFF))
> 
> OLSENP_MAX      <- max(data$OLSENP)*1.1
> OLSENP_MIN      <- min(data$OLSENP)*0.9
> OLSENP_DIFF     <- OLSENP_MAX - OLSENP_MIN
> data$OLSENPt    <-
log(((data$OLSENP-OLSENP_MIN)/OLSENP_DIFF)/(1-(data$OLSENP-OLSENP_MIN)/OLSENP_DIFF))
> 
> EXTK_MAX      <- max(data$EXTK)*1.1
> EXTK_MIN      <- min(data$EXTK)*0.9
> EXTK_DIFF     <- EXTK_MAX - EXTK_MIN
> data$EXTKt    <- log(((data$EXTK-EXTK_MIN)/EXTK_DIFF)/(1-(data$EXTK-EXTK_MIN)/EXTK_DIFF))
> 
> #SOM
>    # plot
>    X11(width=8, height=6.5)
>    par(mfrow=c(1,1))
>    spplot(aSHAPE["BSM"], scales=list(draw=TRUE, cex=0.7), xlab="Easting (m)", ylab="Northing (m)",
+           sp.layout=list("sp.points", pch="+", col="black", data), col.regions=FALSE,
colorkey=FALSE)
> 
>    X11(width=8, height=6.5)
>    par(mfrow=c(1,2))
>    hist(data$SOMt, main="", xlab=VarName[1], n=25)
>    boxplot(data$SOMt, xlab=VarName[1])
> 
>    #estimate semivariogram model
>    soil.v <- variogram(SOMt~X+Y+X*Y+X^2+Y^2, data=data, cutoff=1050, width=90)
>    vmodel <- vgm(0.5, "Sph", 1000, 0.1)
>    v.sph <- fit.variogram(object=soil.v, model=vmodel, fit.sills=TRUE, fit.ranges=TRUE,
fit.method=7)
>    v.sph
  model             psill            range
1   Nug 0.269063860535796   0.000000000000
2   Sph 0.217370623520010 815.994438576304
>    X11(width=8, height=6.5)
>    par(mfrow=c(1,1))
>    plot(soil.v, model=v.sph, col="red", cex=1, lwd=2.0, main=VarName[1], xlab="Separation distance
(m)")
> 
> #Univeral kriging
> data.gstat <- gstat(id="SOMt", formula=SOMt~X+Y+X*Y+X^2+Y^2, data=data, nmax=15, model=v.sph)
> 
> #Predicting surface
> data.grid <- predict(data.gstat, data.grid)
Error in eval(expr, envir, enclos) : object "X" not found
> 
> 


	[[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo