--- title: "Introduction to geographiclib" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Introduction to geographiclib} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` ```{r setup} library(geographiclib) ``` ## Contents - [MGRS - Military Grid Reference System](#mgrs---military-grid-reference-system) - [UTM/UPS - Universal Transverse Mercator](#utmups---universal-transverse-mercator) - [Geohash](#geohash) - [GARS - Global Area Reference System](#gars---global-area-reference-system) - [Georef - World Geographic Reference System](#georef---world-geographic-reference-system) - [Lambert Conformal Conic Projection](#lambert-conformal-conic-projection) - [Azimuthal Equidistant Projection](#azimuthal-equidistant-projection) - [Cassini-Soldner Projection](#cassini-soldner-projection) - [Gnomonic Projection](#gnomonic-projection) - [OSGB - Ordnance Survey National Grid](#osgb---ordnance-survey-national-grid) - [Local Cartesian (ENU) Coordinates](#local-cartesian-enu-coordinates) - [Geocentric (ECEF) Coordinates](#geocentric-ecef-coordinates) - [WGS84 Ellipsoid](#wgs84-ellipsoid) - [Geodesic Calculations](#geodesic-calculations) - [Rhumb Lines (Loxodromes)](#rhumb-lines-loxodromes) - [Polygon Area](#polygon-area) - [Polar Regions](#polar-regions) - [Performance](#performance) - [Further Reading](#further-reading) ## MGRS - Military Grid Reference System package provides R bindings to the [GeographicLib](https://geographiclib.sourceforge.io/) C++ library for precise geodetic calculations. All functions are fully vectorized and return rich metadata. ## MGRS - Military Grid Reference System Convert geographic coordinates to MGRS codes and back: ```{r mgrs} # Forward conversion (code <- mgrs_fwd(c(147.325, -42.881))) # Reverse returns rich metadata mgrs_rev(code) ``` Variable precision from 100km (0) to 1m (5): ```{r mgrs-precision} mgrs_fwd(c(147.325, -42.881), precision = 0:5) ``` ## UTM/UPS - Universal Transverse Mercator Direct access to UTM projections with automatic zone selection: ```{r utmups} pts <- cbind(lon = c(147.325, -74.006, 0), lat = c(-42.881, 40.7128, 88)) utmups_fwd(pts) ``` Reverse conversion requires zone and hemisphere: ```{r utmups-rev} utm <- utmups_fwd(c(147.325, -42.881)) utmups_rev(utm$x, utm$y, utm$zone, utm$northp) ``` ## Geohash Encode locations as compact strings with hierarchical precision: ```{r 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) ``` Resolution table for different lengths: ```{r geohash-res} geohash_resolution(c(4, 6, 8, 10, 12)) ``` ## GARS - Global Area Reference System Military grid system with 3 precision levels: ```{r 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 - World Geographic Reference System Aviation-oriented grid system: ```{r 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") ``` ## Lambert Conformal Conic Projection The LCC projection is widely used for aeronautical charts and regional coordinate systems: ```{r 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) ``` ## Azimuthal Equidistant Projection Project relative to any center point - preserves distances from center: ```{r 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-Soldner Projection Historical projection used for large-scale topographic mapping: ```{r cassini} pts <- cbind(lon = c(-100, -99, -101), lat = c(40, 41, 39)) cassini_fwd(pts, lon0 = -100, lat0 = 40) ``` ## Gnomonic Projection Geodesics appear as straight lines - useful for great circle route planning: ```{r 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 - Ordnance Survey National Grid Grid reference system for Great Britain (requires OSGB36 datum coordinates): ```{r 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") ``` ## Local Cartesian (ENU) Coordinates Convert to East-North-Up coordinates relative to a local origin: ```{r 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 (ECEF) Coordinates Convert between geodetic and Earth-Centered Earth-Fixed coordinates: ```{r geocentric} geocentric_fwd(c(-0.1, 51.5)) # With height geocentric_fwd(c(-0.1, 51.5), h = 1000) ``` ## WGS84 Ellipsoid Access ellipsoid parameters and derived quantities: ```{r ellipsoid} ellipsoid_params() ellipsoid_curvature(c(0, 45, 90)) ``` ## Geodesic Calculations Solve geodesic problems on the WGS84 ellipsoid. ### Inverse problem: distance and bearing between points ```{r geodesic-inverse} # London to New York geodesic_inverse(c(-0.1, 51.5), c(-74, 40.7)) ``` ### Direct problem: destination from start, bearing, and distance ```{r geodesic-direct} # 1000km east from London geodesic_direct(c(-0.1, 51.5), azi = 90, s = 1000000) ``` ### Geodesic paths Generate points along the shortest path between two locations: ```{r geodesic-path} path <- geodesic_path(c(-0.1, 51.5), c(-74, 40.7), n = 5) path ``` ### Distance matrix ```{r 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 Lines (Loxodromes) Rhumb lines are paths of constant bearing. They are not the shortest path (geodesics are), but they maintain a constant compass heading - useful for navigation. ```{r 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) ``` Key properties: east-west rhumb lines maintain constant latitude, north-south rhumb lines maintain constant longitude: ```{r rhumb-ew} # Heading due east maintains latitude rhumb_direct(c(0, 45), azi = 90, s = 1000000) ``` ## Polygon Area Compute accurate polygon area and perimeter on the ellipsoid: ```{r 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 ``` Multiple polygons with the `id` argument: ```{r 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 Regions All functions handle polar regions automatically using UPS (Universal Polar Stereographic) when appropriate: ```{r 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 All functions are implemented in C++ and fully vectorized: ```{r 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,])) ``` ## Further Reading - [GeographicLib documentation](https://geographiclib.sourceforge.io/) - [MGRS on Wikipedia](https://en.wikipedia.org/wiki/Military_Grid_Reference_System) - [Geohash on Wikipedia](https://en.wikipedia.org/wiki/Geohash) - [Geodesics on an ellipsoid](https://en.wikipedia.org/wiki/Geodesics_on_an_ellipsoid)