1. Summarizing gridded population data

Introduction

This vignette describes the use of exactextractr to summarize population and elevation data from the Gridded Population of the World and EU-DEM datasets. The exactextractr package includes samples of both of these datasets, cropped to the extent of São Miguel, the largest and most populous island of the Azores archipelago.

This example uses the following packages:

library(dplyr)
library(exactextractr)
library(raster)
library(sf)

Loading the sample data

To begin, we load the population count file from GPW. This raster provides the total population in each pixel for the calendar year 2020. On top of the population grid, we plot boundaries for the six municipalities, or concelhos, into which the island is divided. We can see that the population is concentrated along the coastlines, with smaller communities located in volcanic calderas inland.

pop_count <- raster(system.file('sao_miguel/gpw_v411_2020_count_2020.tif',
                                package = 'exactextractr'))
concelhos <- st_read(system.file('sao_miguel/concelhos.gpkg',
                                 package = 'exactextractr'),
                     quiet = TRUE)
plot(pop_count, axes = FALSE)
plot(st_geometry(concelhos), add = TRUE)

Calculating total population

Because the population count raster has been cropped and contains no land area outside of São Miguel, we can calculate the total population of the island using the cellStats function from the raster package.

cellStats(pop_count, 'sum')
#> [1] 145603

exact_extract(pop_count, concelhos, 'sum', progress = FALSE)
#> [1] 14539.875  4149.851 66866.711  5293.968 31920.496  9093.449

concelho_pop <- exact_extract(pop_count, concelhos, 'sum', 
                              append_cols = 'name', progress = FALSE) %>%
  rename(pop = sum) %>%
  arrange(desc(pop)) %>%
  bind_rows(summarize(., name = 'Total', pop = sum(pop)))

plot.new()
plot.window(xlim = c(-25.75, -25.55), ylim = c(37.70, 37.77))
rect(par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4], col = "#9cc0f9")
plot(pop_count, add = TRUE, axes = FALSE, alpha = 0.8, horizontal = TRUE)
plot(st_geometry(concelhos), add = TRUE)

pop_density <- raster(system.file('sao_miguel/gpw_v411_2020_density_2020.tif',
                                  package = 'exactextractr'))

concelho_pop2 <- exact_extract(pop_density, concelhos, 
                               function(density, frac, area) {
                                 sum(density * frac * area)
                               }, 
                               weights = raster::area(pop_density),
                               append_cols = 'name',
                               progress = FALSE)

exact_extract(pop_density, concelhos, 'weighted_sum', weights = 'area')

Population-weighted statistics

Suppose that we are interested in calculating the average elevation of a residence in each of the six concelhos. Loading the EU-DEM elevation data for the island, we can see that each concelho is at least partly occupied by interior mountains, indicating that the results of a simple mean would be unrepresentative of the primarily coastal population.

elev <- raster(system.file('sao_miguel/eu_dem_v11.tif', package = 'exactextractr'))

plot(elev, axes = FALSE, box = FALSE)
plot(st_geometry(concelhos), add = TRUE)

As in the previous section, we avoid working with the population count raster to avoid losing population along the coastline. We can formulate the population-weighted average elevation as in terms of population density and pixel areas as:

\[ \bar{x}_\mathrm{pop} = \frac{ \Sigma_{i=0}^n {x_ic_id_ia_i}}{\Sigma_{i=0}^n{c_id_ia_i}} \] where \(x_i\) is the population of pixel \(i\), \(c_i\) is the fraction of pixel \(i\) that is covered by a polygon, \(d_i\) is the population density of pixel \(i\), and \(a_i\) is the area of pixel \(i\).

If we are working with projected data, or geographic data over a small area such as São Miguel, we can assume all pixel areas to be equivalent, in which case the \(a_i\) components cancel each other out and we are left with the direct usage of population density as a weighting raster:

exact_extract(elev, concelhos,  'weighted_mean', weights = pop_density)

What if pixel areas do vary across the region of our analysis?

One option is to create a scaled population count raster by multiplying the population density and the pixel area. For pixels that are partly covered by water, this inflates the pixel population such that we obtain the correct population when only the land area is covered by a polygon. This requires that we create and maintain a separate raster data set.

Another option is to create a RasterStack of pop_density and area(pop_density), and then write a summary function to handle the necessary processing. We use the summarize_df = TRUE argument to combine the elevation, population density, pixel area, and pixel coverage fraction into a single data frame that is passed to the summary function.

exact_extract(elev, concelhos,
              function(df) {
                weighted.mean(x = df$value,
                              w = df$coverage_fraction * df$pop_density * df$area,
                              na.rm = TRUE)},
              weights = stack(list(pop_density = pop_density,
                                   area = area(pop_density))),
              summarize_df = TRUE,
              progress = FALSE)

This solution shares the same limitations with the previous example using an R summary function with raster::area(): we must precompute an area raster and store it in memory, and we must load all raster values intersecting a given polygon into memory at a single time.

A better solution is to use the coverage_area argument to exact_extract, which specifies that all calculations use the area of each cell that is covered by the polygon instead of the fraction of each cell that is covered by the polygon.

concelho_mean_elev <- exact_extract(elev, concelhos, c('mean', 'weighted_mean'), 
                                    weights = pop_density,
                                    coverage_area = TRUE, 
                                    append_cols = 'name', progress = FALSE)

Here we also calculate the unweighted mean for comparison. We can see that the population-weighted mean elevation is substantially lower than the mean elevation in all concelhos.