A dasymetric map of social vulnerability in Belgium
Very often demographic and socio-economic features are cartographically visualized by making use of choropleth mapping techniques. Choropleth maps spatially aggregate data into (administrative) geographical units which implies that study attributes (e.g. population density) are assumed to be distributed homogeneously throughout the geographical unit. Choropleth maps can be misleading since the geographical areas have varying levels of internal homogeneity and therefore don’t take into account variations in the spatial organization and land use. This can lead to interpretation issues, such as the ecological fallacy and the modifiable areal unit problem.
To improve area homogeneity, dasymetric mapping can be used. Dasymetric mapping is a technique in which data is more accurately distributed by using ancillary information to design internally homogeneous zones. The most common type of ancillary data is represented by land cover classifications with zones ranging from uninhabited wilderness to urban development areas. Regions known to be unpopulated (e.g. water bodies) are cropped out of the original choropleth districts and appear empty on the map. Population characteristics are reallocated to appropriate ancillary zones, resulting in a more accurate portrayal of where people live within an administrative boundary.
In this study we use dasymetric mapping to visualize social vulnerability relative to land use. Based on median income and the number of people entitled to a living wage in Belgian municipalities on the one hand and Corine land cover data on the other hand, we create an improved representation of social vulnerability. The attributes we study can be distributed over new and smaller areas to create a more realistic visualization than is the case with choropleth mapping.
Data
The data we use to produce a dasymetric map are the same as decribed in a previous blog post about mapping bivariate relations in R. We downloaded the data of median income in Belgium and the number of persons entitled to a living wage (http://stat.mi-is.be/nl/dashboard/ris_scatter?menu=Scatterplots) and use data from the period of august 2019 to july 2020. The data are aggregated by municipality (N = 581) and report the median income and the number of persons per 1000 inhabitants that receive a minimal subsistence income payment. We create a bivariate map of the relationship between both variables and as decribed above we use land cover data to geographically represent a more realistic visualization of this relationship (several code lines have been taken from the excellent paper written by D. Royé *).
First, we downloaded Corine land cover data for 2018 (CLC 2018) from the Copernicus website (https://land.copernicus.eu/pan-european/corine-land-cover/clc2018). We use a shapefile of the Belgian municipalities and crop the land cover raster data to the extent of the shapefile.
# packages
library(tidyverse)
library(sf)
library(readxl)
library(biscale)
library(patchwork)
library(raster)
library(sysfonts)
library(showtext)
library(raster)
setwd("xxxxx")
# raster of Corine Land Cover 2018
urb <- raster("U2018_CLC2018_V2020_20u1.tif")
# shapefile Belgium
belgie <- read_sf("AD_2_Municipality.shp")
belgie_prj <- st_transform(belgie, projection(urb))
B <- st_zm(belgie_prj)
bel <- crop(urb,B) %>%
mask(B)
plot(bel,bty="n",axes="F",legend=FALSE,box=FALSE,font.main=1)
title(main="Corine landcover types Belgium",col.main="darkgreen",cex.main=0.9,line=0.35)
The Corine Land Cover nomenclature is a 3-level hierarchical classification system and has 44 classes at (the third and) the most detailed level. We follow the assumption that people only live in areas with code 111 and 112 (continuouis urban fabric and discontinuous urban fabric).
bel[!bel %in% 1:2] <- NA
plot(bel,bty="n",axes="F",legend=FALSE,box=FALSE,font.main=1)
title(main="Corine landcover types Belgium after removing unpopulated areas",col.main="darkgreen", cex.main=0.9,line=0.35)
Geographical disaggregation allows the attributes we study to be distributed over new, smaller areas to create a more realistic visualization.
# change projection
bel <- projectRaster(bel, crs=CRS("+proj=longlat +datum=WGS84 +no_defs"))
# convert raster into a sf object
bel <- as.data.frame(bel, xy=TRUE, na.rm = TRUE) %>%
st_as_sf(coords = c("x","y"), crs=4326)
# add the columns of the coordinates
bel <- bel %>% rename(urb = 1) %>%
cbind(st_coordinates(bel))
# read shapefile with data about income and the number of inhabitants entitled to receive a living wage
belgium_medincome <- read_sf("belgium_medincome.shp")
2. Create bivariate classification and dasymetric map of social vulnerability
# create bivariate classification
mapbivar <- bi_class(belgium_medincome, md_nkmn, leflnrs, style = "quantile", dim = 3) %>%
mutate(bi_class = ifelse(str_detect(bi_class, "NA"), NA, bi_class))
# results
head(dplyr::select(mapbivar, md_nkmn, leflnrs, bi_class))
Simple feature collection with 6 features and 3 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: 4.17209 ymin: 50.11018 xmax: 6.034387 ymax: 50.77812
geographic CRS: WGS 84
# A tibble: 6 x 4
md_nkmn leflnrs bi_class geometry
1 22694 26.0 1-3 (((5.44713 50.207, 5.447241 50.20663, 5.447534 50.20611, 5.447693 50.20591, 5.447774.~
2 25913 7.54 2-2 (((4.272817 50.7386, 4.272835 50.73858, 4.27368 50.73767, 4.27387 50.73749, 4.273564.~
3 23539 13.6 1-3 (((4.785421 50.22539, 4.785649 50.22537, 4.785734 50.22472, 4.785852 50.22395, 4.785.~
4 19442 44.1 1-3 (((4.968377 50.21704, 4.968252 50.21642, 4.968216 50.21636, 4.968203 50.21636, 4.968.~
5 23141 8.13 1-2 (((6.021366 50.18216, 6.021107 50.18208, 6.020423 50.18178, 6.020207 50.18172, 6.019.~
6 22271 20.6 1-3 (((5.069918 50.17642, 5.069292 50.17615, 5.068089 50.17559, 5.066567 50.1749, 5.0664.~
# redistribute urban pixels to inequality
mapdasi <- st_join(bel, st_transform(mapbivar, 4326))
# bivariate legend
legend2 <- bi_legend(pal = "DkViolet",
dim = 3,
xlab = "higher no. people\nwith living wage",
ylab = "higher income",
size = 9)
# download font
font_add_google("Montserrat", "Montserrat")
showtext_auto()
# dasymetric map
dmap <- ggplot(mapdasi) +
geom_tile(aes(X, Y,
fill = bi_class),
show.legend = FALSE) +
geom_sf(data = belgie,
color = "grey80",
fill = NA,
size = 0.2) +
annotation_custom(ggplotGrob(legend2),
xmin = 3, xmax = 4,
ymin = 49.5, ymax = 50.25) +
bi_scale_fill(pal = "DkViolet",
dim = 3,
na.value = "grey90") +
labs(title = "Dasymetric map of social vulnerability in Belgium",x = "", y ="") +
bi_theme() +
theme(plot.title = element_text(family = "Montserrat", size = 18, face = "bold")) +
coord_sf(crs = 4326)
dmap
# choropleth map
chmap <- ggplot(mapbivar) +
geom_sf(aes(fill = bi_class),
colour = NA,
size = .1,
show.legend = FALSE) +
geom_sf(data = belgie,
color = "white",
fill = NA,
size = 0.2) +
bi_scale_fill(pal = "DkViolet",
dim = 3,
na.value = "grey90") +
labs(title = "Choropleth map of social vulnerability in Belgium", x = "", y ="") +
bi_theme() +
theme(plot.title = element_text(family = "Montserrat", size = 18, face = "bold")) +
coord_sf(crs = 4326)
# combine choropleth and dasymetric map
p <- chmap | dmap
p
3. Conclusion
We used land cover data to redistribute the bivariate relationship between income and the number of people with a living wage. In zones that are heavily urbanized, the dasymetric mapping method does not show much difference from the choropleth method. However, the mapping of social deprivation relative to land use shows that the dasymetric method results in a more accurate visualization of relative levels of social deprivation, especially in the southern part of Belgium. This shows that ancillary data should be used to interpret levels of habitation.