Weighted means with velox

Question

I'd like to calculate population-weighted average measures from a raster over a set of polygons. I'd prefer to use velox for performance reasons, but I can't figure out how to incorporate weights into the polygon averaging. Below is an MWE demonstrating weighted averaging using raster.

library(raster)
library(sf)

rm(list = ls())

## Make matrix
dim <- c(5, 5)
set.seed(0)
data.mat <- matrix(runif(prod(dim), 0, 100), dim[1], dim[2])
extent <- c(0,1,0,1)
res <- 1/dim
vx <- velox(data.mat, extent, res, crs="")
rast <- vx$as.RasterLayer() # Save rast for comparison to raster::extract() and plotting

## Create sf polygon
pol <-
  st_sfc(st_polygon(list(cbind(
    c(.1, .4, .7, .1), c(.1, .8, .1, .1)
  ))))

## Weighted extract using raster
pol_sp <- as(pol, "Spatial")
wts <- raster::extract(rast, pol_sp, weights = T, normalizeWeights = T, cellnumbers = T, df = T)
weighted.mean(wts$layer, wts$weight) # Weighted average
# [1] 60.43645

Answer 1

The following code will return the weighted mean from velox after running the code above. For large polygons, I believe this could be many times faster than raster::extract. This answer inspired by https://github.com/hunzikp/velox/issues/16.

## Weighted extract using velox
vx_get_weights <- function(rast, poly, normalizeWeights = T) {
  rast$cell <- 1:ncell(rast)
  brk_100 <- disaggregate(rast, fact = 10) 
  brk_100_vx <- velox(brk_100) 
  vx_raw_dt <- setDT(brk_100_vx$extract(poly, fun = NULL, df = TRUE))
  setnames(vx_raw_dt, c("poly_id", "x", "cell"))

  weights <- vx_raw_dt[, .(w = .N / 100), by = .(poly_id, cell, x)]
  if (normalizeWeights) {
    weights[, w := w / sum(w), by = poly_id]
  }
  setorder(weights, poly_id, cell)
  weights
}
weights <- vx_get_weights(rast, pol, normalizeWeights = T)
weighted.mean(rast[weights$cell], weights$w)
# [1] 60.43645