[R-sig-Geo] Multiple predictors for external drift kriging

[Roger Bivand]

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

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

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no