Population map of NHS

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