[R-sig-Geo] Universal Block Kriging with SpatialPolygons

[French, Joshua]

Is it possible to use the predict or gstat function in the gstat package to do universal block kriging with SpatialPolygons where the universal kriging covariates are spatial coordinates?  The relevant covariates must be available for any point within the SpatialPolygon, which effectively limits to covariates to being the spatial coordinates.  Even with the restriction of coordinate-based covariates, I can't get this to work for *universal* kriging (it works fine for ordinary kriging).

The example below illustrates the problem using a reproducible example.  I think the issue is that the coordinates within each Polygons object inside the SpatialPolygons object do not having coordinate names matching those of the observed data set.  However, I can't seem to figure out a way to give them the desired name.  Is there a simple way to overcome this problem?  I do have a hack to get the predictions I need, but it certainly feels like there is a simpler approach that eludes me.

Thanks!

library(gstat)
library(sp)
data(meuse)
coordinates(meuse) = ~x+y
g = gstat(formula = log(zinc) ~ x + y, data = meuse, model = vgm(1, "Sph", 300, 1))

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

# typical prediction on a grid (no problem)
p = predict(g, newdata = meuse.grid)

# create a SpatialPolygons object for prediction
sp_polys = SpatialPolygons
p1 = cbind(c(179000, 179000, 179500, 179600),
           c(330000, 331000, 331000, 330000))
Poly1 = Polygon(p1)
p2 = cbind(c(179600, 179500, 180100, 180000),
           c(330000, 331000, 331000, 330000))
Poly2 = Polygon(p2)
ps1 = Polygons(list(Poly1), 1)
ps2 = Polygons(list(Poly2), 2)
sps = SpatialPolygons(list(ps1, ps2))

predict(g, newdata = sps)
# Error in eval(predvars, data, env) : object 'x' not found
--
Joshua French
Department of Mathematical and Statistical Sciences
University of Colorado Denver


	[[alternative HTML version deleted]]