PART 1: Mapping Geospatial Point Data with R

Overview

Proportional symbol maps (also known as graduate symbol maps) are a class of maps that use the visual variable of size to represent differences in the magnitude of a discrete, abruptly changing phenomenon, e.g. counts of people. Like choropleth maps, you can create classed or unclassed versions of these maps. The classed ones are known as range-graded or graduated symbols, and the unclassed are called proportional symbols, where the area of the symbols are proportional to the values of the attribute being mapped. In this hands-on exercise, you will learn how to create a proportional symbol map showing the number of wins by Singapore Pools’ outlets using an R package called tmap.

Learning outcome

By the end of this exercise, we will acquire the following skills by using appropriate R packages: - To import an aspatial data file into R. - To convert it into simple point feature data frame and at the same time, to assign an appropriate projection reference to the newly create simple point feature data frame. - To plot interactive proportional symbol maps.

Import the packages required

Ensure that tmap package of R and other related R packages have been installed and loaded onto R

packages = c('sf', 'tmap', 'tidyverse')
for (p in packages){
  if(!require(p, character.only = T)){
    install.packages(p)
  }
  library(p,character.only = T)
}
## Loading required package: sf
## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
## Loading required package: tmap
## Loading required package: tidyverse
## ── Attaching packages ────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.0     ✓ purrr   0.3.3
## ✓ tibble  2.1.3     ✓ dplyr   0.8.5
## ✓ tidyr   1.0.2     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0
## ── Conflicts ───────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

Geospatial Data Wrangling

The data

The data set use for this hands-on exercise is called SGPools_svy21. The data is in csv file format.

Figure below shows the first 15 records of SGPools_svy21.csv. It consists of seven columns. The XCOORD and YCOORD columns are the x-coordinates and y-coordinates of SingPools outlets and branches. They are in Singapore SVY21 Projected Coordinates System.

Data Import and Preparation

The code chunk below uses read_csv() function of readr package to import SGPools_svy21.csv into R as a tibble data frame called sgpools

sgpools <- read_csv("data/aspatial/SGPools_svy21.csv")
## Parsed with column specification:
## cols(
##   NAME = col_character(),
##   ADDRESS = col_character(),
##   POSTCODE = col_double(),
##   XCOORD = col_double(),
##   YCOORD = col_double(),
##   `OUTLET TYPE` = col_character(),
##   `Gp1Gp2 Winnings` = col_double()
## )

After importing the data file into R, it is important for us to examine if the data file has been imported correctly. Use list() to check

list(sgpools) 
## [[1]]
## # A tibble: 306 x 7
##    NAME      ADDRESS       POSTCODE XCOORD YCOORD `OUTLET TYPE` `Gp1Gp2 Winning…
##    <chr>     <chr>            <dbl>  <dbl>  <dbl> <chr>                    <dbl>
##  1 Livewire… 2 Bayfront A…    18972 30842. 29599. Branch                       5
##  2 Livewire… 26 Sentosa G…    98138 26704. 26526. Branch                      11
##  3 SportsBu… Lotus Lounge…   738078 20118. 44888. Branch                       0
##  4 SportsBu… 1 Selegie Rd…   188306 29777. 31382. Branch                      44
##  5 Prime Se… Blk 542B Ser…   552542 32239. 39519. Branch                       0
##  6 Singapor… 1A Woodlands…   731001 21012. 46987. Branch                       3
##  7 Singapor… Blk 64 Circu…   370064 33990. 34356. Branch                      17
##  8 Singapor… Blk 88 Circu…   370088 33847. 33976. Branch                      16
##  9 Singapor… Blk 308 Anch…   540308 33910. 41275. Branch                      21
## 10 Singapor… Blk 202 Ang …   560202 29246. 38943. Branch                      25
## # … with 296 more rows

Notice that the sgpools data in tibble data frame and not the common R data frame

Creating a sf data frame from an aspatial data frame

The code chunk below converts sgpools data frame into a simple feature data frame by using st_as_sf() of sf packages

sgpools_sf <- st_as_sf(sgpools, coords = c("XCOORD", "YCOORD"), crs= 3414)

Things to learn from the arguments above:

The coords argument requires you to provide the column name of the x-coordinates first then followed by the column name of the y-coordinates. The crs argument required you to provide the coordinates system in epsg format. EPSG: 3414 is Singapore SVY21 Projected Coordinate System. You can search for other country’s epsg code by refering to epsg.io. Geometry is now added into the data frame -> do View(sgpools_sf)

