Roger Bivand wrote:
On Tue, 10 Nov 2009, ageel bushara wrote:
Dear list members,
I'm doing spatial interpolation of soil moisture using kriging with
external drift following the example given by Tom Hengl
(http://spatial-analyst.net/wiki/index.php?title=Regression-kriging_guide).
When I used one predictor (in my case gradient or cos_aspect, see
below) it works fine for me, but when I used multiple predictors
(gradient +cos_aspect), still I'm able to have the the experimental
variogram as well as the fitted theoretical variogram, however, I'm
not able to have the kriged soil moisture. here is the code
# import and define the input variables
library(sp)
library(lattice)
trellis.par.set(sp.theme()) # plots the final predictions using
blue-pink-yellow legend
water <- read.csv("water_content.txt") # table with x,y
coordinates,and z is? soil moisture
str(water)
coordinates(water)=~x+y?? # this makes depth a SpatialPointsDataFrame
str(water)
cos_aspect = read.asciigrid("cos_aspect.asc")? # reads ArcInfo Ascii
raster map
str(cos_aspect)
spplot(cos_aspect, scales=list(draw=T), sp.layout=list("sp.points",
water, pch="+"))
gradient= read.asciigrid("gradient.asc")
str(gradient)
spplot(gradient, scales=list(draw=T), sp.layout=list("sp.points",
water, pch="+"))
#Plot the xy graph target versus predictor:
cos_aspect.ov = overlay(cos_aspect, water)? # create grid-points overlay
str(cos_aspect...@data)
water$cos_aspect.asc = cos_aspect.ov$cos_aspect.asc?
gradient.ov = overlay(gradient, water)
water$gradient.asc= gradient.ov$gradient.asc
lm.depth <- lm(z~ gradient.asc+cos_aspect.asc, as.data.frame(water))
summary(lm.depth)
plot(z~ gradient.asc+cos_aspect.asc, as.data.frame(water))
abline(lm(z~ gradient.asc, as.data.frame(water)))
abline(lm(z~ cos_aspect.asc, as.data.frame(water)))
#Fit the variogram model of the residuals:
library(gstat)
null.vgm <- vgm(var(water$z), "Sph",
sqrt(areaSpatialGrid(cos_aspect))/4, nugget=0) #initial parameters
vgm_depth_r <- fit.variogram(variogram(z~ cos_aspect.asc+
gradient.asc, water), model=null.vgm)
plot(variogram(z~ gradient.asc+cos_aspect.asc,water), vgm_depth_r,
main="fitted by gstat")
# It works fine till here
# Run uk in gstat:
depth_uk <- krige(z~ cos_aspect.asc+gradient.asc, locations=water,
newdata=cos_aspect, model=vgm_depth_r)
the problem it seems mainly in newdata, here is the error given by R:
Error in eval(expr, envir, enclos) : object 'gradient.asc' not found.
if i assigned the? newdata= gradient then the R error message is :
Error in eval(expr, envir, enclos) : object 'cos_aspect.asc' not found
You probably need to think through what you are doing. You should use a
Spatial*DataFrame for newdata that includes "cos_aspect.asc" and
"gradient.asc" among the values returned by the names() method. Since
each of your two SpatialGridDataFrame objects have single variables, you
should review what a SpatialGridDataFrame is - a data.frame with a
description of the GridTopology. Assuming that your two objects share
their GridTopology values, maybe cbind() them? Then check that the
names() output includes the names of variables on the RHS in the
formula. Don't think of objects as grids, think of them as data.frames,
then you see what is going on.
Try with lm() and predict() first, then with gstat() and predict() - the
newdata= argument works in exactly the same way.
Hope this helps,
Roger
Hi Ageel,
Roger is right. You should read all rasters to a single (multilayer)
grid e.g.:
> grids <- readGDAL("cos_aspect.asc")
> names(grids) <- "cos_aspect.asc"
> grids$gradient.asc <- readGDAL("gradient.asc")$band1
...
Then overlay, filter the NA's, and then you can fit a variogram and run
predictions e.g.:
depth_uk <- krige(z~ cos_aspect.asc+gradient.asc, locations=water,
newdata=grids, model=vgm_depth_r)
...
Here are some more examples:
http://geomorphometry.org/content/some-examples-rsaga
http://spatial-analyst.net/scripts/meuse.R
all the best,
T. Hengl
http://home.medewerker.uva.nl/t.hengl/
as I understand that newdata is the grid where it does the interpolations
How can I obtained the kriged soil moisture?, what is wrong here?,
Thanks,
Ageel
[[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
