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@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
