library(tidyverse)
library(rnaturalearth)
library(sf)
library(ggrepel)
library(glue)

north_pole_file <- "NP.xy"
if (!file.exists(north_pole_file)) {
  download.file("https://www.ngdc.noaa.gov/geomag/data/poles/NP.xy",
                north_pole_file)
}

np <- read_delim(north_pole_file, 
                 delim = " ", 
                 col_names = c("x", "y", "year")) |> 
  filter(year <= year(Sys.Date())) |> 
  mutate(lon = if_else(x <= 180, x, x - 360),
         lat = y) |> 
  st_as_sf(coords = c("lon", "lat"), 
           crs = "EPSG:4326") |> 
  st_transform("EPSG:3995")

countries <- ne_countries(scale = 50, returnclass = "sf") |> 
  st_transform("EPSG:3995")

bbox <- np |> 
  st_buffer(1.2e6) |> 
  st_bbox(np)
  
np |> 
  ggplot() +
  geom_sf(data = countries, fill = "wheat3", color = "wheat3") +
  geom_sf(aes(color = year), size = 1.2) +
  geom_text_repel(aes(label = if_else(!year %% 100, year, NA), geometry = geometry),
                  bg.color = "white",
                  bg.r = 0.1,
                  stat = "sf_coordinates", 
                  size = 3) +
  scale_color_viridis_c() +
  coord_sf(xlim = bbox[c(1, 3)], 
           ylim = bbox[c(2, 4)]) +
  labs(title = "Wandering of the Geomagnetic North Pole",
       subtitle = glue("1590-{max(np$year)}"),
       x = "",
       y = "", 
       caption = glue("Data: NOAA - IGRF Pole Locations
                       Basemap: Natural Earth
                       http://r.iresmi.net/ {Sys.Date()}")) +
  theme(legend.position = "none",
        plot.caption = element_text(size = 6),
        panel.background = element_rect(fill = "azure2"),
        panel.grid = element_line(color = "paleturquoise3"))