Bivariate choropleth

Bivariate choropleth maps visualise the spatial relationship between two quantitative variables. The biscale package has a handy set of functions that allow you to create bivariate thematic maps in R with ease. The recipe below uses the biscale package to visualise the relationship between commuting to work by bike and overall physical activity in Manchester’s MSOAs.

Ingredients

Data sources
QS701EW, 2011 Census
Active Lives Survey, Sport England

R packages

Bivariate choropleth map

Instructions

Load the necessary R packages

library(tidyverse) ; library(httr) ; library(jsonlite) ; library(readxl) ; library(sf) ; library(biscale) ; library(cowplot)

Create a string object with the name of your local authority, e.g. Manchester

la <- "Manchester"

Retrieve the GSS codes for MSOAs in your chosen local authority

codes <- fromJSON(paste0("https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/OA11_LSOA11_MSOA11_LAD11_EW_LUv2/FeatureServer/0/query?where=LAD11NM%20like%20'%25", URLencode(toupper(la), reserved = TRUE), "%25'&outFields=MSOA11CD,LAD11NM&outSR=4326&f=json"), flatten = TRUE) %>% 
  pluck("features") %>% 
  as_tibble() %>% 
  distinct(attributes.MSOA11CD) %>% 
  pull(attributes.MSOA11CD)

Identify the percentage of residents that cycle to work

cycling <- read_csv(paste0("http://www.nomisweb.co.uk/api/v01/dataset/NM_568_1.data.csv?date=latest&geography=", paste0(codes, collapse=","), "&rural_urban=0&cell=0,9&measures=20100&select=geography_code,cell_name,obs_value")) %>% 
  spread(CELL_NAME, OBS_VALUE) %>%
  mutate(pct_cycling = round((Bicycle/`All categories: Method of travel to work`)*100,1)) %>%
  select(msoa11cd = GEOGRAPHY_CODE, pct_cycling)

Identify the percentage of resident doing at least 150 minutes of physical activity a week

tmp <- tempfile(fileext = ".xls")
GET(url = "https://www.sportengland.org/media/13326/als-2018-sae.xlsx",
    write_disk(tmp))
activity <- read_xlsx(tmp) %>% 
  filter(grepl(la, `Local Authority`, ignore.case = TRUE)) %>% 
  mutate(pct_active = round(`Active (%)`*100, 1)) %>% 
  select(msoa11cd = MSOA11, pct_active)

Retrieve the MSOA vector boundaries

msoa <- st_read(paste0("https://ons-inspire.esriuk.com/arcgis/rest/services/Census_Boundaries/Middle_Super_Output_Areas_December_2011_Boundaries/MapServer/2/query?where=", 
                        URLencode(paste0("msoa11cd IN (", paste(shQuote(codes), collapse = ", "), ")")), 
                        "&outFields=msoa11cd,msoa11nm&outSR=4326&f=geojson"))

Join data to the vector boundary layer

sf <- left_join(msoa, cycling, by = "msoa11cd") %>% 
  left_join(activity, by = "msoa11cd")

Create classes

data <- bi_class(sf, x = pct_cycling, y = pct_active, style = "quantile", dim = 3)

Create map

map <- ggplot() +
  geom_sf(data = data, mapping = aes(fill = bi_class), color = "white", size = 0.1, show.legend = FALSE) +
  bi_scale_fill(pal = "DkViolet", dim = 3) +
  labs(title = paste0("Cycling to work and physical activity by MSOA in ", la),
       caption = "Source: 2011 Census | Sport England") +
  coord_sf(crs = st_crs(4326), datum = NA) +
  theme_void() +
  theme(plot.margin = unit(c(1,1,1,1), "cm"),
        plot.title = element_text(size = 18, face = "bold", hjust = 0.5),
        plot.caption = element_text(size = 9, color = "grey50", hjust = 1, margin = margin(t = 15)))

Create legend

legend <- bi_legend(pal = "DkViolet", dim = 3, xlab = "Higher % cycling", ylab = "Higher % active", size = 8)

Combine the map with the legend. You might need to adjust the placement of the legend

plot <- ggdraw() +
  draw_plot(map, 0, 0, 1, 1) +
  draw_plot(legend, 0.5, .2, 0.2, 0.2)

Output the map as a PNG file

save_plot("bivariate_choropleth.png", plot, base_height = 8, base_aspect_ratio = 1.8)

Bivariate choropleth

Ingredients

Instructions

Further reading