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

[Robert J. Hijmans]

Or essentially the same with using the polygonValues function in
raster (with weights=TRUE)

library(raster)
r <- raster(test)
r[] <- 1:ncell(r)
# finding out the fraction covered by each polygon
v <- polygonValues(srdf, r, weights=TRUE)

# below could possibly done more elegantly in an sapply
vv <- matrix(nrow=0, ncol=2)
polyvar <- 2  # second variable of the SpPolDF
for (i in 1:length(v)) {
	a <- v[[i]]
	vv <- rbind(vv, cbind(a[,1], a[,2] * srdf at data[i,polyvar]))
}

# sum over cells
vvv <- tapply(vv[,2], vv[,1], sum)
rr <- r
r[] <- NA
# apply cell values to raster
r[as.integer(rownames(vvv))] <- vvv
plot(r)
plot(srdf, add=T)


Robert

On Thu, May 20, 2010 at 12:17 PM, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> 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
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>