Let's Make a Map in R - Medium post (20 February 2018)

script.R

## Let’s make a map in R ##

# Setup --------------------------- 
install.packages("sf", "tidyverse", "classInt", "viridis")
library(sf) ; library(tidyverse) ; library(classInt) ; library(viridis)

# Reading spatial data --------------------------- 
sf_gb <- st_read("https://opendata.arcgis.com/datasets/07194e4507ae491488471c84b23a90f2_3.geojson", quiet = TRUE)
glimpse(sf_gb)
st_geometry(sf_gb)
plot(st_geometry(sf_gb))

# Filtering spatial data --------------------------- 
lookup <- read_csv("https://opendata.arcgis.com/datasets/046394602a6b415e9fe4039083ef300e_0.csv") %>%
  filter(LAD17NM %in% c("Bolton","Bury","Manchester","Oldham","Rochdale","Salford","Stockport","Tameside","Trafford","Wigan")) %>%
  pull(WD17CD)
sf_gm <- sf_gb %>%
  filter(wd17cd %in% lookup)
plot(st_geometry(sf_gm))
sf_gm <- sf_gm %>%
  select(area_code = wd17cd, area_name = wd17nm)

# Joining spatial data --------------------------- 
df <- read_csv("data/bulk.csv")
df_census <- df %>%
  select(area_code = `geography code`,
         n = `Qualification: No qualifications; measures: Value`,
         total = `Qualification: All usual residents aged 16 and over; measures: Value`) %>%
  mutate(percent = (n/total)*100)
sf_gm_census <- left_join(sf_gm, df_census, by = "area_code")

# Creating a choropleth map --------------------------- 
# equal intervals
ggplot() +
  geom_sf(data = sf_gm_census,
          aes(fill = cut_number(percent, 5)),
          alpha = 0.8,
          colour = 'white',
          size = 0.3) +
  scale_fill_brewer(palette = "PuBu",
                    name = "No qualifications (%)") +
  labs(x = NULL, y = NULL,
       title = "Residents with no qualifications in Greater Manchester, 2011",
       subtitle = "Source: Table QS502EW, Census 2011",
       caption = "Contains OS data © Crown copyright and database right (2018)") +
  theme(line = element_blank(),
        axis.text = element_blank(),
        axis.title = element_blank(),
        panel.background = element_blank()) +
  coord_sf(datum = NA)

# natural breaks
classes <- classIntervals(sf_gm_census$percent, n = 5, style = "jenks")
classes$brks
sf_gm_census <- sf_gm_census %>%
  mutate(percent_class = cut(percent, classes$brks, include.lowest = T))

ggplot() +
  geom_sf(data = sf_gm_census,
          aes(fill = percent_class),
          alpha = 0.8,
          colour = 'white',
          size = 0.3) +
  scale_fill_brewer(palette = "PuBu",
                    name = "No qualifications (%)") +
  labs(x = NULL, y = NULL,
       title = "Residents with no qualifications in Greater Manchester, 2011",
       subtitle = "Source: Table QS502EW, Census 2011",
       caption = "Contains OS data © Crown copyright and database right (2018)") +
  theme(line = element_blank(),
        axis.text = element_blank(),
        axis.title = element_blank(),
        panel.background = element_blank()) +
  coord_sf(datum = NA)

# unclassed
ggplot(sf_gm_census, aes(fill = percent)) +
  geom_sf(alpha = 0.8, colour = 'white', size = 0.3) +
  scale_fill_viridis(discrete = F,
                     name = "No qualifications (%)", 
                     direction = -1, 
                     guide = guide_colourbar(
                       direction = "horizontal",
                       barheight = unit(2, units = "mm"),
                       barwidth = unit(50, units = "mm"),
                       draw.ulim = F,
                       title.position = 'top',
                       title.hjust = 0.5,
                       label.hjust = 0.5)) +
  labs(x = NULL, y = NULL, 
       title = "Residents with no qualifications in Greater Manchester, 2011",
       subtitle = "Source: Table QS502EW, Census 2011",
       caption = "Contains OS data © Crown copyright and database right (2018)") +
  coord_sf(datum = NA) +
  theme(line = element_blank(),
        axis.text = element_blank(),
        axis.title = element_blank(),
        panel.background = element_blank(),
        legend.position = c(0.2, 0.09),
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 8))