Hands-on Exercise 7A

Author

Janet Toa

Published

February 26, 2024

Modified

February 26, 2024

21 Choropleth Mapping

21.1 Overview and Learning Outcomes

This hands-on exercise is based on Chapter 21 of the R for Visual Analytics book.

Choropleth mapping involves the symbolisation 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 based on the Master Plan 2014 subzone boundaries.

The learning outcome is to plot functional and truthful choropleth maps using the tmap package.

21.2 Getting Started

21.2.1 Installing and Loading Required Libraries

In this hands-on exercise, the following R packages are used:

tidyverse (i.e. readr, tidyr, dplyr) for performing data science tasks such as importing, tidying, and wrangling data;
sf for handling geospatial data; and
tmap for thematic mapping.

The code chunk below uses the p_load() function in the pacman package to check if the packages are installed. If yes, they are then loaded into the R environment. If no, they are installed, then loaded into the R environment.

pacman::p_load(tidyverse, sf, tmap)

21.2.2 Importing and Preparing Aspatial Data

The aspatial dataset for this hands-on exercise is imported into the R environment using the read_csv() function in the readr package and stored as the R object, popdata.

The data contains information from the Department of Statistics regarding Singapore Residents by Planning Area / Subzone, Age Group, Sex and Type of Dwelling, June 2011-2020 in a csv format (i.e. respopagesextod2011to2020.csv). This is an aspatial data fie. Although it does not contain any coordinates values, but its PA and SZ fields can be used as unique identifiers to geocode to the geospatial shapefile.

popdata = read_csv("data/aspatial/respopagesextod2011to2020.csv")

The tibble data frame, popdata, has 7 columns and 984,656 rows.

Using popdata, a data table consisting of 2020 values would need to be prepared to include the variables “PA”, “SZ”, “YOUNG”, “ECONOMY ACTIVE”, “AGED”, and “TOTAL, DEPENDENCY”.

“YOUNG”: age group 0 to 4 until age groyup 20 to 24;
“ECONOMY ACTIVE”: age group 25-29 until age group 60-64;
“AGED”: age group 65 and above;
“TOTAL”: all age group; and
“DEPENDENCY”: the ratio between young and aged against economy active group.

The pivot_wider() function in the tidyr package, and the mutate(), filter(), group_by(), and select() functions in the dplyr package are used.

popdata2020 = popdata %>%
  filter(Time == 2020) %>%
  group_by(PA, SZ, AG) %>%
  summarise(`POP` = sum(`Pop`)) %>%
  ungroup() %>%
  pivot_wider(names_from=AG, 
              values_from=POP) %>%
  mutate(YOUNG = rowSums(.[3:6])
         +rowSums(.[12])) %>%
mutate(`ECONOMY ACTIVE` = rowSums(.[7:11])+
rowSums(.[13:15]))%>%
mutate(`AGED`=rowSums(.[16:21])) %>%
mutate(`TOTAL`=rowSums(.[3:21])) %>%  
mutate(`DEPENDENCY` = (`YOUNG` + `AGED`)
/`ECONOMY ACTIVE`) %>%
  select(`PA`, `SZ`, `YOUNG`, 
       `ECONOMY ACTIVE`, `AGED`, 
       `TOTAL`, `DEPENDENCY`)

The values in the “PA” and “SZ” fields are converted to uppercase to facilitate joining later on. This is because the values of the “PA” and “SZ” fields are in uppercase and lowercase , while the “SUBZONE_N” and “PLN_AREA_N” are in uppercase.

popdata2020 = popdata2020 %>%
  mutate_at(.vars = vars(PA, SZ), 
          .funs = funs(toupper)) %>%
  filter(`ECONOMY ACTIVE` > 0)

21.2.3 Importing and Preparing Geospatial Data

The geospatial dataset for this hands-on exercise is imported into the R environment using the st_read() function in the sf package and stored as the R object, mpsz.

The data contains information regarding the Master Plan 2014 Subzone Boundary (i.e. MP14_SUBZONE_WEB_PL) in the ESRI shapefile format. This is a geospatial data. It consists of the geographical boundary of Singapore at the planning subzone level. The data is based on URA Master Plan 2014.

mpsz = st_read(dsn = "data/geospatial", 
                layer = "MP14_SUBZONE_WEB_PL")

Reading layer `MP14_SUBZONE_WEB_PL' from data source 
  `C:\jmphosis\ISSS608\Hands-on_Ex\Hands-on_Ex07\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21