list(sgpools_sf)
## [[1]]
## Simple feature collection with 306 features and 5 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 7844.194 ymin: 26525.7 xmax: 45176.57 ymax: 47987.13
## CRS:            EPSG:3414
## # A tibble: 306 x 6
##    NAME  ADDRESS POSTCODE `OUTLET TYPE` `Gp1Gp2 Winning…
##    <chr> <chr>      <dbl> <chr>                    <dbl>
##  1 Live… 2 Bayf…    18972 Branch                       5
##  2 Live… 26 Sen…    98138 Branch                      11
##  3 Spor… Lotus …   738078 Branch                       0
##  4 Spor… 1 Sele…   188306 Branch                      44
##  5 Prim… Blk 54…   552542 Branch                       0
##  6 Sing… 1A Woo…   731001 Branch                       3
##  7 Sing… Blk 64…   370064 Branch                      17
##  8 Sing… Blk 88…   370088 Branch                      16
##  9 Sing… Blk 30…   540308 Branch                      21
## 10 Sing… Blk 20…   560202 Branch                      25
## # … with 296 more rows, and 1 more variable: geometry <POINT [m]>

The output shows that sgppols_sf is in point feature class. It’s epsg ID is 3414. The bbox provides information of the extend of the geospatial data.

Drawing proportional symbol map

To create an interactive proportional symbol map in R, the view mode of tmap will be used

tmap_mode("view")
## tmap mode set to interactive viewing

To create an interactive point symbol map

tm_shape(sgpools_sf)+
tm_bubbles(col = "pink",
           size = 1,
           border.col = "black",
           border.lwd = 1)

Make it proportional, bubbles of different sizes

Assign a numerical variable to the size visual attribute

tm_shape(sgpools_sf)+
tm_bubbles(col = "pink",
           size = "Gp1Gp2 Winnings",
           border.col = "black",
           border.lwd = 0.1)
## Legend for symbol sizes not available in view mode.

Different colour, one colour for each

tm_shape(sgpools_sf)+
tm_bubbles(col = "OUTLET TYPE", 
          size = "Gp1Gp2 Winnings",
          alpha=0.5,
          border.col = "black",
          border.lwd = 1)
## Legend for symbol sizes not available in view mode.

Twin brothers

Produce multiple maps with synchronised zoom and pan settings

tm_shape(sgpools_sf) +
  tm_bubbles(col = "OUTLET TYPE", 
          size = "Gp1Gp2 Winnings",
          border.col = "black",
          border.lwd = 1) +
  tm_facets(by= "OUTLET TYPE",
            nrow = 1,
            sync = TRUE)
## Legend for symbol sizes not available in view mode.
tmap_mode("plot")
## tmap mode set to plotting

PART 2: Choropleth Mapping

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 exercise, we will learn how to perform choropleth mapping using tmap package in R.

Getting Started

Beside tmap package, four other R packages will be used, they are: readr, tidyr, dplyr and sf. Among the four packages, readr, tidyr and dplyr are part of tidyverse package. Notice that, we only need to install tidyverse instead of readr, tidyr and dplyr individually.

packages = c('sf', 'tmap', 'tidyverse')
for (p in packages){
  if(!require(p, character.only = T)){
    install.packages(p)
  }
  library(p,character.only = T)
}

Importing geospatial data into R

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). 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.

Singapore Residents by Planning Area/Subzone, Age Group and Sex, June 2000 - 2018 in csv format (i.e. respopagsex2000to2018.csv). This is an aspatial data fie. Although it does not contain any coordinates values, but it’s PA and SZ fields can be used as unique identifiers to georeference to MP14_SUBZONE_WEB_PL shapefile.

mpsz <- st_read(dsn = "data/geospatial", 
                layer = "MP14_SUBZONE_WEB_PL")
## Reading layer `MP14_SUBZONE_WEB_PL' from data source `/Users/theodora/Desktop/IS415/Lesson 3/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
## 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

Examine 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
## 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...

Importing attribute data into R

Next, we will import respopagsex2000to2018.csv file into RStudio and save the file into an R dataframe called popagsex.

popagsex <- read_csv("data/aspatial/respopagsex2000to2018.csv")
## Parsed with column specification:
## cols(
##   PA = col_character(),
##   SZ = col_character(),
##   AG = col_character(),
##   Sex = col_character(),
##   Pop = col_double(),
##   Time = col_double()
## )

Data Preparation

