#30DayMapChallenge 2022 - Day 06 - Networks

pacific_network.R

# Map of South Pacific islands, drawing a line from each to the nearest larger island
# David Friggens, November 2022

# Adapting https://gist.github.com/dakvid/982a017f46556806e33fd3b520877ce2

library(dplyr)
library(tidyr)
library(sf)
library(ggplot2)
library(sysfonts)
library(showtext)
showtext_auto()
library(ragg)

font_add_google("Pacifico", "Pacifico")
CHART_FONT <- "Pacifico"

# Read data ---------------------------------------------------------------

# from https://pacificdata.org/data/dataset/pacific-island-region-spatial-data48d9036d-04cb-4cb0-8aa8-b16197985b6a
pacific_islands <-
  st_read("Pacific islands region land.geojson")

pacific_points <- 
  pacific_islands |> 
  # essential to not split the Pacific in the middle!
  st_shift_longitude() |> 
  mutate(land_area = geometry |> st_area()) |> 
  st_centroid() |> 
  arrange(land_area)

pacific_nodes <- 
  pacific_points |> 
  st_drop_geometry()


# Compute connections -----------------------------------------------------

# NB The naive "cartesion product and then filter" approach I've used previously
#    is not really practical here. I made it work with some decomposition, but
#    I've got a lot of RAM on this machine, so more effort may be needed elsewhere.

pacific_connections <- 
  # pairwise joining of all islands (non-geo) ... in chunks ...
  inner_join(
    pacific_nodes |> 
      slice_head(n = 3000) |> 
      select(small_fid = fid,
             small_area = land_area),
    pacific_nodes |> 
      slice_head(n = 3000) |> 
      select(large_fid = fid,
             large_area = land_area),
    by = character(0)
  ) |> 
  # reduce to ordered pairs
  filter(small_area < large_area) |> 
  # ... and keep repeating to cut down the permutations - probably a better way to do this...
  bind_rows(
    inner_join(
      pacific_nodes |> 
        slice_head(n = 3000) |> 
        select(small_fid = fid,
               small_area = land_area),
      pacific_nodes |> 
        slice_tail(n = -3000) |> 
        select(large_fid = fid,
               large_area = land_area),
      by = character(0)
    ) |> 
      filter(small_area < large_area)
  ) |> 
  bind_rows(
    inner_join(
      pacific_nodes |> 
        slice_head(n = 6000) |> slice_tail(n = 3000) |> 
        select(small_fid = fid,
               small_area = land_area),
      pacific_nodes |> 
        slice_head(n = 6000) |> slice_tail(n = 3000) |> 
        select(large_fid = fid,
               large_area = land_area),
      by = character(0)
    ) |> 
      filter(small_area < large_area)
  ) |> 
  bind_rows(
    inner_join(
      pacific_nodes |> 
        slice_head(n = 6000) |> slice_tail(n = 3000) |> 
        select(small_fid = fid,
               small_area = land_area),
      pacific_nodes |> 
        slice_tail(n = -6000) |> 
        select(large_fid = fid,
               large_area = land_area),
      by = character(0)
    ) |> 
      filter(small_area < large_area)
  ) |> 
  bind_rows(
    inner_join(
      pacific_nodes |> 
        slice_head(n = 9000) |> slice_tail(n = 3000) |> 
        select(small_fid = fid,
               small_area = land_area),
      pacific_nodes |> 
        slice_head(n = 9000) |> slice_tail(n = 3000) |> 
        select(large_fid = fid,
               large_area = land_area),
      by = character(0)
    ) |> 
      filter(small_area < large_area)
  ) |> 
  bind_rows(
    inner_join(
      pacific_nodes |> 
        slice_head(n = 9000) |> slice_tail(n = 3000) |> 
        select(small_fid = fid,
               small_area = land_area),
      pacific_nodes |> 
        slice_tail(n = -9000) |> 
        select(large_fid = fid,
               large_area = land_area),
      by = character(0)
    ) |> 
      filter(small_area < large_area)
  ) |> 
  bind_rows(
    inner_join(
      pacific_nodes |> 
        slice_head(n = 12000) |> slice_tail(n = 3000) |> 
        select(small_fid = fid,
               small_area = land_area),
      pacific_nodes |> 
        slice_head(n = 12000) |> slice_tail(n = 3000) |> 
        select(large_fid = fid,
               large_area = land_area),
      by = character(0)
    ) |> 
      filter(small_area < large_area)
  ) |> 
  bind_rows(
    inner_join(
      pacific_nodes |> 
        slice_head(n = 12000) |> slice_tail(n = 3000) |> 
        select(small_fid = fid,
               small_area = land_area),
      pacific_nodes |> 
        slice_tail(n = -12000) |> 
        select(large_fid = fid,
               large_area = land_area),
      by = character(0)
    ) |> 
      filter(small_area < large_area)
  ) |> 
  bind_rows(
    inner_join(
      pacific_nodes |> 
        slice_tail(n = -12000) |> 
        select(small_fid = fid,
               small_area = land_area),
      pacific_nodes |> 
        slice_tail(n = -12000) |> 
        select(large_fid = fid,
               large_area = land_area),
      by = character(0)
    ) |> 
      filter(small_area < large_area)
  )
  