The geospatial objects are multipolygon features. There are a total of 323 features and 15 fields in mpsz simple feature data frame. mpsz is in svy21 projected coordinate system. The bounding box provides the x extend and y extend of the data.

mpsz

Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21
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...

21.2.4 Joining Aspatial and Geospatial Data

The left_join() function in the dplyr package is used to join the geographical data and attribute (aspatial) table using planning subzone names, i.e., “SUBZONE_N” and “SZ” as the common identifier.

The function is used with the mpsz simple feature data frame as the left data table is to ensure that the output will be a simple feature data frame.

mpsz_pop2020 = left_join(mpsz, popdata2020,
                          by = c("SUBZONE_N" = "SZ"))

21.3 Plotting Choropleth Maps Using tmap

Two approaches can be used to prepare the thematic map using the tmap package:

Plott a thematic map quickly using the qtm() function; and
Plot a highly customisable thematic map using tmap elements.

21.3.1 Plotting Choropleth Map: qtm()

The easiest and quickest way to draw a choropleth map using the tmap package is by using the qtm() function. It is concise and provides a good default visualisation in many cases.

A cartographic standard choropleth map is plotted below.

The tmap_mode() function with the value “plot” is used to produce a static map. For an interactive mode, the “view” value should be used.
The “fille”argument is used to map the attribute (i.e. DEPENDENCY).

tmap_mode("plot")
qtm(mpsz_pop2020, 
    fill = "DEPENDENCY")

21.3.2 Plotting Choropleth Map: tmap Package

Despite its usefulness in drawing a choropleth map quickly and easily, the disadvantge of the qtm() function is that the aesthetics of individual layers are harder to control. To draw a high quality cartographic choropleth map, tmap’s drawing elements should be used.

tm_shape(mpsz_pop2020)+
  tm_fill("DEPENDENCY", 
          style = "quantile", 
          palette = "Blues",
          title = "Dependency ratio") +
  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() +
  tm_grid(alpha =0.2) +
  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"))

21.3.2.1 Drawing Base Map

The basic building block of tmap is the tm_shape() function, followed by one or more layer elements such as tm_fill() and tm_polygons() functions.

The tm_shape() function is used to define the input data (i.e mpsz_pop2020) and the tm_polygons() function is used to draw the planning subzone polygons

tm_shape(mpsz_pop2020) +
  tm_polygons()

21.3.2.2 Drawing Choropleth Map: 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 the tm_polygons() function.

tm_shape(mpsz_pop2020)+
  tm_polygons("DEPENDENCY")

Note:

The default interval binning used is called “pretty”.

The default colour scheme used is YlOrRd of ColorBrewer.

By default, missing values will be shaded in grey.

21.3.2.3 Drawing Choropleth Map: tm_fill() and tm_border()

The tm_polygons() function is a wraper of the tm_fill() and tm_border() functions. The tm_fill() function shades the polygons by using the default colour scheme and the tm_borders() function adds the borders of the shapefile onto the choropleth map.

tm_shape(mpsz_pop2020)+
  tm_fill("DEPENDENCY")

The “alpha” argument is used in the tm_borders() function to define transparency number between 0 (totally transparent) and 1 (not transparent). By default, the alpha value of the col is used (normally 1).

Besides the “alpha” argument, there are three other arguments for the tm_borders() function:

“col” for border colour;
“lwd” for border line width (default is 1); and
“lty” for border line type (default is “solid”).

tm_shape(mpsz_pop2020)+
  tm_fill("DEPENDENCY") +
  tm_borders(lwd = 0.1,  alpha = 1)

21.3.3 Data Classification Methods of tmap

Most choropleth maps employ some methods of data classification. The aim is to take a large number of observations and group them into data ranges or classes.

The tmap package provides a total ten data classification methods: “fixed”, “sd”, “equal”, “pretty (default)”, “quantile”, “kmeans”, “hclust”, “bclust”, “fisher”, and “jenks”. The “style” argument in the tm_fill() or tm_polygons() functions is used to indicate the data classification method.

21.3.3.1 Built-in Classification Methods

A “jenks” data classification that used 5 classes is plotted.

tm_shape(mpsz_pop2020)+
  tm_fill("DEPENDENCY",
          n = 5,
          style = "jenks") +
  tm_borders(alpha = 0.5)

An “equal” data classification that used 5 classes is plotted.

tm_shape(mpsz_pop2020)+
  tm_fill("DEPENDENCY",
          n = 5,
          style = "equal") +
  tm_borders(alpha = 0.5)

