#30DayMapChallenge 2021 - Day 04 - Hexagons

hex_earthquakes.R

# Create a heatmap grid of earthquakes in New Zealand
# For #30DayMapChallenge 2021 - Day 04 - Hexagons
#  -- David Friggens, November 2021


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

library(readr)
library(dplyr)
library(sf)
library(rmapshaper)
library(ggplot2)

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") %>% 
  select(longitude, latitude, magnitude)

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(cell) %>%
  summarise(magnitude = log(sum(10^magnitude))) %>% 
  ungroup()

grid_heat <- 
  geo_grid %>% 
  inner_join(grid_eq,
             by = "cell")


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

gg_eq <- 
  ggplot() +
  geom_sf(data = grid_heat,
          aes(fill = magnitude),
          colour = NA) +
  scale_fill_viridis_c(na.value = "white") +
  geom_sf(data = grid_nz,
          colour = "black", size = 1,
          fill = NA) +
  coord_sf(datum = NA) +
  labs(title = "NZ Earthquakes, 4M+, 2006-2021",
       subtitle = "#30DayMapChallenge 2021 - Day 04 - Hexagons",
       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_eq,
       "Day_04/Day_04_earthquakes_log.png",
       width = 20, height = 26, dpi = 72)