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

[Bruin, Sytze de]

On 13/01/16 15:01, 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]<http://r-sig-geo.2731867.n2.nabble.com/user/SendEmail.jtp?type=node&node=7589380&i=0> <mailto:[hidden email]<http://r-sig-geo.2731867.n2.nabble.com/user/SendEmail.jtp?type=node&node=7589380&i=1>>>
> 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]<http://r-sig-geo.2731867.n2.nabble.com/user/SendEmail.jtp?type=node&node=7589380&i=2> <mailto:[hidden email]<http://r-sig-geo.2731867.n2.nabble.com/user/SendEmail.jtp?type=node&node=7589380&i=3>>
>     > 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<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


	[[alternative HTML version deleted]]