Search code examples
rterra

Doing zonal statistics using exact_extract package


library(exactextractr)
library(sf)
library(terra)

Create sample data:

# sample polygon
p <- system.file("ex/lux.shp", package="terra")
p <- vect(p)
p <- p[1, ]
p_sf <- st_as_sf(p)

# sample raster
r <- system.file("ex/elev.tif", package="terra")
r <- rast(r)
r <- crop(r, p , snap = "out", mask = T)

To calculate the mean value of raster for the polygon, I do this:

exact_extract(r, p_sf, 'mean')
# 467.3792

Based on the documentation, this is a area weighted mean where the cell intersection area of the raster is used as a weight to do the averaging.

However, if I want to use another raster layer as weight instead of using the default cell area weight, I did this:

# Create a raster weight 
r_wt <- r / global(r, "sum", na.rm = T)$sum

exact_extract(r, p_sf, weights = r_wt, 'weighted_mean')
# 469.8506

Why are the values different in both cases? I suspect even if I provide an external weight, the external weight is still being multiplied by the cell area weight and the final weight is being used to calculate the weighted mean. Is there any way I can switch off the multiplying with the cell area weight and only use the external weight to do the weighted_mean?


Solution

  • You can write a custom summary function that ignores the fraction of each cell that is covered by the polygon, e.g.:

    exact_extract(r, p_sf, weights = r_wt, fun = function(values, cov, weights) {
      weighted.mean(values, weights)
    })
    

    Some more information is available at https://github.com/isciences/exactextractr#weighted-summary-functions