In JaFro96/dasymetric: Create Dasymetric Maps

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 7,
  fig.align="center"
)

library(areal)
library(dasymetric)
library(sf)
library(tmap)
library(dplyr)

Dasymetric mapping is a mapping technique to redistribute spatially extensive data using ancillary information to better represent an underlying statistical surface (Petrov, 2012).

This vignette:

• describes how to prepare land use data and building footprints as ancillary information
• creates dasymetric maps of the population in Münster
• compares the outcome with the true population counts and with a naive area-weighted interpolation

Data preparation

Population data for the districts in Münster can be obtained by their open data platform.

data(population_counts)
map_census = tmap::tm_shape(population_counts) + tmap::tm_polygons("population", title = "Census data (2018)")
map_census

The European Union's Earth Observation Programme - Copernicus - offers freely and openly accessible services that are based on in situ and on satellite data. One of the six thematic services is land monitoring, which provides information on land cover in 44 classes. The latest product is based on Sentinel-2 satellite data from 2018 (European Environment Agency, 2018). The German Bundesamt für Kartographie und Geodäsie provides data with the same nomenclature based on Landbedeckungsmodell Deutschland 2018 (LBM-DE2018) product with a geometric detail of 5 hectare minimum mapping unit.

data(corine_18)
# group polygons by their class ids
grouped_cor = corine_18 |> dplyr::group_by(CLC18) |> dplyr::summarize(geometry = sf::st_union(geometry))
class_names = c("Continuous urban fabric","Discontinuous urban fabric",
                "Industrial or commercial units","Road and rail networks and associated land",
                "Port areas","Airports"," Mineral extraction sites","Dump sites",
                "Green urban areas","Sport and leisure facilities",
                "Non-irrigated arable land","Fruit trees and berry plantations",
                "Pastures","Broad-leaved forest","Coniferous forest","Mixed forest",
                "Natural grassland","Transitional woodland-scrub","Inland marshes",
                "Peat bogs","Water courses","Water bodies")
# add class names
grouped_cor$class_names = class_names
# visualize land use
tmap::tmap_mode("view")
tmap::tm_shape(grouped_cor) + tmap::tm_polygons("class_names") +
  tmap::tm_shape(population_counts) + 
  tmap::tm_borders()

For the following analysis we assume that people live in the two urban fabric land cover classes ( Continuous urban fabric , 111) and Discontinuous urban fabric , 112).

urban_fabric = prep_landuse(grouped_cor)
tmap::tmap_mode("plot")
tmap::tm_shape(population_counts) + 
  tmap::tm_borders() + 
  tmap::tm_shape(urban_fabric) + 
  tmap::tm_polygons("class_names") + 
  tmap::tm_layout(frame.lwd = 3, legend.position = c("right", "top"))

An alternative to land use information as ancillary information is to use building footprints. Buildings are filtered to the classes where we assume that people live there like residential buildings, dormitories, mixed-use buildings, retirement homes, farm houses, houseboats.

data("buildings")
tmap::tm_shape(population_counts) + 
  tmap::tm_borders() + 
  tmap::tm_shape(buildings) + 
  tmap::tm_polygons("funktion") + 
  tmap::tm_layout(frame.lwd = 3, legend.position = c("right", "top"))

Create dasymetric map

# create one unified polygon
source_geom = sf::st_union(population_counts)
# Create source object: Münster with its population
source = sf::st_sf(ID = 1, pop_prediction = sum(population_counts["population"]$population), source_geom)
# Dasymetric map with land use
dm_pop = dasymetric_map(population_counts, source, urban_fabric, extensive = "pop_prediction")
# Dasymetric map with building footprints
dm_pop_buildings = dasymetric_map(population_counts, source, buildings, extensive = "pop_prediction")

# do areal-weighted interpolation for comparison
aw_pop = areal::aw_interpolate(population_counts,NR_STATIST,source = source, sid = ID,weight = "sum", extensive = "pop_prediction", output = "sf")
map_aw = tmap::tm_shape(aw_pop) + tmap::tm_polygons("pop_prediction", title = "AW-interpolation")

# visualize output
map_dm_lu = tmap::tm_shape(dm_pop) + tmap::tm_polygons("pop_prediction", title = "DM land use")
map_dm_bu = tmap::tm_shape(dm_pop_buildings) + tmap::tm_polygons("pop_prediction", title = "DM buildings")
tmap_arrange(map_dm_bu, map_dm_lu, map_aw, map_census, ncol = 2)

Compare the outcome

# compute absolute errors of dasymetric mapping (land use)
dm_pop$pred_inaccuracy = abs(dm_pop$population - dm_pop$pop_prediction)
# compute absolute errors of dasymetric mapping (buildings)
dm_pop_buildings$pred_inaccuracy = abs(dm_pop_buildings$population - dm_pop_buildings$pop_prediction)
# create tibble output to see errors for each district
dm_pop|> select(District = NAME_STATI, Prediction_Inaccuracy_Dasymetric_Mapping = pred_inaccuracy)
dm_pop_buildings|> select(District = NAME_STATI, Prediction_Inaccuracy_Dasymetric_Mapping = pred_inaccuracy)

# compute errors of aw-interpolation
aw_pop$pred_inaccuracy = abs(aw_pop$population - aw_pop$pop_prediction)
# errors = st_join(aw_pop|> select(District = NAME_STATI, Prediction_Inaccuracy_AW_Interpolation = pred_inaccuracy),dm_pop|> select(District = NAME_STATI, Prediction_Inaccuracy_Dasymetric_Mapping = pred_inaccuracy), largest = TRUE)
# output absolute errors
cat("Absolute error of aw-interpolation: ",sum(aw_pop$pred_inaccuracy))
cat("Absolute error of dasymetric mapping using land use data: ",sum(dm_pop$pred_inaccuracy))
cat("Absolute error of dasymetric mapping using building footprints: ",sum(dm_pop_buildings$pred_inaccuracy))

err1 = tmap::tm_shape(dm_pop_buildings) + tmap::tm_polygons("pred_inaccuracy", title = "Absolute Errors (DM buildings)")
err2 = tmap::tm_shape(dm_pop) + tmap::tm_polygons("pred_inaccuracy", title = "Absolute Errors (DM land use)")
err3 = tmap::tm_shape(aw_pop) + tmap::tm_polygons("pred_inaccuracy", title = "Absolute Errors (AW)")
true_vals = tmap::tm_shape(population_counts) + tmap::tm_polygons("population", title = "Census data (2018)")
tmap_arrange(err1,err2,err3, true_vals, ncol = 2)

In this case study more accurate predictions of an extensive variable ( population counts ) were computed. Compared to an absolute error of 305181 for areal-weighted interpolation, dasymetric mapping with land use information reduced the absolute error to 103122 and dasymetric mapping based on building footprints further reduced the absolute error to 71146.

JaFro96/dasymetric documentation built on Dec. 18, 2021, 12:28 a.m.