Compass Declination - GSMNP

# Load necessary libraries
pacman::p_load(pacman, readr, tidyverse, ggplot2, dplyr, gridExtra, sf, akima, rnaturalearth, rnaturalearthdata)

# Load the dataset
data <- read.csv("GRSM_COMPASS_DECLINATION_5646692726181302872.csv")

# Convert the data into an sf object
coordinates <- st_as_sf(data, coords = c("LON", "LAT"), crs = 4326)

# Get natural earth data for basemaps (filter for North America region)
world <- ne_countries(scale = "medium", continent = "North America", returnclass = "sf")

# Visualise the magnetic declination for GSMNP; North American Map
ggplot(data = world) +
  geom_sf(fill = "#cccccc", color = "#ffffff") +
  geom_sf(data = coordinates, aes(color = DECLINATION), size = 1, alpha = 0.9) +
  scale_color_viridis_c(option = "C", name = "Magnetic Declination (°)") +
  theme_minimal() +
  labs(title = "Magnetic Declination in Great Smoky Mountains National Park (GSMNP); North America",
       subtitle = "Data visualised using rnaturalearthdata and ggplot2",
       caption = "Source: GRSM GIS; National Park Service; 27 February 2017; Data updated 6 May 2021") +
  coord_sf(xlim = c(-130, -60), ylim = c(20, 60)) # Focus on North America

# Visualise the annual drift for GSMNP; North American Map
ggplot(data = world) +
  geom_sf(fill = "#cccccc", color = "#ffffff") +
  geom_sf(data = coordinates, aes(color = ANNUAL_DRIFT), size = 1, alpha = 0.9) +
  scale_color_viridis_c(option = "D", name = "Annual Drift (°/yr)") +
  theme_minimal() +
  labs(title = "Annual Drift of Magnetic Declination in Great Smoky Mountains National Park (GSMNP); North America",
       subtitle = "Data visualised using rnaturalearthdata and ggplot2",
       caption = "Source: GRSM GIS; National Park Service; 27 February 2017; Data updated 6 May 2021") +
  coord_sf(xlim = c(-130, -60), ylim = c(20, 60)) # Focus on North America

# Get natural earth data for basemaps
world <- ne_countries(scale = "medium", returnclass = "sf")

# Define latitude and longitude range for GSMNP
lon_min <- -84
lon_max <- -82
lat_min <- 35.4
lat_max <- 36.1

# Visualise magnetic declination for GSMNP
ggplot(data = world) +
  geom_sf(fill = "#cccccc", color = "#ffffff") +
  geom_sf(data = coordinates, aes(color = DECLINATION), size = 1.5, alpha = 0.9) +
  scale_color_viridis_c(option = "C", name = "Magnetic Declination (°)") +
  theme_minimal() +
  labs(title = "Magnetic Declination in Great Smoky Mountains National Park (GSMNP)",
       subtitle = "Data visualised using rnaturalearthdata and ggplot2",
       caption = "Source: GRSM GIS; National Park Service; 27 February 2017; Data updated 6 May 2021") +
  coord_sf(xlim = c(lon_min, lon_max), ylim = c(lat_min, lat_max)) # Zoom in on GSMNP

# Visualise annual drift for GSMNP
ggplot(data = world) +
  geom_sf(fill = "#cccccc", color = "#ffffff") +
  geom_sf(data = coordinates, aes(color = ANNUAL_DRIFT), size = 1.5, alpha = 0.9) +
  scale_color_viridis_c(option = "D", name = "Annual Drift (°/yr)") +
  theme_minimal() +
  labs(title = "Annual Drift of Magnetic Declination in Great Smoky Mountains National Park (GSMNP) ",
       subtitle = "Data visualised using rnaturalearthdata and ggplot2",
       caption = "Source: GRSM GIS; National Park Service; 27 February 2017; Data updated 6 May 2021") +
  coord_sf(xlim = c(lon_min, lon_max), ylim = c(lat_min, lat_max)) # Zoom in on GSMNP

# Interpolate magnetic declination for contours
declination_interp <- with(data, interp(x = LON, y = LAT, z = DECLINATION, duplicate = "mean"))

# Convert interpolated data to a data frame
declination_interp_df <- as.data.frame(expand.grid(x = declination_interp$x, y = declination_interp$y))
declination_interp_df$DECLINATION <- as.vector(declination_interp$z)

# Remove rows with missing declination values
declination_interp_df <- declination_interp_df[!is.na(declination_interp_df$DECLINATION), ]

# Plot with contours
ggplot() +
  # Base map
  geom_sf(data = world, fill = "#cccccc", color = "#ffffff") +
  # Contour lines for magnetic declination
  geom_contour(data = declination_interp_df, aes(x = x, y = y, z = DECLINATION), color = "#3333ff", bins = 10) +
  # Magnetic declination points with elevation
  geom_sf(data = coordinates, aes(color = DECLINATION, size = ELEVATION), alpha = 0.8) +
  # Colour scale for declination
  scale_color_viridis_c(option = "C", name = "Magnetic Declination (°)") +
  # Size scale for elevation
  scale_size_continuous(name = "Elevation (m)", range = c(0.5, 3)) +
  # Minimal theme
  theme_minimal() +
  # Labels and captions
  labs(
    title = "Magnetic Declination and Elevation in Great Smoky Mountains National Park; GSMNP",
    subtitle = "Contour lines represent interpolated magnetic declination",
    caption = "Source: GRSM GIS; National Park Service; 27 February 2017; Data updated 6 May 2021"
  ) +
  # Coordinate limits for GSMNP
  coord_sf(xlim = c(lon_min, lon_max), ylim = c(lat_min, lat_max))

## Warning: Removed 90 rows containing missing values or values outside the scale range
## (`geom_sf()`).