#30DayMapChallenge 2021 - Day 20 - Movement

hex_earthquake_animation.R

# Create an animated heatmap grid of earthquakes in New Zealand
# For #30DayMapChallenge 2021 - Day 20 - Movement
#  -- David Friggens, November 2021


# Setup -------------------------------------------------------------------

library(glue)
library(readr)
library(dplyr)
library(tidyr)
library(lubridate)
library(purrr)
library(sf)
library(rmapshaper)
library(ggplot2)
library(av)

FONT = "Wallpoet"

# Data --------------------------------------------------------------------

# from https://datafinder.stats.govt.nz
nz <- 
  read_sf("statsnz/territorial-authority-2021-clipped-generalised.gpkg") %>% 
  filter(! TA2021_V1_00 %in% c("067", "999")) %>% 
  summarise(n = n()) %>% 
  ms_simplify()

# from https://quakesearch.geonet.org.nz
eq <-
  read_csv("data/earthquakes.csv") %>% 
  filter(eventtype == "earthquake") %>% 
  transmute(longitude,
            latitude,
            magnitude,
            year = year(origintime),
            month = month(origintime, label = TRUE, abbr = TRUE))

data_dates <- 
  eq %>% 
  distinct(year, month)

geo_eq <- 
  eq %>% 
  st_as_sf(coords = c("longitude", "latitude"), crs = 4326L) %>% 
  st_transform(crs = st_crs(nz))


# Grids -------------------------------------------------------------------

geo_grid <- 
  geo_eq %>% 
  st_make_grid(n = c(50, 60),
               square = FALSE,
               flat_topped = TRUE) %>% 
  st_as_sf() %>% 
  mutate(cell = row_number())

grid_nz <- 
  geo_grid %>% 
  filter(st_intersects(geo_grid, nz, sparse = FALSE)) %>% 
  summarise(n = n())

grid_eq <- 
  geo_grid %>% 
  st_join(geo_eq) %>% 
  st_drop_geometry() %>% 
  group_by(year, month, cell) %>%
  summarise(magnitude = log(sum(10^magnitude))) %>% 
  ungroup()

grid_heat <- 
  geo_grid %>% 
  inner_join(data_dates,
             by = character()) %>% 
  left_join(grid_eq,
             by = c("year", "month", "cell"))

max_magnitude <- 
  grid_heat %>% 
  pull(magnitude) %>% 
  max(na.rm = TRUE) %>% 
  round(1)


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

grids_by_month <- 
  grid_heat %>% 
  nest(earthquakes = c(magnitude, cell, x))

plot_earthquakes <- 
  function(year, month, earthquakes) {
    gg_one_month <- 
      ggplot() +
      geom_sf(data = earthquakes,
              aes(fill = magnitude),
              colour = NA) +
      scale_fill_viridis_c(na.value = "white",
                           limits = c(0, max_magnitude + 0.1),
                           breaks = c(5, 10, 15)) +
      geom_sf(data = grid_nz,
              colour = "black", size = 1,
              fill = NA) +
      coord_sf(datum = NA) +
      labs(title = glue("NZ Earthquakes, 4M+, {year} {month}"),
           subtitle = "#30DayMapChallenge 2021 - Day 20 - Movement",
           caption = "Data: GeoNet. Map: David Friggens, @dakvid",
           fill = "Combined\nMagnitude") +
      theme_minimal(base_family = FONT) +
      theme(legend.position = c(0.85, 0.15),
            plot.background = element_rect(fill = "white", colour = "white"),
            panel.background = element_rect(fill = "white", colour = "white"),
            plot.title = element_text(size = 56, margin = margin(t = 30), hjust = 0.1),
            plot.subtitle = element_text(size = 36, hjust = 0.1),
            plot.caption = element_text(size = 32, margin = margin(b = 20), hjust = 0.88),
            legend.title = element_text(size = 28),
            legend.text = element_text(size = 28),
            legend.key.height = unit(60, "points"),
            legend.key.width = unit(20, "points"))
    
    ggsave(plot = gg_one_month,
           glue("Day_20/earthquakes_{year}_{if_else(month < 'Oct', '0', '')}{as.integer(month)}.png"),
           width = 20, height = 26, dpi = 72)
  }

grids_by_month %>%
  arrange(year, month) %>%
  pwalk(plot_earthquakes)


# Empty grid to go at the end, for looping

gg_empty <- 
  ggplot() +
  geom_sf(data = geo_grid %>% mutate(magnitude = NA_real_),
          aes(fill = magnitude),
          colour = NA) +
  scale_fill_viridis_c(na.value = "white",
                       limits = c(0, max_magnitude + 0.1),
                       breaks = c(5, 10, 15)) +
  geom_sf(data = grid_nz,
          colour = "black", size = 1,
          fill = NA) +
  coord_sf(datum = NA) +
  labs(title = "NZ Earthquakes, 4M+,",
       subtitle = "#30DayMapChallenge 2021 - Day 20 - Movement",
       caption = "Data: GeoNet. Map: David Friggens, @dakvid",
       fill = "Combined\nMagnitude") +
  theme_minimal(base_family = FONT) +
  theme(legend.position = c(0.85, 0.15),
        plot.background = element_rect(fill = "white", colour = "white"),
        panel.background = element_rect(fill = "white", colour = "white"),
        plot.title = element_text(size = 56, margin = margin(t = 30), hjust = 0.1),
        plot.subtitle = element_text(size = 36, hjust = 0.1),
        plot.caption = element_text(size = 32, margin = margin(b = 20), hjust = 0.88),
        legend.title = element_text(size = 28),
        legend.text = element_text(size = 28),
        legend.key.height = unit(60, "points"),
        legend.key.width = unit(20, "points"))

ggsave(plot = gg_empty,
       "Day_20/earthquakes_9999.png",
       width = 20, height = 26, dpi = 72)

av_encode_video(input = list.files(path = "Day_20", pattern = "earthquakes.*png", full.names = TRUE),
                output = "Day_20/Day_20_earthquakes.mp4",
                framerate = 2)