Area-weighted average raster values within each SpatialPolygonsDataFrame polygon (R)

I'm working in R. I wish to find the average value of all rasters cntrystack overlaid by each polygon within a SpatialPolygonsDataFrame cntryshp. The raster object is actually a RasterBrick. So I want to:

1. extract all raster cell values, by individual polygons within cntryshp, for EACH layer in the RasterBrick,
2. average values ACROSS raster cells associated w/ each polygon in cntryshp, for each layer in the RasterBrick, but with the average weighted by the proportion of each raster that falls within the relevant polygon,
3. bind those data to the cntryshp polygon IDs that identify each polygon, to make a data.frame with 1 row per polygon, and one column per layer of the RasterBrick.

A picture of the polygon over the top layer of my raster brick: At the moment, I'm successfully doing this, except without the weighting. So my code goes as thus:

myfunc <- function(x) { x.mean <- colMeans(x, na.rm = FALSE, dims = 1) }

r.vals <- extract(cntrystack, cntryshp)

r.mean <- lapply(r.vals, FUN=myfunc)

mean.matrix <- do.call(rbind, r.mean)

countrydata <- data.frame(cntryshp@data\$ID,mean.matrix)

However, I wish instead to perform a weighted average, in that function, the weight being the proportion of the raster that falls within each polygon. So to be absolutely clear, consider the bottom right polygon (call it poly1) in my posted graphic above. It falls over a number of rasters: [2,1], [3,1], [4,1], [4,2], [5,1],[5,2]. The raster that holds the largest chunk of poly1 is [4,1]. Then [3,1], which holds slightly less of poly1. [2,1] and [5,2] both hold quite small chunks of poly1. [4,2] and [5,1] both do hold a TINY bit of poly1, but the areas are close to negligible. So I want to average the values of each of these 6 rasters, but weighting them according to proportion under poly1, such that the value of [4,1] is weighted most heavily, and the values of [4,2] and [5,1] are weighted very, very little.

One way to do this would be to use instead the weighted mean function:

myfunc <- function(x) { colWeightedMeans(x, w=wghts, na.rm=FALSE, dims = 1)  }

But that requires extracting the proportion/area of each raster w/in each polygon ("arealist"), I think ideally as a list, so that arealist[[i]] gives the vector of raster areas for rasters that fall within polygon i, ordered in the same way that the rasters values are ordered in r.vals, after extracting. I'm struggling with this (perhaps a raster layer needs to be transformed into a polygon, and intersections taken?), and worried that it might not even be the best approach.

An addendum: I figured out that a built in R functionality does this task perfectly. I had simply overlooked it. So now my code is as follows (actually shorter than the original code):

polymeans <- extract(cntrystack, cntryshp_rdf, weights=TRUE, fun=mean) countrydata <- data.frame(cntryshp@data\$ID,polymeans)

You can iterate through each polygon, mask the raster to the subset polygon, coerce to a vector and then calculate weighted mean using rgeos::gArea to return the area of the subset polygon.

Create some data and plot (ignore projection assignment error).

library(raster)
library(sp)
library(rgeos)
x <- raster(xmn=-110, xmx=-90, ymn=40, ymx=60, ncols=10, nrows=10)
x[] <- runif(ncell(x)) * 10
x <- rasterToPolygons(x, fun=function(x){x > 9})
proj4string(x) <- "+proj=lcc +lat_1=48 +lat_2=33 +lon_0=-100 +ellps=WGS84"
y <- raster(xmn=-110, xmx=-90, ymn=40, ymx=60, nrow=100, ncol=100)
y[] <- runif(ncell(y))
proj4string(y) <- "+proj=lcc +lat_1=48 +lat_2=33 +lon_0=-100 +ellps=WGS84"
plot(y)

Create a for loop for processing each polygon. The empty results vector is used to accumulate the weighted.mean values. The if else condition checks that there are values associated with the subset and if not assigns a NA.

results <- vector()
for(j in 1:nrow(x)) {
lsub <- x[j, ]
cr <- raster::crop(y, raster::extent(lsub), snap = "out")
fr <- raster::rasterize(lsub, cr)
r <- na.omit(raster::values(raster::mask(x = cr, mask = fr)))
if (length(r) < 1) {
results <- append(results, NA)
} else {
results <- append(results, stats::weighted.mean(r, rep(gArea(lsub),length(r))))
}
}
results

If you want to weight the mean by the area of intersecting raster cells then you can use a similar for loop. The difference here is that we coerce the subset raster to a SpatialPolygonsDataFrame and then use raster::intersect to merge the raster cell polygons with the polygon. We can then use rgeos::gArea to return areas but in this case it is the areas of the raster cell polygons and not the source polygon providing the subset. We have to handle NA values a bit differently.

results <- vector()
for(j in 1:nrow(x)) {
lsub <- x[j, ]
cr <- raster::crop(y, raster::extent(lsub), snap = "out")
r <- as(cr, "SpatialPolygonsDataFrame")
names(r@data) <- "raster.value"
r <- intersect(lsub, r)
na.idx <- which(is.na(r\$raster.value))
if(length(na.idx) > 0) { r <- r[-na.idx,] }
if (nrow(r) < 1) {
results <- append(results, NA)
} else {
results <- append(results, weighted.mean(r@data[,"raster.value"], gArea(r, byid=TRUE)))
}
}
results

For weighting by cell proportions, you can also use the raster extract function directly with weights = TRUE and normalizeWeights=TRUE. This is certainly the simple approach however, in benchmark tests it is ~5 times slower that the above for loop. I imagine that this will eventually change as the raster package is consistently under development and this is one of the workhorse functions.

( results <- extract(y, x, weights = TRUE, normalizeWeights=TRUE, fun=mean) )

The results vector is the same length as your SpatialPolygonsDataFrame object and is also ordered so can be joined directly back to the polygons.

x@data <- data.frame(x@data, means=results)
str(x@data)
spplot(x, "means")
• Hi, Thanks so much for the code! However, I think you slightly misunderstood my question, or I wasn't totally clear in stating it. – Leah Bevis Oct 9 '16 at 18:47
• Apologies --- I was in an airport as I wrote that, and just then my "free wifi" ran out. So, I think I wasn't totally clear, but I've edited above and added an addendum to my question, to be more precise. Your code weights each raster w/in each polygon item identically --- by the area of the polygon in question, given by rep(gArea(lsub),length(r)). Instead, for each polygon, I wish to weight each raster value by the proportion of the raster that falls within the polygon. I would love your further thoughts, if you can figure out how to do this! (To me, it's much harder.) – Leah Bevis Oct 10 '16 at 6:09
• I expanded my answer to include weighting the mean by the area(s) of the intersecting raster cells. – Jeffrey Evans Oct 10 '16 at 16:47
• Actually, I've figured out a much simpler way to do this, using a built in R functionality. I can simply skip the myfunc I had above, and use this command: polymeans <- extract(cntrystack, cntryshp, weights=TRUE, fun=mean). Then join as before: countrydata <- data.frame(cntryshp@data\$ID,mean.matrix). I hadn't realized that the extract command allowed for functions, nor that it could output weights by raster area. A good lesson in reading command help better... – Leah Bevis Oct 11 '16 at 17:19
• Interesting observation, I bench marked the two approaches, using these small example datasets, using raster::extract with weights = TRUE took 5.0278 seconds whereas my for loop approach took 0.7203 seconds. This notable difference in speed was very surprising because for loops are notoriously slow. Both approaches produced exactly the same results. – Jeffrey Evans Oct 11 '16 at 17:47