Before a thematic map can be prepared, you need to preform the following data preparation. 1. Extracting 2018 records only. 2. Extracting Males records only. 3. Deriving three new variables, namely: Young, Economic Active and Aged.

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") %>% #interested in Male
  filter(Time == 2018) %>% #interested in 2018
  spread(AG, Pop) %>% #like a pivot
  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`) %>% ##calculate dependency ratio
mutate_at(.vars = vars(PA, SZ), .funs = funs(toupper)) %>% ##change all the words in the column PA and SZ to upper case
select(`PA`, `SZ`, `YOUNG`, `ECONOMY ACTIVE`, `AGED`, 
       `TOTAL`, `DEPENDENCY`) %>%
  filter(`ECONOMY ACTIVE` > 0)

Joining attribute data and geospatial data

left_join() of dplyr is used to join the geographical data and attribute table using planning subzone name eg SUBZONE_N and SZ as the common identifier

mpsz_agemale2018 <- left_join(mpsz, popagsex2018_male, 
                              by = c("SUBZONE_N" = "SZ"))
## Warning: Column `SUBZONE_N`/`SZ` joining factor and character vector, coercing
## into character vector

Choropleth Mapping Geospatial Data using tmap

Two approaches can be used to prepare thematic map using tmap: 1. Plotting a thematic map quickly by using qtm() 2. Plotting highly customisable thematic map by using tmap elements

Plot choropleth map using qtm()

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

qtm(mpsz_agemale2018, fill = "DEPENDENCY")
## Warning: The shape mpsz_agemale2018 is invalid. See sf::st_is_valid

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.

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

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.

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(lwd = 0.1,  alpha = 1)
## 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”.

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.

Plotting choropleth maps with built-in classification methods

tm_shape(mpsz_agemale2018)+
  tm_fill("DEPENDENCY",
          n = 8,
          style = "jenks") +
  tm_borders(alpha = 0.5)
## Warning: The shape mpsz_agemale2018 is invalid. See sf::st_is_valid

Equal data classification method

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.

Plot choropleth map

For all the built-in styles, the category breaks are computed internally. In order to overide these defaults, the breakpoints can be set explicitly by means of the breaks argument to the tm_fill(). 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).

Before we get started, it is always a good practice to get some descriptive statistics on the variable before setting the break points. Code chunk below will be used to compute and display the descriptive statistics of DEPENDENCY field.

summary(mpsz_agemale2018$DEPENDENCY)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##  0.0000  0.6480  0.6847  0.6708  0.7393  0.9559      90

With reference to the results above, we set break point at 0.60, 0.70, 0.80, and 0.90. In additional, we also need to include a minimum and maximum, which we set at 0 and 100. Our breaks vector is thus c(0, 0.60, 0.70, 0.80, 0.90, 1.00)

Now, we will plot the choropleth map by using the code chunk below.

tm_shape(mpsz_agemale2018)+
  tm_fill("DEPENDENCY",
          breaks = c(0, 0.60, 0.70, 0.80, 0.90, 1.00)) +
  tm_borders(alpha = 0.5)
## Warning: The shape mpsz_agemale2018 is invalid. See sf::st_is_valid

Colour Scheme

tmap supports colour ramps either defined by the user or a set of predefined colour ramps from the RColorBrewer package ### Using colour brewer palatte To change the colour, we assign the preferred colour to palette argument of tm_fill()

tm_shape(mpsz_agemale2018)+
  tm_fill("DEPENDENCY",
          n = 6,
          style = "quantile",
          palette = "Blues") +
  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.

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.

Map Legends

tm_shape(mpsz_agemale2018)+
  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)
## Warning: The shape mpsz_agemale2018 is invalid. See sf::st_is_valid

# Map style

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

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(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"))
## Warning: The shape mpsz_agemale2018 is invalid. See sf::st_is_valid

To reset data style

tmap_style("white") 
## tmap style set to "white"
## other available styles are: "gray", "natural", "cobalt", "col_blind", "albatross", "beaver", "bw", "classic", "watercolor"

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: 1. By assigning multiple values to at least one of the asthetic arguments 2. By defining a group-by variable in tm_facets() 3. By creating multiple stand-alone maps with tmap_arrange().

By assigning multiple values to at least one of the aesthetic arguments

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","AGED"),
          style = c("equal", "quantile"), 
          palette = list("Blues","Greens")) +
  tm_layout(legend.position = c("right", "bottom"))
## Warning: The shape mpsz_agemale2018 is invalid. See sf::st_is_valid

By defining a group-by variable in 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)
## Warning: The argument drop.shapes has been renamed to drop.units, and is
## therefore deprecated
## Warning: The shape mpsz_agemale2018 is invalid. See sf::st_is_valid

By creating multiple stand alone maps with tp_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

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
## Warning in preprocess_gt(x, interactive = interactive, orig_crs =
## gm$shape.orig_crs): legend.width controls the width of the legend within a map.
## Please use legend.outside.size to control the width of the outside legend