Hands-on Exercise 10: Choropleth Mapping with R
Dr. Kam Tin Seong
Assoc. Professor of Information Systems
2019/10/23 (updated: 2020-03-20)
1.0 Overview
Choropleth mapping involves the symbolization of enumeration units, such as countries, provinces, states, counties or census units, using area patterns or graduated colors. For example, a social scientist may need to use a choropleth map to portray the spatial distribution of aged population of Singapore by Master Plan 2014 Subzone Boundary.
In this hands-on exercise, you will learn how to perform choropleth mapping using an R package called tmap.
2.0 Getting Started
Before we get started, we need to ensure that tmap package of R and other related R packages have been installed and loaded into R.
packages = c('sf', 'tmap', 'tidyverse')
for (p in packages){
if(!require(p, character.only = T)){
install.packages(p)
}
library(p,character.only = T)
}3.0 Geospatial Data Wrangling
3.1 The Data
Two data set will be used to create the choropleth map, they are: - URA Master Plan subzone boundary in shapefile format (i.e. MP14_SUBZONE_WEB_PL) - Singapore Residents by Planning Area/Subzone, Age Group and Sex, June 2000 - 2018 in csv format (i.e. respopagsex2000to2018.csv)
3.2 Importing Geospatial Data into R
The code chunk below use the st_read() function of sf package to import MP14_SUBZONE_WEB_PL shapfile into R as a simple feature data frame called mpsz.
mpsz <- st_read(dsn = "data/geospatial",
layer = "MP14_SUBZONE_WEB_PL")## Reading layer `MP14_SUBZONE_WEB_PL' from data source `D:\IS428_AY2018-20T2\03-In-class Exercise\In-class_Ex10-Mapping Geospatial Data in R\data\geospatial' using driver `ESRI Shapefile'
## Simple feature collection with 323 features and 15 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
## epsg (SRID): NA
## proj4string: +proj=tmerc +lat_0=1.366666666666667 +lon_0=103.8333333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +datum=WGS84 +units=m +no_defsUsing the code below to check the content of mpsz
mpsz## Simple feature collection with 323 features and 15 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
## epsg (SRID): NA
## proj4string: +proj=tmerc +lat_0=1.366666666666667 +lon_0=103.8333333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +datum=WGS84 +units=m +no_defs
## First 10 features:
## OBJECTID SUBZONE_NO SUBZONE_N SUBZONE_C CA_IND PLN_AREA_N
## 1 1 1 MARINA SOUTH MSSZ01 Y MARINA SOUTH
## 2 2 1 PEARL'S HILL OTSZ01 Y OUTRAM
## 3 3 3 BOAT QUAY SRSZ03 Y SINGAPORE RIVER
## 4 4 8 HENDERSON HILL BMSZ08 N BUKIT MERAH
## 5 5 3 REDHILL BMSZ03 N BUKIT MERAH
## 6 6 7 ALEXANDRA HILL BMSZ07 N BUKIT MERAH
## 7 7 9 BUKIT HO SWEE BMSZ09 N BUKIT MERAH
## 8 8 2 CLARKE QUAY SRSZ02 Y SINGAPORE RIVER
## 9 9 13 PASIR PANJANG 1 QTSZ13 N QUEENSTOWN
## 10 10 7 QUEENSWAY QTSZ07 N QUEENSTOWN
## PLN_AREA_C REGION_N REGION_C INC_CRC FMEL_UPD_D X_ADDR
## 1 MS CENTRAL REGION CR 5ED7EB253F99252E 2014-12-05 31595.84
## 2 OT CENTRAL REGION CR 8C7149B9EB32EEFC 2014-12-05 28679.06
## 3 SR CENTRAL REGION CR C35FEFF02B13E0E5 2014-12-05 29654.96
## 4 BM CENTRAL REGION CR 3775D82C5DDBEFBD 2014-12-05 26782.83
## 5 BM CENTRAL REGION CR 85D9ABEF0A40678F 2014-12-05 26201.96
## 6 BM CENTRAL REGION CR 9D286521EF5E3B59 2014-12-05 25358.82
## 7 BM CENTRAL REGION CR 7839A8577144EFE2 2014-12-05 27680.06
## 8 SR CENTRAL REGION CR 48661DC0FBA09F7A 2014-12-05 29253.21
## 9 QT CENTRAL REGION CR 1F721290C421BFAB 2014-12-05 22077.34
## 10 QT CENTRAL REGION CR 3580D2AFFBEE914C 2014-12-05 24168.31
## Y_ADDR SHAPE_Leng SHAPE_Area geometry
## 1 29220.19 5267.381 1630379.3 MULTIPOLYGON (((31495.56 30...
## 2 29782.05 3506.107 559816.2 MULTIPOLYGON (((29092.28 30...
## 3 29974.66 1740.926 160807.5 MULTIPOLYGON (((29932.33 29...
## 4 29933.77 3313.625 595428.9 MULTIPOLYGON (((27131.28 30...
## 5 30005.70 2825.594 387429.4 MULTIPOLYGON (((26451.03 30...
## 6 29991.38 4428.913 1030378.8 MULTIPOLYGON (((25899.7 297...
## 7 30230.86 3275.312 551732.0 MULTIPOLYGON (((27746.95 30...
## 8 30222.86 2208.619 290184.7 MULTIPOLYGON (((29351.26 29...
## 9 29893.78 6571.323 1084792.3 MULTIPOLYGON (((20996.49 30...
## 10 30104.18 3454.239 631644.3 MULTIPOLYGON (((24472.11 29...3.3 Importing Attribute Data into R
The codes below used read_csv() function of readr package to import the Residential Population by Age and Sex 2000 to 2017 csv file into R.
popagsex <- read_csv("data/aspatial/respopagsex2000to2018.csv")- It will be in R data frame.
3.4 Data Preparation
Before a thematic map can be prepared, you need to preform the following data preparation.
- Extracting 2018 records only.
- Extracting Males records only.
- Deriving three new variables, namely: Young, Economic Active and Aged.
3.4.1 Data wrangling
The following data wrangling and transformation functions will be used:
- spread() of tidyr package, and
- mutate(), filter(), and select() of dplyr package
popagsex2018_male <- popagsex %>%
filter(Sex == "Males") %>%
filter(Time == 2018) %>%
spread(AG, Pop) %>%
mutate(YOUNG = `0_to_4`+`5_to_9`+`10_to_14`+
`15_to_19`+`20_to_24`) %>%
mutate(`ECONOMY ACTIVE` = rowSums(.[9:13])+
rowSums(.[15:17]))%>%
mutate(`AGED`=rowSums(.[18:22])) %>%
mutate(`TOTAL`=rowSums(.[5:22])) %>%
mutate(`DEPENDENCY` = (`YOUNG` + `AGED`)
/`ECONOMY ACTIVE`) %>%
mutate_at(.vars = vars(PA, SZ), .funs = funs(toupper)) %>%
select(`PA`, `SZ`, `YOUNG`, `ECONOMY ACTIVE`, `AGED`,
`TOTAL`, `DEPENDENCY`) %>%
filter(`ECONOMY ACTIVE` > 0)## Warning: funs() is soft deprecated as of dplyr 0.8.0
## Please use a list of either functions or lambdas:
##
## # Simple named list:
## list(mean = mean, median = median)
##
## # Auto named with `tibble::lst()`:
## tibble::lst(mean, median)
##
## # Using lambdas
## list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## This warning is displayed once per session.3.4.2 Joining the attribute data and geospatial data
Next, left_join() of dplyr is used to join the geographical data and attribute table using planning subzone name e.g. SUBZONE_N and SZ as the common identifier.
mpsz_agemale2018 <- left_join(mpsz, popagsex2018_male,
by = c("SUBZONE_N" = "SZ"))4.0 Choropleth Mapping Geospatial Data Using tmap
Two approaches can be used to prepare thematic map using tmap, they are:
- Plotting a thematic map quickly by using qtm().
- Plotting highly customisable thematic map by using tmap elements.
4.1 Plotting a choropleth map quickly by using qtm()
The easiest and quickest to draw a choropleth map using tmap is using qtm(). It is concise and provides a good default visualisation in many cases.
The code chunk below will draw and cartographic standard choropleth map as shown below.
qtm(mpsz_agemale2018, fill = "DEPENDENCY")## Warning: The shape mpsz_agemale2018 is invalid. See sf::st_is_valid
4.2 Creating a choropleth map by using tmap’s elements
Despite its usefulness of drawing a choropleth map quickly and easily, the disadvantge of qtm() is that it makes aesthetics of individual layers harder to control. To draw a high quality cartographic choropleth map as shown in the figure below tmap’s drawing elements should be used.
## Warning: The shape mpsz_agemale2018 is invalid. See sf::st_is_valid
4.2.1 Drawing a base map
The basic building block of tmap is tm_shape() followed by one or more layer elemments such as tm_fill() and tm_polygons().
In the code chunk below, tm_shpae() is used to define the input data (i.e mpsz_agmale2018) and tm_polygons() is used draw the planning subzone polygons
tm_shape(mpsz_agemale2018) +
tm_polygons()## Warning: The shape mpsz_agemale2018 is invalid. See sf::st_is_valid
4.2.2 Drawing a choropleth map using tm_polygons()
To draw a choropleth map showing the geographical distribution of a selected variable by planning subzone, we just need to assign the target variable such as Dependency to tm_polygons().
tm_shape(mpsz_agemale2018)+
tm_polygons("DEPENDENCY")## Warning: The shape mpsz_agemale2018 is invalid. See sf::st_is_valid
Things to learn from tm_polygons():
- The default interval binning used to draw the choropleth map is called “pretty”. A detail disucssion of the data classification methods supported by tmap will be discussed in sub-section 4.3.
- The default colour scheme used is “YlOrRd” of ColorBrewer. You will learn more about the color pallete in sub-section 4.4.
- By default, Missing value will be shaded in grey.
4.2.3 Drawing a choropleth map using tm_fill() and *tm_border()**
Actually, tm_polygons() is a wraper of tm_fill() and tm_border(). tm_fill() shades the polygons by using the default colour scheme and tm_borders() adds the borders of the shapefile onto the choropleth map.
THe code chunk below draw a choropleth map by using tm_fill() alone.
tm_shape(mpsz_agemale2018)+
tm_fill("DEPENDENCY")## Warning: The shape mpsz_agemale2018 is invalid. See sf::st_is_valid
Notice that the planning subzones are shared according to the respective dependecy values
To add the boundary of the planning subzones, tm_borders will be used as shown in the code chunk below.
tm_shape(mpsz_agemale2018)+
tm_fill("DEPENDENCY") +
tm_borders(alpha = 0.5)## Warning: The shape mpsz_agemale2018 is invalid. See sf::st_is_valid
Notice that light-gray border lines have been added on the choropleth map.
The alpha argument is used to define transparency number between 0 (totally transparent) and 1 (not transparent). By default, the alpha value of the col is used (normally 1).
Beside alpha argument, there are three other arguments for tm_borders(), they are:
- col = border colour,
- lwd = border line width. The default is 1, and
- lty = border line type. The default is “solid”.
4.3 Data classification methods of tmap
Most choropleth maps employ some method of data classification. The point of classification is to take a large number of observations and group them into data ranges or classes.
tmap provides a total ten data classification methods, namely: fixed, sd, equal, pretty (default), quantile, kmeans, hclust, bclust, fisher, and jenks.
To define a data classification method, the style argument of tm_fill() or tm_polygons() will be used.
The code chunk below shows a quantile data classification with 5 classes are used.
tm_shape(mpsz_agemale2018)+
tm_fill("DEPENDENCY",
n = 5,
style = "quantile") +
tm_borders(alpha = 0.5)## Warning: The shape mpsz_agemale2018 is invalid. See sf::st_is_valid
In the code chunk below, euql data classification method is used.
tm_shape(mpsz_agemale2018)+
tm_fill("DEPENDENCY",
n = 5,
style = "equal") +
tm_borders(alpha = 0.5)## Warning: The shape mpsz_agemale2018 is invalid. See sf::st_is_valid
Notice that the distribution of quantile data classification method are more evenly distributed then equal data classification method.
4.4 Colour Scheme
tmap supports colour ramps either defined by the user or a set of predefined colour ramps from the RColorBrewer package.
4.4.1 Using ColourBrewer palette
To change the colour, we assign the prefered colour to palette argument of tm_fill() as shown in the code chunk below.
tm_shape(mpsz_agemale2018)+
tm_fill("DEPENDENCY",
style = "quantile",
palette = "Greens") +
tm_borders(alpha = 0.5)## Warning: The shape mpsz_agemale2018 is invalid. See sf::st_is_valid
Notice that the choropleth map is shaded in green.
To reverse the colour shading, add a “-” prefix.
tm_shape(mpsz_agemale2018)+
tm_fill("DEPENDENCY",
style = "quantile",
palette = "-Greens") +
tm_borders(alpha = 0.5)## Warning: The shape mpsz_agemale2018 is invalid. See sf::st_is_valid
Notice that the colour scheme has been reversed.
4.5 Map Layouts
Map layout refers to the combination of all map elements into a cohensive map. Map elements include among others the objects to be mapped, the title, the scale bar, the compass, margins and aspects ratios, while the colour settings and data classification methods covered in the previous section relate to the palette and break-points used to affect how the map looks.
4.5.1 Map Legend
In tmap, several legend options are provided to change the placement, format and appearance of the legend.
tm_shape(mpsz_agemale2018)+
tm_fill("DEPENDENCY",
style = "quantile",
palette = "Blues",
legend.hist = TRUE,
legend.is.portrait = TRUE,
legend.hist.z = 0.1) +
tm_layout(legend.height = 0.45,
legend.width = 0.35,
legend.outside = FALSE,
legend.position = c("right", "bottom"),
frame = FALSE) +
tm_borders(alpha = 0.5)## Warning: The shape mpsz_agemale2018 is invalid. See sf::st_is_valid
4.5.2 Map style
tmap allows a wide variety of layout settings to be changed. They can be called by using tmap_style().
The code chunk below shows the classic style is used.
tm_shape(mpsz_agemale2018)+
tm_fill("DEPENDENCY",
style = "quantile",
palette = "-Greens") +
tm_borders(alpha = 0.5) +
tmap_style("classic")## tmap style set to "classic"## other available styles are: "white", "gray", "natural", "cobalt", "col_blind", "albatross", "beaver", "bw", "watercolor"## Warning: The shape mpsz_agemale2018 is invalid. See sf::st_is_valid
4.5.3 Cartographic Furniture
Beside map style, tmap also also provides arguments to draw other map furniture such as compass, scale bar and grid lines.
In the code chunk below, tm_compass(), tm_scale_bar() and tm_grid() are used to add compass, scale bar and grid lines onto the choropleth map.
tm_shape(mpsz_agemale2018)+
tm_fill("DEPENDENCY",
style = "quantile",
palette = "Blues",
title = "No. of persons") +
tm_layout(main.title = "Distribution of Dependency Ratio \nby planning subzone",
main.title.position = "center",
main.title.size = 1.2,
legend.height = 0.45,
legend.width = 0.35,
frame = TRUE) +
tm_borders(alpha = 0.5) +
tm_compass(type="8star", size = 2) +
tm_scale_bar(width = 0.15) +
tm_grid() +
tm_credits("Source: Planning Sub-zone boundary from Urban Redevelopment Authorithy (URA)\n and Population data from Department of Statistics DOS",
position = c("left", "bottom"))## Warning: The shape mpsz_agemale2018 is invalid. See sf::st_is_valid
To reset the default style, the code chunk use the code chunk below.
tmap_style("white")## tmap style set to "white"## other available styles are: "gray", "natural", "cobalt", "col_blind", "albatross", "beaver", "bw", "classic", "watercolor"4.6 Drawing Small Multiple Choropleth Maps
Small multiple maps, also referred to facet maps, are composed of many maps arrange side-by-side, and sometimes stacked vertically. Small multiple maps enable the visualisation of how spatial relationships change with respect to another variable, such as time.
In tmap, small multiple maps can be plotted in three ways:
- by assigning multiple values to at least one of the asthetic arguments,
- by defining a group-by variable in tm_facets(), and
- by creating multiple stand-alone maps with tmap_arrange().
4.6.1 By assigning multiple values to at least one of the aesthetic arguments
In this example, small multiple choropleth maps are created by defining ncols in tm_fill()
tm_shape(mpsz_agemale2018)+
tm_fill(c("YOUNG", "AGED"),
style = "equal",
palette = "Blues") +
tm_layout(legend.position = c("right", "bottom")) +
tm_borders(alpha = 0.5) +
tmap_style("white")## tmap style set to "white"## other available styles are: "gray", "natural", "cobalt", "col_blind", "albatross", "beaver", "bw", "classic", "watercolor"## Warning: The shape mpsz_agemale2018 is invalid. See sf::st_is_valid
In this example, small multiple choropleth maps are created by assigning multiple values to at least one of the aesthetic arguments
tm_shape(mpsz_agemale2018)+
tm_polygons(c("DEPENDENCY","DEPENDENCY"),
style = c("equal", "quantile"),
palette = list("Blues","Blues")) +
tm_layout(legend.position = c("right", "bottom"))## Warning: The shape mpsz_agemale2018 is invalid. See sf::st_is_valid
4.6.2 By defining a group-by variable in tm_facets()
In this example, multiple small choropleth maps are created by using tm_facets().
tm_shape(mpsz_agemale2018) +
tm_fill("DEPENDENCY",
style = "quantile",
palette = "Blues",
thres.poly = 0) +
tm_facets(by="REGION_N",
free.coords=TRUE,
drop.shapes=TRUE) +
tm_layout(legend.show = FALSE,
title.position = c("center", "center"),
title.size = 20) +
tm_borders(alpha = 0.5)
4.6.3 By creating multiple stand-alone maps with tmap_arrange()
In this example, multiple small choropleth maps are created by creating multiple stand-alone maps with tmap_arrange().
youngmap <- tm_shape(mpsz_agemale2018)+
tm_polygons("YOUNG",
style = "quantile",
palette = "Blues")
agedmap <- tm_shape(mpsz_agemale2018)+
tm_polygons("AGED",
style = "quantile",
palette = "Blues")
tmap_arrange(youngmap, agedmap, asp=1, ncol=2)## Warning: The shape mpsz_agemale2018 is invalid. See sf::st_is_valid
## Warning: The shape mpsz_agemale2018 is invalid. See sf::st_is_valid
4.7 Mappping Spatial Object Meeting a Selection Criterion
Instead of creating small multiple choropleth map, you can also use selection funtion to map spatial objects meeting the selection criterion.
tm_shape(mpsz_agemale2018[mpsz_agemale2018$REGION_N=="CENTRAL REGION", ])+
tm_fill("DEPENDENCY",
style = "quantile",
palette = "Blues",
legend.hist = TRUE,
legend.is.portrait = TRUE,
legend.hist.z = 0.1) +
tm_layout(legend.outside = TRUE,
legend.height = 0.45,
legend.width = 5.0,
legend.position = c("right", "bottom"),
frame = FALSE) +
tm_borders(alpha = 0.5)## Warning: The shape mpsz_agemale2018[mpsz_agemale2018$REGION_N == "CENTRAL
## REGION", is invalid. See sf::st_is_valid
4.8 Creating An Interactive Map Application
In this section, you will learn how to create an interactive mapping application by integrating leaflet and tmap
tmap_mode("view")## tmap mode set to interactive viewingmpsz_wgs84 <- st_transform(mpsz_agemale2018, 4326)
tm_shape(mpsz_wgs84)+
tm_fill("DEPENDENCY",
n = 6,
style = "quantile",
palette = "Blues") +
tm_borders(alpha = 0.5)## Warning: The shape mpsz_wgs84 is invalid. See sf::st_is_validTo switch tmap’s viewer back to plotting mode, the code chunk below will be used.
tmap_mode("plot")## tmap mode set to plotting