projections.R
In geographiclib: Access to 'GeographicLib'

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(geographiclib)

## ----locations----------------------------------------------------------------
# Australian cities
australia <- cbind(
  lon = c(151.21, 144.96, 153.02, 115.86, 138.60),
  lat = c(-33.87, -37.81, -27.47, -31.95, -34.93)
)
rownames(australia) <- c("Sydney", "Melbourne", "Brisbane", "Perth", "Adelaide")

# Antarctic stations
antarctic <- cbind(
  lon = c(166.67, 77.97, 39.58, -64.05, 0),
  lat = c(-77.85, -68.58, -67.60, -64.25, -90)
)
rownames(antarctic) <- c("McMurdo", "Davis", "Mawson", "Palmer", "South Pole")

# World cities
world_pts <- cbind(
  lon = c(-0.13, -74.01, 139.69, 151.21, -43.17),
  lat = c(51.51, 40.71, 35.69, -33.87, -22.91)
)
rownames(world_pts) <- c("London", "New York", "Tokyo", "Sydney", "Rio")

## ----utm-basic----------------------------------------------------------------
# Convert Australian cities
utmups_fwd(australia)

## ----utm-zones----------------------------------------------------------------
# Points at different longitudes show different zones
lon_transect <- cbind(
  lon = seq(-180, 180, by = 30),
  lat = -45
)

result <- utmups_fwd(lon_transect)
data.frame(
  lon = result$lon,
  zone = result$zone,
  northp = result$northp,
  crs = result$crs
)

## ----utm-polar----------------------------------------------------------------
# Antarctic stations
utmups_fwd(antarctic)

## ----utm-roundtrip------------------------------------------------------------
fwd <- utmups_fwd(australia)
rev <- utmups_rev(fwd$x, fwd$y, fwd$zone, fwd$northp)

# Verify accuracy
max(abs(rev$lon - australia[,1]))
max(abs(rev$lat - australia[,2]))

## ----tm-custom----------------------------------------------------------------
# Custom TM centered on Tasmania
tm_fwd(australia, lon0 = 147, k0 = 1.0)

# Compare series approximation vs exact
pts <- cbind(lon = c(147, 148, 149), lat = c(-42, -43, -44))
tm_fwd(pts, lon0 = 147)       # Fast (~5nm accuracy)
tm_exact_fwd(pts, lon0 = 147) # Exact (slower)

## ----lcc-single---------------------------------------------------------------
# Project Australia using a single standard parallel at -35
lcc_fwd(australia, lon0 = 135, stdlat = -35)

## ----lcc-two------------------------------------------------------------------
# Project Australia using two standard parallels
result <- lcc_fwd(australia, lon0 = 135, stdlat1 = -18, stdlat2 = -36)
result

## ----lcc-antarctic------------------------------------------------------------
# Antarctic projection centered on the pole
lcc_fwd(antarctic, lon0 = 0, stdlat1 = -71, stdlat2 = -89)

## ----albers-two---------------------------------------------------------------
# Albers Equal Area for Australia
albers_fwd(australia, lon0 = 132, stdlat1 = -18, stdlat2 = -36)

## ----albers-conus-------------------------------------------------------------
# Continental US configuration
conus <- cbind(
  lon = c(-122.42, -74.01, -87.63, -104.99, -118.24),
  lat = c(37.77, 40.71, 41.88, 39.74, 34.05)
)
rownames(conus) <- c("San Francisco", "New York", "Chicago", "Denver", "Los Angeles")

albers_fwd(conus, lon0 = -96, stdlat1 = 29.5, stdlat2 = 45.5)

## ----albers-antarctic---------------------------------------------------------
# Antarctic equal-area projection
albers_fwd(antarctic, lon0 = 0, stdlat1 = -72, stdlat2 = -60)

## ----albers-vs-lcc------------------------------------------------------------
# Compare Albers (equal-area) vs LCC (conformal)
albers_result <- albers_fwd(australia, lon0 = 132, stdlat1 = -18, stdlat2 = -36)
lcc_result <- lcc_fwd(australia, lon0 = 132, stdlat1 = -18, stdlat2 = -36)

data.frame(
  city = rownames(australia),
  albers_scale = round(albers_result$scale, 4),
  lcc_scale = round(lcc_result$scale, 4)
)