Note: The distribution by the quantile data classification method is more even than the equal data classification method. Warning: Maps Lie!

21.3.3.2 [DIY!] Choropleth Maps Using Different Classification Methods

Using different styles:

tm_shape(mpsz_pop2020)+
  tm_fill("DEPENDENCY",
          n = 5,
          style = "quantile") +
  tm_borders(alpha = 0.5)

tm_shape(mpsz_pop2020)+
  tm_fill("DEPENDENCY",
          n = 5,
          style = "sd") +
  tm_borders(alpha = 0.5)

tm_shape(mpsz_pop2020)+
  tm_fill("DEPENDENCY",
          n = 5,
          style = "pretty") +
  tm_borders(alpha = 0.5)

tm_shape(mpsz_pop2020)+
  tm_fill("DEPENDENCY",
          n = 5,
          style = "kmeans") +
  tm_borders(alpha = 0.5)

tm_shape(mpsz_pop2020)+
  tm_fill("DEPENDENCY",
          n = 5,
          style = "hclust") +
  tm_borders(alpha = 0.5)

tm_shape(mpsz_pop2020)+
  tm_fill("DEPENDENCY",
          n = 5,
          style = "bclust") +
  tm_borders(alpha = 0.5)

Committee Member: 1(1) 2(1) 3(1) 4(1) 5(1) 6(1) 7(1) 8(1) 9(1) 10(1)
Computing Hierarchical Clustering

tm_shape(mpsz_pop2020)+
  tm_fill("DEPENDENCY",
          n = 5,
          style = "log10_pretty") +
  tm_borders(alpha = 0.5)

tm_shape(mpsz_pop2020)+
  tm_fill("DEPENDENCY",
          n = 5,
          style = "fisher") +
  tm_borders(alpha = 0.5)

tm_shape(mpsz_pop2020)+
  tm_fill("DEPENDENCY",
          n = 5,
          style = "dpih") +
  tm_borders(alpha = 0.5)

tm_shape(mpsz_pop2020)+
  tm_fill("DEPENDENCY",
          n = 5,
          style = "headtails") +
  tm_borders(alpha = 0.5)

Using different number of bins:

tm_shape(mpsz_pop2020)+
  tm_fill("DEPENDENCY",
          n = 2,
          style = "jenks") +
  tm_borders(alpha = 0.5)

tm_shape(mpsz_pop2020)+
  tm_fill("DEPENDENCY",
          n = 5,
          style = "jenks") +
  tm_borders(alpha = 0.5)

tm_shape(mpsz_pop2020)+
  tm_fill("DEPENDENCY",
          n = 6,
          style = "jenks") +
  tm_borders(alpha = 0.5)

tm_shape(mpsz_pop2020)+
  tm_fill("DEPENDENCY",
          n = 10,
          style = "jenks") +
  tm_borders(alpha = 0.5)

tm_shape(mpsz_pop2020)+
  tm_fill("DEPENDENCY",
          n = 20,
          style = "jenks") +
  tm_borders(alpha = 0.5)

21.3.3.3 Custom Break

For all the built-in styles, the category breaks are computed internally. In order to override these defaults, the breakpoints can be set explicitly by means of the “breaks“ argument in the tm_fill()function. It is important to note that, in tmap, the breaks include a minimum and maximum. As a result, in order to end up with n categories, n+1 elements must be specified in the breaks option (the values must be in increasing order).

summary(mpsz_pop2020$DEPENDENCY)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
 0.1111  0.7147  0.7866  0.8585  0.8763 19.0000      92

With reference to the summary results above, the break points are set at 0.60, 0.70, 0.80, and 0.90, with 0 set as minimum, and 1.00 set as maximum. The choropleth map is then plotted using the custom breaks.

tm_shape(mpsz_pop2020)+
  tm_fill("DEPENDENCY",
          breaks = c(0, 0.60, 0.70, 0.80, 0.90, 1.00)) +
  tm_borders(alpha = 0.5)

21.3.4 Colour Scheme

The tmap package supports colour ramps either defined by the user or a set of predefined colour ramps from the RColorBrewer package.

21.3.4.1 ColourBrewer Palette

The preferred colour is assigned to the “palette“ argument of the tm_fill() function.

tm_shape(mpsz_pop2020)+
  tm_fill("DEPENDENCY",
          n = 6,
          style = "quantile",
          palette = "Greens") +
  tm_borders(alpha = 0.5)

To reverse the colour shading, a “-” prefix is added.

