[R-sig-Geo] overlay spatialgrid on spatial polygons data frame

[Roger Bivand]

On Thu, 20 May 2010, Agustin Lobo wrote:

> I have an SpPolDF and a coarse SpatialGrid, i.e.,
> srdf=SpatialPolygonsDataFrame(sr, data.frame(cbind(1:3,5:3),
> row.names=c("r1","r2","r3")))
> #the one used in help(overlay);
>> srdf at data
>   X1 X2
> r1  1  5
> r2  2  4
> r3  3  3
>
>> bbox(srdf)
>     min    max
> x 178544 181477
> y 329665 333676
>
>> test <- SpatialGrid(GridTopology(cellcentre.offset=c(178440,329600), cellsize=c(500,500), cells.dim=c(8,10)))
>> plot(srdf)
>> plot(test,add=T)
>
> My goal is to "overlay" test over srdf and tranfer the data in
> srdf at data to a data frame associated
> to test (=> test must become a SpatialGridDF) so that each cell has
> the average value of the overlaid srdf at data polygon(s) weighted
> by area: in case 25% of a given cell in test would overlay polygon r1
> and 75% polygon r2, the value of X1 for that cell would be
> 0.25*1 + 0.75*2

This involves a topology operation on two sets of polygons, and is not 
supported. You can approximate it by creating a finer raster and using 
that subsequently aggregate, for example using the polygons in ?overlay 
to form srdf, then:

test2 <- SpatialGrid(GridTopology(cellcentre.offset=c(178440,329600),
  cellsize=c(50,50), cells.dim=c(80,100)))
o <- overlay(test2, srdf)
o1 <- SpatialPointsDataFrame(coordinates(test2), data=data.frame(o=o))
o2 <- o1[!is.na(o1$o),]
o2$X1 <- srdf$X1[o2$o]
o2$X2 <- srdf$X2[o2$o]
u <- overlay(test, o2)
uu <- aggregate(slot(o2, "data"), list(u), mean)
X1 <- rep(as.numeric(NA), nrow(coordinates(test)))
X2 <- rep(as.numeric(NA), nrow(coordinates(test)))
X1[uu$Group.1] <- uu$X1
X2[uu$Group.1] <- uu$X2
testdf <- 
SpatialGridDataFrame(GridTopology(cellcentre.offset=c(178440,329600),
cellsize=c(500,500), cells.dim=c(8,10)), data=data.frame(X1, X2))
spl <- list("sp.lines", as(srdf, "SpatialLines"), col="grey30")
spplot(testdf, c("X1", "X2"), col.regions=bpy.colors(20), sp.layout=spl)

The finer resolution should be fine enough to capture the relative areas 
of the polygons adequately.

Roger


>
> According to help(overlay),
> "Value
> a numerical array of indices of x on locations of y, or a data.frame
> with (possibly aggregate) properties of x in units of y"
>
> I've tried this:
>> overlay(srdf, test, fn=mean, na.rm)
>
>     X1 X2
> NA   NA NA
> NA.1 NA NA
>
> Cannot go beyond, perhaps overlay() is not suited for this: the
> definition of overlay in its help page
> "overlay combines points (or grids) and polygons by performing
> point-in-polygon operation on all point-polygons combinations"
> is a bit confusing, the first part of the sentence gives me some hope
> but the second one seems to exclude what I'm trying to do.
>
> Any help appreciated (or perhaps pointing to more doc?)
>
>
> Agus
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no