International Geomagnetic Reference Field (IGRF); Movement of Magnetic North and South Pole Locations (1590–2025); Globe

# Load necessary libraries
pacman::p_load(pacman, readr, tidyverse, ggplot2, viridis, maps, mapproj)

# Function to load, process, and split the dataset
load_and_process_data <- function(file_path) {
  # Load the dataset
  data <- read.csv(file_path, header = FALSE, col.names = "Values")
  
  # Split the 'Values' column into Longitude, Latitude, and Year
  data <- data %>%
    mutate(
      Longitude = as.numeric(sub(" .*", "", Values)),
      Latitude = as.numeric(sub(".* (-?\\d+\\.\\d+) .*", "\\1", Values)),
      Year = as.numeric(sub(".* (-?\\d+\\.\\d+)$", "\\1", Values))
    ) %>%
    select(-Values)
  
  return(data)
}

# Function to create a polar plot for North and South Pole data
create_polar_plot <- function(pole_data, pole_name, orientation) {
  world <- map_data("world")  # Get world map data
  
  ggplot() +
    geom_polygon(data = world, aes(x = long, y = lat, group = group), fill = "#cccccc", color = "#000000") +  # Add world map
    geom_point(data = pole_data, aes(x = Longitude, y = Latitude, color = Year), alpha = 0.7, size = 2) +  # Add pole locations
    scale_color_viridis(option = "turbo", guide = guide_colorbar(title = "Year")) +  # Viridis colour scale
    labs(
      title = paste("International Geomagnetic Reference Field (IGRF); Movement of Magnetic", pole_name, "Pole Locations (1590–2025); Geographic context"),
      subtitle = "Observed north and south dip poles during 1831–2007. Modelled pole locations from 1590–2025",
      caption = "Source: International Association of Geomagnetism and Aeronomy (IAGA)"
    ) +
    theme_minimal() +
    theme(
      legend.position = "right",  # Move legend to the right
      axis.text.x = element_blank(),  # Remove longitude labels
      axis.text.y = element_blank(),  # Remove latitude labels
      axis.ticks = element_blank(),  # Remove axis ticks
      panel.grid.major = element_blank(),  # Remove major gridlines
      panel.grid.minor = element_blank()   # Remove minor gridlines
    ) +
    coord_map("ortho", orientation = orientation)  # Orthographic projection centered on the pole
}

# Load and process North Pole data
north_pole_data <- load_and_process_data("IGRF_North_Pole_Locations_(1590–2025).csv")

# Plot North Pole data with a polar projection
create_polar_plot(north_pole_data, "North", orientation = c(90, 0, 0))

# Load and process South Pole data
south_pole_data <- load_and_process_data("IGRF_South_Pole_Locations_(1590–2025).csv")

# Plot South Pole data with a polar projection
create_polar_plot(south_pole_data, "South", orientation = c(-90, 0, 0))