2. Summarizing categorical data

Introduction

This vignette builds upon the sample data for São Miguel (introduced in the previous vignette) to demonstrate the use of exactextractr with categorical land cover data. A sample of the CORINE 2018 raster dataset is included in exactextractr.

As in the previous vignette, the following packages are used:

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

Loading the sample data

First, we load the CORINE land cover data and the concelho boundaries.

clc <- raster(system.file('sao_miguel/clc2018_v2020_20u1.tif',
                     package = 'exactextractr'))
concelhos <- 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 a `dplyr join.

factorValues(clc, c(2, 18, 24))
#>                      landcov
#> 1 Discontinuous urban fabric
#> 2                   Pastures
#> 3          Coniferous forest

Summarizing land cover classifications

One of the most basic questions we might ask is which land cover classification is predominant in each concelho. We can do this with the built-in 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 <- exact_extract(clc, concelhos, 'mode', 
                              append_cols = 'name', progress = FALSE) %>%
  inner_join(clc_classes, by=c(mode = 'value'))

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

head(landcov_fracs)
#> # A tibble: 6 × 3
#> # Groups:   name [1]
#>   name  value    freq
#>   <chr> <dbl>   <dbl>
#> 1 Lagoa     2 0.0630 
#> 2 Lagoa     3 0.0136 
#> 3 Lagoa     7 0.00659
#> 4 Lagoa     9 0.00637
#> 5 Lagoa    12 0.170  
#> 6 Lagoa    18 0.154

landcov_fracs %>%
  inner_join(clc_classes, by = 'value') %>%
  group_by(name) %>%
  arrange(desc(freq)) %>%
  slice_head(n = 3) %>%
  mutate(freq = sprintf('%0.1f%%', 100*freq)) %>%
  knitr::kable()

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

Summarizing population land cover

One extension of the analysis above is to see which land covers are associated with human population in a given concelho. Is the population primary urban or rural?

As described in the previous vignette, the population density raster provides the most robust results in the presence of partially-covered pixels.

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

We are able to perform this analysis because the CORINE sample distributed with exactextractr has been reprojected from its native Lambert Equal Area projection into geographic coordinates consistent with GPW. Otherwise, working with multiple rasters in different projections requires transformation to a common grid using tools such as raster::projectRaster or the gdalwarp command-line utility.

Very little about the call to exact_extract requires changing to incorporate population. 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 <- exact_extract(clc, concelhos, function(df) {
  df %>%
    group_by(name, value) %>%
    summarize(pop = sum(coverage_area * weight / 1e6)) %>%
    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.