[STATSGRASS] gstat: how to do something like this

[Edzer J. Pebesma]

Jaroslaw, in R this example is obtained as follows. Assume part.a is a 
field in meuse.grid that tells part a apart from part b:

require(gstat)
data(meuse)
coordinates(meuse) = ~x+y
data(meuse.grid)
gridded(meuse.grid) = ~x+y

# now get the partitioning for meuse:
meuse$part.a = meuse.grid$part.a[overlay(meuse.grid, meuse)]

# next do the kriging for each sub-region:
x1 = krige(log(zinc)~1, meuse[meuse$part.a == 0,],
meuse.grid[meuse.grid$part.a == 0,], model = vgm(.548, "Sph", 900, .0654),
nmin = 20, nmax = 40, maxdist = 1000)

# next...
x2 = krige(log(zinc)~1, meuse[meuse$part.a == 1,],
meuse.grid[meuse.grid$part.a == 1,], model = vgm(.716, "Sph", 900),
nmin = 20, nmax = 40, maxdist = 1000)

# combine results...
x = rbind(as.data.frame(x1), as.data.frame(x2))
gridded(x) = ~x+y
spplot(x["var1.pred"], main = "stratified kriging predictions")

Hth,
--
Edzer


Jarosław Jasiewicz wrote:
> Hi
>
>
> It means deffernts models for different areasin one grid. Now I can do 
> this with different areas (next merged in grass). As I look on meuse 
> dataset it is probably connected with factor in meuse.grid data frame. 
> I need this to analize different geomorphological units.
>
> regards
> Jarek
>
>
> ------------------------------------------------------------------------
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> statsgrass mailing list
> statsgrass at grass.itc.it
> http://grass.itc.it/mailman/listinfo/statsgrass
>