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

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

## ----mgrs---------------------------------------------------------------------
# Forward conversion
(code <- mgrs_fwd(c(147.325, -42.881)))

# Reverse returns rich metadata
mgrs_rev(code)

## ----mgrs-precision-----------------------------------------------------------
mgrs_fwd(c(147.325, -42.881), precision = 0:5)

## ----utmups-------------------------------------------------------------------
pts <- cbind(lon = c(147.325, -74.006, 0), lat = c(-42.881, 40.7128, 88))
utmups_fwd(pts)

## ----utmups-rev---------------------------------------------------------------
utm <- utmups_fwd(c(147.325, -42.881))
utmups_rev(utm$x, utm$y, utm$zone, utm$northp)

## ----geohash------------------------------------------------------------------
# Default length 12 (~19mm precision)
(gh <- geohash_fwd(c(147.325, -42.881)))

# Truncation preserves containment
substr(gh, 1, c(4, 6, 8))

# Reverse shows resolution
geohash_rev(gh)

## ----geohash-res--------------------------------------------------------------
geohash_resolution(c(4, 6, 8, 10, 12))

## ----gars---------------------------------------------------------------------
gars_fwd(c(-74, 40.7), precision = 0)
gars_fwd(c(-74, 40.7), precision = 1)
gars_fwd(c(-74, 40.7), precision = 2)

gars_rev("213LR29")

## ----georef-------------------------------------------------------------------
georef_fwd(c(-0.1, 51.5), precision = 0)
georef_fwd(c(-0.1, 51.5), precision = 1)
georef_fwd(c(-0.1, 51.5), precision = 2)

georef_rev("GJPJ3230")

## ----lcc----------------------------------------------------------------------
# Single standard parallel (tangent cone)
pts <- cbind(lon = c(-100, -99, -98), lat = c(40, 41, 42))
lcc_fwd(pts, lon0 = -100, stdlat = 40)

# Two standard parallels (secant cone)
lcc_fwd(pts, lon0 = -96, stdlat1 = 33, stdlat2 = 45)

## ----azeq---------------------------------------------------------------------
pts <- cbind(lon = c(-74, 139.7, 151.2), lat = c(40.7, 35.7, -33.9))
result <- azeq_fwd(pts, lon0 = -0.1, lat0 = 51.5)
result

# Distance from London = sqrt(x^2 + y^2)
sqrt(result$x^2 + result$y^2) / 1000

## ----cassini------------------------------------------------------------------
pts <- cbind(lon = c(-100, -99, -101), lat = c(40, 41, 39))
cassini_fwd(pts, lon0 = -100, lat0 = 40)

## ----gnomonic-----------------------------------------------------------------
# Project relative to a center point
gnomonic_fwd(cbind(lon = c(-74, 2.3), lat = c(40.7, 48.9)), lon0 = -37, lat0 = 46)

## ----osgb---------------------------------------------------------------------
# Convert OSGB36 coordinates to grid
london <- c(-0.127, 51.507)
osgb_fwd(london)

# Get alphanumeric grid reference
osgb_gridref(london, precision = 3)  # 100m precision

# Parse a grid reference
osgb_gridref_rev("TQ308080")

## ----localcartesian-----------------------------------------------------------
# Set up local system centered on London
pts <- cbind(lon = c(-0.1, -0.2, 0.0), lat = c(51.5, 51.6, 51.4))
localcartesian_fwd(pts, lon0 = -0.1, lat0 = 51.5)

## ----geocentric---------------------------------------------------------------
geocentric_fwd(c(-0.1, 51.5))

# With height
geocentric_fwd(c(-0.1, 51.5), h = 1000)

## ----ellipsoid----------------------------------------------------------------
ellipsoid_params()

ellipsoid_curvature(c(0, 45, 90))

## ----geodesic-inverse---------------------------------------------------------
# London to New York
geodesic_inverse(c(-0.1, 51.5), c(-74, 40.7))

## ----geodesic-direct----------------------------------------------------------
# 1000km east from London
geodesic_direct(c(-0.1, 51.5), azi = 90, s = 1000000)

## ----geodesic-path------------------------------------------------------------
path <- geodesic_path(c(-0.1, 51.5), c(-74, 40.7), n = 5)
path

## ----geodesic-matrix----------------------------------------------------------
cities <- cbind(
  lon = c(-0.1, -74, 139.7, 151.2),
  lat = c(51.5, 40.7, 35.7, -33.9)
)
rownames(cities) <- c("London", "NYC", "Tokyo", "Sydney")

# Distance matrix in km
round(geodesic_distance_matrix(cities) / 1000)

## ----rhumb--------------------------------------------------------------------
# Rhumb distance vs geodesic: London to New York
geodesic_inverse(c(-0.1, 51.5), c(-74, 40.7))$s12 / 1000  # km
rhumb_inverse(c(-0.1, 51.5), c(-74, 40.7))$s12 / 1000     # km (longer!)

# Generate rhumb line path
rhumb_path(c(-0.1, 51.5), c(-74, 40.7), n = 5)

## ----rhumb-ew-----------------------------------------------------------------
# Heading due east maintains latitude
rhumb_direct(c(0, 45), azi = 90, s = 1000000)

## ----polygon-area-------------------------------------------------------------
# Triangle: London - New York - Sydney
triangle <- cbind(
  lon = c(-0.1, -74, 151.2),
  lat = c(51.5, 40.7, -33.9)
)
result <- polygon_area(triangle)
result

# Area in millions of km²
abs(result$area) / 1e12

## ----polygon-multiple---------------------------------------------------------
pts <- cbind(
  lon = c(0, 1, 1, 10, 11, 11),
  lat = c(0, 0, 1, 0, 0, 1)
)
polygon_area(pts, id = c(1, 1, 1, 2, 2, 2))

## ----polar--------------------------------------------------------------------
polar <- cbind(c(0, 180), c(89, -89))

# MGRS uses polar grid zones
mgrs_fwd(polar)

# UTM/UPS shows zone 0 for polar
utmups_fwd(polar)

## ----performance, eval=TRUE---------------------------------------------------
# 40,000+ points in milliseconds
x <- cbind(runif(40000, -180, 180), runif(40000, -80, 80))

system.time(mgrs_fwd(x))

system.time(geohash_fwd(x))

system.time(geodesic_distance_matrix(x[1:100,]))