[R-sig-Geo] Universal Block Kriging covariate definition

[Bruin, Sytze de]

Antonio Manuel Moreno Ródenas wrote: 

> Thanks a lot Edzer, 
> 
> I'm not sure that would work. 
> In that way I would transfer to the kriging function the averaged value 
> of the covariate in the block. I'm not sure that would make the kriging 
> behave correctly. 
> 
> All the points calculated with the prediction inside the block (and 
> later averaged to give the block kriging prediction) will have as 
> "drift" the average of the covariate in the block. Instead of getting 
> the correct spatial variability inside the block (given by the 
> covariate). At first sight it doesn't seems correct to me. Am I wrong? 
I think so, when the operations (computing the drift, and block 
averaging) are both linear, it does not matter in which order they are 
carried out: f(g(x)) = g(f(x)). 

It would be easy to verify by computing the universal point kriging 
values and aggregating those. Try 

> library(sp) 
> demo(meuse, ask = FALSE, echo = FALSE) 
> library(gstat) 
>
> v = vgm(.5, "Sph", 900, .1) 
> kr1 = krige(log(zinc)~dist, meuse, meuse.grid, v) 
[using universal kriging] 
>
> meuse.area$dist = aggregate(meuse.grid["dist"], meuse.area)[[1]] 
> kr2 = krige(log(zinc)~dist, meuse, meuse.area, v) 
[using universal kriging] 
> kr2$kr1 = aggregate(kr1["var1.pred"], meuse.area)[[1]] 
>
> kr2$var1.pred 
[1] 5.687753 
> kr2$kr1 
[1] 5.685026 
>
> kr2$var1.pred / kr2$kr1 
[1] 1.00048 
>

My guess is that the difference can be attributed to how the area is 
discretized (see ?predict.gstat) 

> 
> kind regards 
> 
> Antonio Manuel Moreno Rodenas 
> 
> /Marie Curie Early Stage Researcher/ 
> /PhD Candidate/ 
> 
> *T**U **Delft / Section Sanitary Engineering, office 4.64*____ 
> 
> *Civil Engineering and Geoscience Faculty * 
> 
> T +31 15 278 14 62 
> 
> 
> On 13 January 2016 at 14:43, Edzer Pebesma 
> <[hidden email] <mailto:[hidden email]>>
> wrote: 
> 
> 
> 
>     On 13/01/16 14:16, Antonio Manuel Moreno Ródenas wrote: 
>     > Hello, I would like to rise a question on the use of predict {gstat}, 
>     > 
>     > I'm trying to perform the estimation of a spatially distributed variable at 
>     > the support scale of a particular area (Block kriging). I have access to an 
>     > additional variable, it is known that the variable of interest is 
>     > correlated to the new variable. So I would be interested on updating my 
>     > estimation by the use of this new information. This could be done by the 
>     > use of a kriging with external drift (KED), but with a block support 
>     > (Universal Block kriging). Theoretically this is included in the gstat 
>     > library as mentioned in the documentation. 
>     > 
>     > The issue comes when I try to perform the prediction: 
>     > 
>     > blockprediction <- predict(gstat(formula=Variabletopredict~additionalVariable, 
>     > data=Observed, model=vgm), newdata = shapefile) 
>     > 
>     > The newdata argument should contain the prediction location. In a normal 
>     > KED we would include a dataframe with a grid (coordinates in which to 
>     > predict) and the values of the covariate (additionalVariable). As I'm 
>     > trying to use a universal block kriging, I understood the newdata should be 
>     > the region in which I'm interested to know the prediction, hence a polygon. 
>     > How could I include in newdata the values of the covariate if its 
>     > resolution is finer than my block? 
>     > 
>     > As far as I know, what block kriging does is to predict point values inside 
>     > the region (which I could specified with the argument sps.args 
>     > discretization), and later average them. But I don't know how to attach the 
>     > covariate values to the block of interest (shapefile). 
> 
>     maybe by 
> 
>     shapefile = aggregate(additionalVariable, shapefile, mean) 
> 
>     > 
>     > Thanks in advance, 
>     > I hope I could explain it properly, but I will give more details if 
>     > necessary. 
>     > Kind regards, 
>     > Antonio 
>     > 
>     >       [[alternative HTML version deleted]] 
>     > 
>     > _______________________________________________ 
>     > R-sig-Geo mailing list 
>     > [hidden email] <mailto:[hidden email]>
>     > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>     > 
> 
>     -- 
>     Edzer Pebesma 
>     Institute for Geoinformatics  (ifgi),  University of Münster 
>     Heisenbergstraße 2, 48149 Münster, Germany; +49 251 83 33081 
>     <tel:%2B49%20251%2083%2033081> 
>     Journal of Statistical Software:   http://www.jstatsoft.org/
>     Computers & Geosciences:   http://elsevier.com/locate/cageo/
>     Spatial Statistics Society http://www.spatialstatistics.info
> 
> 

I believe the residual variogram should then be computed using the covariate data at block support. 

Sytze de Bruin
Wageningen University
Laboratory of Geo-Information Science and Remote Sensing