## ----polarstereo-antarctic----------------------------------------------------
# Antarctic stations with UPS-standard scale
polarstereo_fwd(antarctic, northp = FALSE, k0 = 0.994)

## ----polarstereo-arctic-------------------------------------------------------
# Arctic circle of points
arctic <- cbind(lon = seq(0, 315, by = 45), lat = 85)
polarstereo_fwd(arctic, northp = TRUE)

# All points at same latitude have same distance from pole
result <- polarstereo_fwd(arctic, northp = TRUE)
sqrt(result$x^2 + result$y^2)  # All equal

## ----polarstereo-pole---------------------------------------------------------
# South pole
polarstereo_fwd(c(0, -90), northp = FALSE)

# North pole
polarstereo_fwd(c(0, 90), northp = TRUE)

## ----azeq-sydney--------------------------------------------------------------
# Project world cities relative to Sydney
sydney <- c(151.21, -33.87)
result <- azeq_fwd(world_pts, lon0 = sydney[1], lat0 = sydney[2])
result

# Distance from Sydney (in km) = sqrt(x^2 + y^2) / 1000
distances <- sqrt(result$x^2 + result$y^2) / 1000
data.frame(
  city = rownames(world_pts),
  distance_km = round(distances)
)

## ----azeq-pole----------------------------------------------------------------
# Distance from South Pole to Antarctic stations
result <- azeq_fwd(antarctic, lon0 = 0, lat0 = -90)

distances <- sqrt(result$x^2 + result$y^2) / 1000
data.frame(
  station = rownames(antarctic),
  lat = antarctic[,2],
  distance_from_pole_km = round(distances)
)

## ----cassini-basic------------------------------------------------------------
# Tasmania centered on Hobart
tasmania <- cbind(
  lon = c(147.32, 145.49, 146.82, 148.29, 147.13),
  lat = c(-42.88, -40.83, -41.44, -42.15, -43.21)
)
rownames(tasmania) <- c("Hobart", "Launceston", "Devonport", "St Helens", "Dover")

cassini_fwd(tasmania, lon0 = 147, lat0 = -42)

## ----cassini-antarctic--------------------------------------------------------
# McMurdo area survey
mcmurdo_area <- cbind(
  lon = c(166.67, 166.40, 167.00, 166.87, 168.40),
  lat = c(-77.85, -77.55, -78.15, -78.65, -77.18)
)
rownames(mcmurdo_area) <- c("McMurdo", "Marble Point", "Black Island",
                            "Minna Bluff", "Cape Adare")

cassini_fwd(mcmurdo_area, lon0 = 166.67, lat0 = -77.85)

## ----gnomonic-routes----------------------------------------------------------
# Project Sydney-London great circle path
sydney_london <- geodesic_path(c(151.21, -33.87), c(-0.13, 51.51), n = 10)

# Project onto gnomonic centered between them
gnomonic_fwd(cbind(sydney_london$lon, sydney_london$lat),
             lon0 = 75, lat0 = 10)

## ----gnomonic-planning--------------------------------------------------------
# Flights from Sydney - project candidate destinations
destinations <- cbind(
  lon = c(-0.13, -74.01, 139.69, 77.22, -43.17),
  lat = c(51.51, 40.71, 35.69, 28.61, -22.91)
)
rownames(destinations) <- c("London", "New York", "Tokyo", "Delhi", "Rio")

# Gnomonic from Sydney shows great circle routes as straight lines
gnomonic_fwd(destinations, lon0 = 151.21, lat0 = -33.87)

## ----osgb---------------------------------------------------------------------
# British locations (using approximate OSGB36 coordinates)
britain <- cbind(
  lon = c(-0.127, -3.188, -4.251, -1.890, -2.587),
  lat = c(51.507, 55.953, 55.864, 52.486, 51.454)
)
rownames(britain) <- c("London", "Edinburgh", "Glasgow", "Birmingham", "Cardiff")

# Convert to OSGB grid
osgb_fwd(britain)

## ----osgb-gridref-------------------------------------------------------------
# Get alphanumeric grid references
osgb_gridref(britain, precision = 3)  # 100m precision

# Parse a grid reference
osgb_gridref_rev("TQ308080")