Zonal Statistics of Categorical Data with exactextractr

Introduction

This vignette is the second in a two-part series demonstrating the use of the extactextractr package. In the first vignette, we demonstrated several considerations for performing zonal extractions on meandering coastal regions such as São Miguel. In this vignette, we illustrate the the use of exactextractr with categorical land cover data.

We will use the following packages for this vignette.

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

Loading the sample data

First, we load the CORINE land cover data and the concelho boundaries used in the first vignette. A sample of the CORINE 2018 raster dataset is included in exactextractr.

clc <- raster::raster(system.file('sao_miguel/clc2018_v2020_20u1.tif',
 package = 'exactextractr'))
concelhos <- sf::st_read(system.file('sao_miguel/concelhos.gpkg',
 package = 'exactextractr'),
 quiet = TRUE)

The land cover class descriptions are provided in a separate DBF file. We read this in to a data frame, then use levels() to associate the class descriptions with the raster.

clc_classes <- foreign::read.dbf(system.file('sao_miguel/clc2018_v2020_20u1.tif.vat.dbf',
 package = 'exactextractr'),
 as.is = TRUE) %>%
 dplyr::select(value = Value,
 landcov = LABEL3)

levels(clc) <- list(data.frame(ID = clc_classes$value,
 landcov = clc_classes$landcov))

This association provides us with a way to look up the description for a given ID. Alternatively, we can relate the values using merge or one of dplyr’s joins.

raster::factorValues(clc, c(2, 18, 24))

Summarizing Land Cover Classifications

One of the first questions we might have is which land cover classification is predominant in each concelho. We can do this with exactextractr’s mode summary operation. The minority and variety operations are also applicable to categorical data and provide the least-common classification and number of distinct classifications, respectively.

landcov_mode <- exactextractr::exact_extract(clc, concelhos, 'mode', 
 append_cols = 'name', progress = FALSE) %>%
 dplyr::inner_join(clc_classes, by=c(mode = 'value'))

landcov_fracs <- exactextractr::exact_extract(clc, concelhos, function(df) {
 df %>%
 dplyr::mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
 dplyr::group_by(name, value) %>%
 dplyr::summarize(freq = sum(frac_total))
}, summarize_df = TRUE, include_cols = 'name', progress = FALSE)

head(landcov_fracs)

fracs_summary <-
 landcov_fracs %>%
 dplyr::inner_join(clc_classes, by = 'value') %>%
 dplyr::group_by(name) %>%
 dplyr::arrange(dplyr::desc(freq)) %>%
 dplyr::slice_head(n = 3) %>%
 dplyr::mutate(freq = sprintf('%0.1f%%', 100 * freq))
 
 kableExtra::kbl(fracs_summary[, -1], 
 col.names = c("Cover", "Frequency", "Cover Description"),
 align = c('c',"c", "l"),
 full_width = F, 
 position = "center") %>%
 kableExtra::pack_rows(index = table(fracs_summary$name))

landcov_areas <- exactextractr::exact_extract(clc, concelhos, function(df) {
 df %>%
 dplyr::group_by(name, value) %>%
 dplyr::summarize(area_km2 = round(sum(coverage_area) / 1e6))
}, summarize_df = TRUE, coverage_area = TRUE, include_cols = 'name', progress = FALSE)

Summarizing Population Land Cover

Taking this a step further, one may be interested in the correlation between population and land cover class in a given concelho; i.e. is the population primarily urban or rural? As described in the previous vignette, using the population density raster is more accurate when pixels are partially covered by the boundary layer.

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

This analysis is only possible because the sample CORINE dataset packaged with exactextractr has been reprojected from its native Lambert Equal Area projection into geographic coordinates used for the Gridded Population of the World; the projections must match. If you are working with raster data in different projections, you must transform them to a common projection using a tool like raster::projectRaster, the gdalwarp command-line utility, or a projection tool bundled with a standalone GIS product like QGIS.

We can use a similar call to exact_extract to incorporate the population data into our analysis. We set weights = pop_density and, since we are using the summarize_df option, a column called weight will be added to the data frame passed to the summary function. We also set coverage_area = TRUE so that we can multiply the density by the covered pixel area to get a population count.

landcov_pop_areas <- exactextractr::exact_extract(clc, concelhos, function(df) {
 df %>%
 dplyr::group_by(name, value) %>%
 dplyr::summarize(pop = sum(coverage_area * weight / 1e6)) %>%
 dplyr::mutate(pop_frac = pop / sum(pop))
}, weights = pop_density, coverage_area = TRUE,
 summarize_df = TRUE, include_cols = 'name')

Looking at the highest-population land cover type in each concelho, we can can see that the western/central concelhos of Lagoa, Ponta Delgada, Ribeira Grande, and Vila Franca do Campo have a more urban population than Nordeste or Povoação to the east.