Map connecting each New Zealand settlement with it's nearest largest neighbour

connections2.R

# Map connecting each New Zealand settlement with it's nearest largest neighbour
# For #30DayMapChallenge 2020, day 15 (connections)
# Inspired by https://twitter.com/AndriyYaremenko/status/1327885776047239168
# -- David Friggens, December 2020

library(dplyr)
library(stringr)
library(tidyr)
library(sf)
library(ggplot2)
library(hrbrthemes)

urban_rural_settlements <-
  read_sf("estimated-resident-population-at-30-june-2018-by-urban-rural.gpkg") %>% 
  filter(! IUR2020_V1_00_NAME %in% c("Inland water", "Oceanic", "Rural other")) %>% 
  mutate(island = if_else(UR2020_V1_00 %>% str_detect("^1"), "North Island", "South Island")) %>% 
  select(code = UR2020_V1_00,
         name = UR2020_V1_00_NAME,
         population = erp18,
         island) %>% 
  st_centroid()

settlement_connections <- 
  inner_join(
    urban_rural_settlements %>% 
      st_drop_geometry() %>% 
      select(small_code = code,
             small_population = population,
             island),
    urban_rural_settlements %>% 
      st_drop_geometry() %>% 
      select(big_code = code,
             big_population = population,
             island),
    by = character(0)
  ) %>% 
  filter(small_population < big_population) %>% 
  inner_join(
    urban_rural_settlements %>% 
      select(small_code = code,
             small_centre = geom),
    by = "small_code"
  ) %>% 
  inner_join(
    urban_rural_settlements %>% 
      select(big_code = code,
             big_centre = geom),
    by = "big_code"
  ) %>% 
  mutate(distance = st_distance(small_centre, big_centre,
                                by_element = TRUE)) %>%
  group_by(small_code) %>%
  filter(distance == min(distance)) %>%
  ungroup() %>% 
  pivot_longer(cols = c(small_centre, big_centre),
               values_to = "geom") %>% 
  st_as_sf() %>% 
  group_by(small_code, small_population,
           big_code, big_population,
           distance) %>% 
  summarise(do_union = FALSE) %>% 
  ungroup() %>% 
  st_cast("LINESTRING") %>% 
  mutate(distance = as.numeric(distance))

g_settlements <- 
  ggplot(settlement_connections) +
  geom_sf(aes(color = log(distance)),
          size = 1,
          alpha = 0.8,
          show.legend = FALSE) +
  coord_sf(datum = NA) +
  labs(title = "Nearest Bigger Settlement: NZ edition",
       subtitle = "#30DayMapChallenge 2020 - Day 15 - Connections",
       caption = "David Friggens, @dakvid") +
  theme_ipsum_rc()

ggsave(plot = g_settlements,
       path = ".", filename = "15_connections2.png", device = "png",
       width = 10.3, height = 12, dpi = 72)