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.ov@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

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]]

