An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20080414/389cd434/attachment.pl>
Problem with universal kriging using gstat
3 messages · Yong Li, 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
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
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
_______________________________________________ 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