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

[Tomislav Hengl]

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

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