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

# 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 2D plot for North and South Pole data
create_pole_plot <- function(pole_data, pole_name) {
  ggplot(pole_data, aes(x = Longitude, y = Latitude)) +
    geom_point(aes(color = Year), alpha = 0.7) +  # Only colour by Year, no size variation
    scale_color_viridis(option = "turbo", guide = guide_colorbar(title = "Year")) +  # Viridis colour scale
    theme_minimal() +
    labs(
      title = paste("International Geomagnetic Reference Field (IGRF); Movement of Magnetic", pole_name, "Pole Locations (1590–2025); Poles trajectory"),
      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)",
      x = "Longitude", y = "Latitude"
    ) +
    theme(legend.position = "right")  # Move legend to the right for better visibility
}

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

# Plot North Pole data
create_pole_plot(north_pole_data, "North")

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

# Plot South Pole data
create_pole_plot(south_pole_data, "South")