[R-sig-Geo] Problem with universal kriging using gstat
Tomislav Hengl
hengl at science.uva.nl
Mon Apr 14 13:17:04 CEST 2008
sorry, a type error:
> cellsize <- 10 # see also http://spatial-analyst.net/pixel.php
> xLL <- round(slot(aSHAPE,"bbox")[1])+cellsize/2
> xUR <- round(slot(aSHAPE,"bbox")[3])-cellsize/2
> yLL <- round(slot(aSHAPE,"bbox")[2])+cellsize/2
> yUR <- round(slot(aSHAPE,"bbox")[4])-cellsize/2
> xN <- round((xUR-xLL)/cellsize)
> yN <- round((yUR-yLL)/cellsize)
> data.grid <- expand.grid(x=seq(xLL,xUR,xN), y=seq(yLL,yUR,yN))
> names(data.grid) = c("X","Y")
> gridded(data.grid) <- ~X+Y
-----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
Tomislav Hengl
Sent: maandag 14 april 2008 10:37
To: 'Yong Li'; r-sig-geo at stat.math.ethz.ch
Subject: Re: [R-sig-Geo] Problem with universal kriging using gstat
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
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
More information about the R-sig-Geo
mailing list