Skip to content

Instantly share code, notes, and snippets.

@jonathonmellor
Created April 19, 2023 17:26
Show Gist options
  • Save jonathonmellor/a568381d338fa281f5cc9c2a91b6be13 to your computer and use it in GitHub Desktop.
Save jonathonmellor/a568381d338fa281f5cc9c2a91b6be13 to your computer and use it in GitHub Desktop.
wd <- system("echo $(git rev-parse --show-toplevel)", intern = TRUE)
setwd(wd)
source(paste0(wd, "/influenza/scripts/depends.R"))
output_dir <- fs::path(paste0(wd, "/influenza/papers/univariate/outputs"))
# nhs_data path for getting population sizes:
nhs_data_path <- "" # Removed path to influenza data + population catchments per trust
# generate population size for each STP
popsizes <- aws.s3::s3read_using(
vroom::vroom,
object = nhs_data_path
) %>%
select(trust_code, trust_name, stp, nhs_region, population_size) %>%
distinct() %>%
mutate(stp = tolower(stp)) %>%
mutate(
stp = gsub(" \\(stp)", "", stp),
stp = gsub(" stp", "", stp),
stp = gsub("&", "and", stp),
stp = gsub(" \\()", "", stp),
nhs_region = gsub("Of", "of", nhs_region),
nhs_region = gsub("And", "and", nhs_region)
)
# Loading look ups for stp boundaries, regions and trust codes
STP <- s3read_using(
object = "", # Removed path to ONS' NHS STP Boundaries
FUN = unzip,
exdir = tempdir()
) %>%
str_subset(".shp$") %>%
st_read() %>%
rename(stp_name = STP21NM) %>%
select(stp_name, geometry) %>%
arrange(stp_name) %>%
mutate(stp_name = tolower(stp_name))
regions <- s3read_using(
object = "", # Removed path to ONS' NHS Commissioning Region Boundaries
FUN = unzip,
exdir = tempdir()
) %>%
str_subset(".shp$") %>%
st_read() %>%
rename(nhs_region = NHSER21NM) %>%
select(nhs_region, geometry) %>%
arrange(nhs_region)
Trust_coordinates <- aws.s3::s3read_using(
data.table::fread,
object = "", # Removed path to IDM team custom Trust coordinates (ets, etr derived)
na.string = c("NA", "NULL")
) %>% filter(trust_code %in% unique(popsizes$trust_code))
# join spatial frames to population sizes
popsizes_full <- popsizes %>%
group_by(stp) %>%
summarise(nhs_region,
stp_popsizes = sum(population_size)
) %>%
group_by(nhs_region) %>%
summarise(stp, stp_popsizes,
region_popsizes = sum(stp_popsizes)
) %>%
full_join(STP, by = c("stp" = "stp_name")) %>%
full_join(regions, by = c("nhs_region"), suffix = c("_stp", "_region")) %>%
distinct()
# trust coordinates, convert to spatial points
trust_coords <- st_as_sf(Trust_coordinates,
coords = c("longitude", "latitude"),
crs = 4326, agr = "constant"
)
# convert the crs to match stps
trusts_regions <- popsizes %>%
select(trust_code, nhs_region)
# bring crs in line between boundaries and coords
trust_coords <- st_transform(trust_coords, crs = 27700) %>%
left_join(trusts_regions, by = c("trust_code"))
labels_legend <- c("0 - 625,000", "625,000 - 1,250,000", "1,250,000 - 1,876,000", "1,876,000 - 2,501,000", "2,501,000 - 3,126,000")
# dataframe with binned population sizes
binned <- popsizes_full %>%
ungroup() %>%
mutate(
stp_popsizes_1000 = stp_popsizes / 1000,
bins = cut(stp_popsizes_1000, seq(0, max(stp_popsizes_1000), length.out = 6), dig.lab = 4, labels = labels_legend)
)
# ordering bins correctly
binned$bins <- fct_rev(binned$bins)
# plotting London
london_binned <- ggplot(filter(binned, nhs_region == "London")) +
geom_sf(aes(geometry = geometry_stp, fill = bins), alpha = 0.5, linewidth = 0.2) +
geom_sf(aes(geometry = geometry_region), color = "black", fill = NA, linewidth = 0.75) +
scale_fill_viridis_d(
option = "D",
direction = 1,
drop = FALSE
) +
theme_void() +
geom_sf(
# problem with how coordinates have been assigned - moorefield eye appears outside london (many sites)
data = filter(trust_coords, nhs_region == "London", trust_code != "RP6"), aes(geometry = geometry),
fill = "white", colour = "black", shape = 21, size = 1.5
) +
labs(
fill = "Population Size",
title = "London"
) +
theme(
legend.position = "none",
plot.title = element_text(hjust = 0.5, size = 14),
plot.background = element_rect(fill = "white", colour = "black")
)
london_binned <- ggplotGrob(london_binned)
# plotting england + london
binned_map <- ggplot() +
geom_sf(
data = binned, aes(geometry = geometry_stp, fill = bins, linewidth = "Sustainability and \nTransformation Partnership"), alpha = 0.5, # linewidth = 0.2
) +
geom_sf(
data = binned, aes(geometry = geometry_region, linewidth = "NHS Region"), color = "black", fill = NA, # linewidth = 0.75
) +
scale_fill_viridis_d(
option = "D",
direction = 1
) +
scale_discrete_manual("linewidth",
values = c("NHS Region" = 0.75, "Sustainability and \nTransformation Partnership" = 0.2),
guide = guide_legend(override.aes = list(linewidth = c(0.75, 0.2)))
) +
theme_void() +
geom_sf(
data = trust_coords, aes(geometry = geometry, shape = "NHS Trust"),
fill = "white", colour = "black", size = 1.5
) +
scale_discrete_manual("shape",
values = c("NHS Trust" = 21),
guide = guide_legend(override.aes = list(
shape = c(21),
size = 3
))
) +
labs(
fill = "Population Size",
linewidth = NULL,
# title = "Regional, STP, and NHS Trust breakdowns of England",
shape = NULL
) +
theme(
plot.title = element_text(hjust = 0.5),
plot.background = element_rect(fill = "white", colour = "white"),
legend.key.size = unit(1, "cm"),
legend.text = element_text(size = 12),
legend.title = element_text(size = 14)
) +
annotation_custom(
grob = london_binned,
xmin = 550000,
xmax = 750000,
ymin = 500000,
ymax = 670000
)
ggsave(
plot = binned_map,
filename = paste0(output_dir, "/binned_map.png"),
width = 12, height = 10, type = "cairo"
)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment