Skip to content

Instantly share code, notes, and snippets.

@rcatlord
Created March 3, 2023 13:27
Show Gist options
  • Save rcatlord/2e08a004874d46dba046f50d0a8d51cd to your computer and use it in GitHub Desktop.
Save rcatlord/2e08a004874d46dba046f50d0a8d51cd to your computer and use it in GitHub Desktop.
Census data for bespoke geographies
# Returning Census data for custom geographies #
# Load packages ---------------------
library(tidyverse)
library(sf)
library(mapedit)
library(mapview)
library(RColorBrewer)
# Retrieve data ---------------------
# MSOA names
# Source: House of Commons Library
# URL: https://houseofcommonslibrary.github.io/msoanames
msoa_names <- read_csv("https://houseofcommonslibrary.github.io/msoanames/MSOA-Names-Latest2.csv") %>%
select(msoa21cd, msoa21nm, msoa21hclnm, lad21nm = localauthorityname)
# MSOA boundaries
# Source: ONS Open Geography Portal
# URL: https://geoportal.statistics.gov.uk/datasets/ons::middle-layer-super-output-areas-december-2021-boundaries-super-generalised-clipped-ew-bsc
msoa <- st_read("https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/MSOA_2021_EW_BSC/FeatureServer/0/query?outFields=*&where=1%3D1&f=geojson") %>%
select(msoa21cd = MSOA21CD) %>%
left_join(msoa_names, by = "msoa21cd")
# Passports held: No passport held
# Source: 2021 Census, ONS
# URL: hhttps://www.nomisweb.co.uk/sources/census_2021_bulk
url <- "https://www.nomisweb.co.uk/output/census/2021/census2021-ts005.zip"
download.file(url, dest = "ensus2021-ts005.zip")
unzip("ensus2021-ts005.zip")
file.remove("ensus2021-ts005.zip")
df <- read_csv("census2021-ts005-msoa.csv") %>%
select(msoa21cd = `geography code`,
total = `Passports held: Total; measures: Value`,
n = `Passports held: No passport held; measures: Value`) %>%
mutate(percent = round(n/total*100,1)) %>%
select(msoa21cd, percent)
# Join Census data to spatial layer ---------------------
sf <- left_join(msoa, df, by = "msoa21cd")
# Draw a polygon or a rectangle ---------------------
box <- drawFeatures()
# Filter by chosen area and plot ---------------------
sf_use_s2(FALSE) # switch off s2 spherical geometry package
selected <- st_filter(sf, box, .predicate = st_intersects) # or st_within
# Create choropleth map ---------------------
mapview(selected, zcol = "percent",
col.regions = brewer.pal(5, "YlGnBu")) +
mapview(box, lwd = 2, alpha.regions = 0)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment