1. 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. By the end of this hands-on 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.

2. Load Packages

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

3. Aspatial Data Wrangling

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

setwd("D:/Materials of Study/VA Hands-on_Ex09")
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()
## )
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

3.2 Create a sf data frame from an aspatial data frame

sgpools_sf <- st_as_sf(sgpools, 
                       coords = c("XCOORD", "YCOORD"),
                       crs= 3414)
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
## epsg (SRID):    3414
## proj4string:    +proj=tmerc +lat_0=1.366666666666667 +lon_0=103.8333333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +units=m +no_defs
## # 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]>

4. Draw proportional symbol map

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

4.1 Start with an interactive point symbol map

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

4.2 Make it proportional (add a size argument)

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

4.3 Give it color

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

4.4 Twin plot

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) # sync actions on twin plot
## Legend for symbol sizes not available in view mode.

5.Geospatial Data Wranglin

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

mpsz <- st_read(dsn = "data/geospatial", 
                layer = "MP14_SUBZONE_WEB_PL")
## Reading layer `MP14_SUBZONE_WEB_PL' from data source `D:\Materials of Study\VA Hands-on_Ex09\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_defs

5.2 Import ttribute data

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()
## )

5.3 Data Wrangling

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

5.4 Joining the attribute data and geospatial data

mpsz_agemale2018 <- left_join(mpsz, popagsex2018_male, 
                              by = c("SUBZONE_N" = "SZ"))

6. Choropleth Mapping Geospatial Data Useing 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.

6.1 Plot a choropleth map quickly by using qtm()

qtm(mpsz_agemale2018, fill = "DEPENDENCY")

6.2 Create a choropleth map by using tmap’s elements

6.2.1 base map

tm_shape(mpsz_agemale2018) +
  tm_polygons()

6.2.2 add attribute

tm_shape(mpsz_agemale2018)+
  tm_polygons("DEPENDENCY")

6.2.3 try tm_fill() and tm_border()

tm_shape(mpsz_agemale2018)+
  tm_fill("DEPENDENCY")
tm_shape(mpsz_agemale2018)+
  tm_fill("DEPENDENCY") +
  tm_borders(alpha = 0.5)

6.3 Data classification methods of tmap

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

6.3.1 Use ColourBrewer palette

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

6.4 Map layout

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.

6.4.1 Map Legend (the hist is not shows as expected in view mode)

tmap_mode("plot")
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)

6.4.2 Map style

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

6.4.3 Cartographic Furniture ()

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

6.5 Draw small multiple choropleth maps

6.5.1 split within tm_fill()

#tmap_mode("plot")
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")

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

6.5.2 split by 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)

6.5.3 split by 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)

6.6 Map Spatial Object Meeting a 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)

6.7 Create an interactive map application

tmap_mode("view")
mpsz_wgs84 <- st_transform(mpsz_agemale2018, 4326)
tm_shape(mpsz_wgs84)+
  tm_fill("DEPENDENCY",
          n = 6,
          style = "quantile", 
          palette = "Blues") +
  tm_borders(alpha = 0.5)