tm_shape(mpsz_pop2020)+
  tm_fill("DEPENDENCY",
          style = "quantile",
          palette = "-Greens") +
  tm_borders(alpha = 0.5)

21.3.5 Map Layouts

The map layout refers to the combination of all map elements into a cohesive map. Map elements include, among others, the objects to be mapped, title, scale bar, compass, margins and aspects ratios. Colour settings and data classification methods covered in the previous section relate to the palette and break-points are used to affect how the map looks.

21.3.5.1 Legend

In the tmap package, several “legend” argument options are provided to change the placement, format and appearance of the legend.

tm_shape(mpsz_pop2020)+
  tm_fill("DEPENDENCY", 
          style = "jenks", 
          palette = "Blues", 
          legend.hist = TRUE, 
          legend.is.portrait = TRUE,
          legend.hist.z = 0.1) +
  tm_layout(main.title = "Distribution of Dependency Ratio by planning subzone \n(Jenks classification)",
            main.title.position = "center",
            main.title.size = 1,
            legend.height = 0.45, 
            legend.width = 0.35,
            legend.outside = FALSE,
            legend.position = c("right", "bottom"),
            frame = FALSE) +
  tm_borders(alpha = 0.5)

21.3.5.2 Style

The tmap package allows a wide variety of layout settings to be changed using the tmap_style() function.

The “classic” style is used below.

tm_shape(mpsz_pop2020)+
  tm_fill("DEPENDENCY", 
          style = "quantile", 
          palette = "-Greens") +
  tm_borders(alpha = 0.5) +
  tmap_style("classic")

21.3.5.3 Cartographic Furniture

Besides map style, the tmap package also provides arguments to draw other map furniture such as compass, scale bar and grid lines.

The tm_compass(), tm_scale_bar(), and tm_grid() functions are used to add compass, scale bar and grid lines respectively.

tm_shape(mpsz_pop2020)+
  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(lwd = 0.1, alpha = 0.2) +
  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"))

To reset the default style, the tmap_style” function is set back to the “white” value.

tmap_style("white")

21.3.6 Drawing Small Multiple Choropleth Maps

Small multiple maps, also referred to as 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 the tmap package, 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 the tm_facets() function; and
By creating multiple stand-alone maps with the tmap_arrange() function.

21.3.6.1 By Assigning Multiple Values to Aesthetic Argument

Small multiple choropleth maps are created below by defining “ncols” argument in the tm_fill() function.

tm_shape(mpsz_pop2020)+
  tm_fill(c("YOUNG", "AGED"),
          style = "equal", 
          palette = "Blues") +
  tm_layout(legend.position = c("right", "bottom")) +
  tm_borders(alpha = 0.5) +
  tmap_style("white")

Small multiple choropleth maps are created below by assigning multiple values to at least one of the aesthetic arguments.

tm_shape(mpsz_pop2020)+ 
  tm_polygons(c("DEPENDENCY","AGED"),
          style = c("equal", "quantile"), 
          palette = list("Blues","Greens")) +
  tm_layout(legend.position = c("right", "bottom"))

21.3.6.2 By Defining A Group-by Variable in tm_facets()

Small multiple choropleth maps are created using the tm_facets() function to show the different regions.

tm_shape(mpsz_pop2020) +
  tm_fill("DEPENDENCY",
          style = "quantile",
          palette = "Blues",
          thres.poly = 0) + 
  tm_facets(by="REGION_N", 
            free.coords=TRUE, 
            drop.shapes=FALSE) +
  tm_layout(legend.show = FALSE,
            title.position = c("center", "center"), 
            title.size = 20) +
  tm_borders(alpha = 0.5)

21.3.6.3 By Creating Multiple Stand-alone Maps with tmap_arrange()

Small multiple choropleth maps are created as multiple stand-alone maps with the tmap_arrange() function.

youngmap = tm_shape(mpsz_pop2020)+ 
  tm_polygons("YOUNG", 
              style = "quantile", 
              palette = "Blues")

agedmap = tm_shape(mpsz_pop2020)+ 
  tm_polygons("AGED", 
              style = "quantile", 
              palette = "Blues")

tmap_arrange(youngmap, agedmap, asp=1, ncol=2)

21.3.7 Mapping Spatial Object Meeting Selection Criterion

Instead of creating small multiple choropleth maps, the selection function to map spatial objects meeting the selection criterion can also be used. The map below is for the Central Region only.

tm_shape(mpsz_pop2020[mpsz_pop2020$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)