pacific_network <- 
  pacific_connections |> 
  # add coordinates
  inner_join(
    pacific_points |> 
      select(small_fid = fid,
             small_point = geometry),
    by = "small_fid"
  ) |> 
  inner_join(
    pacific_points |> 
      select(large_fid = fid,
             large_point = geometry),
    by = "large_fid"
  ) |> 
  # find closest large neighbour
  mutate(distance = st_distance(small_point, large_point,
                                by_element = TRUE)) |> 
  group_by(small_fid) |> 
  filter(distance == min(distance)) |> 
  ungroup() |> 
  # convert points to lines
  pivot_longer(cols = c(small_point, large_point),
               values_to = "geom") |> 
  st_as_sf() |> 
  group_by(small_fid, small_area,
           large_fid, large_area,
           distance) |> 
  summarise(do_union = FALSE) |> 
  ungroup() |> 
  st_cast("LINESTRING") |> 
  mutate(distance = as.numeric(distance))

# saveRDS(pacific_network, "pacific_network.rds")
# pacific_network <- readRDS("pacific_network.rds")

# Plot --------------------------------------------------------------------

LIGHT_BLUE <- "#009fde"
DARK_BLUE <- "#031f79"

agg_png(filename = "Day_06_Network.png",
        width = 1690, height = 1024,
        background = LIGHT_BLUE)
ggplot() +
  geom_sf(data = pacific_network |> st_shift_longitude(),
          # aes(colour = log(distance)),
          colour = DARK_BLUE,
          size = 1,
          alpha = 0.6,
          show.legend = FALSE) +
  coord_sf(datum = NA) +
  labs(title = "South Pacific Islands: Connections to the nearest larger island",
       subtitle = "#30DayMapChallenge 2022 - Day 06 - Networks",
       caption = "Data: SPREP Environmental Monitoring and Governance. Map: David Friggens, @dakvid") +
  theme_void() +
  theme(plot.background = element_rect(fill = LIGHT_BLUE, colour = LIGHT_BLUE),
        panel.background = element_rect(fill = LIGHT_BLUE, colour = LIGHT_BLUE),
        plot.title = element_text(size = 38, family = CHART_FONT, colour = DARK_BLUE, hjust = 0.1),
        plot.subtitle = element_text(size = 14, family = CHART_FONT, colour = DARK_BLUE, hjust = 0.03),
        plot.caption = element_text(size = 12, family = CHART_FONT, colour = DARK_BLUE, hjust = 0.95))
dev.off()