knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
library(geographiclib)
GeographicLib C++ library for precise geodetic calculations. All functions are fully vectorized and return rich metadata.
Convert geographic coordinates to MGRS codes and back:
# 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):
mgrs_fwd(c(147.325, -42.881), precision = 0:5)
Direct access to UTM projections with automatic zone selection:
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:
utm <- utmups_fwd(c(147.325, -42.881)) utmups_rev(utm$x, utm$y, utm$zone, utm$northp)
Encode locations as compact strings with hierarchical precision:
# 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:
geohash_resolution(c(4, 6, 8, 10, 12))
Military grid system with 3 precision levels:
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")
Aviation-oriented grid system:
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")
The LCC projection is widely used for aeronautical charts and regional coordinate systems:
# 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)
Project relative to any center point - preserves distances from center:
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
Historical projection used for large-scale topographic mapping:
pts <- cbind(lon = c(-100, -99, -101), lat = c(40, 41, 39)) cassini_fwd(pts, lon0 = -100, lat0 = 40)
Geodesics appear as straight lines - useful for great circle route planning:
# Project relative to a center point gnomonic_fwd(cbind(lon = c(-74, 2.3), lat = c(40.7, 48.9)), lon0 = -37, lat0 = 46)
Grid reference system for Great Britain (requires OSGB36 datum coordinates):
# 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")
Convert to East-North-Up coordinates relative to a local origin:
# 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)
Convert between geodetic and Earth-Centered Earth-Fixed coordinates:
geocentric_fwd(c(-0.1, 51.5)) # With height geocentric_fwd(c(-0.1, 51.5), h = 1000)
Access ellipsoid parameters and derived quantities:
ellipsoid_params() ellipsoid_curvature(c(0, 45, 90))
Solve geodesic problems on the WGS84 ellipsoid.
# London to New York geodesic_inverse(c(-0.1, 51.5), c(-74, 40.7))
# 1000km east from London geodesic_direct(c(-0.1, 51.5), azi = 90, s = 1000000)
Generate points along the shortest path between two locations:
path <- geodesic_path(c(-0.1, 51.5), c(-74, 40.7), n = 5) path
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 are paths of constant bearing. They are not the shortest path (geodesics are), but they maintain a constant compass heading - useful for navigation.
# 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:
# Heading due east maintains latitude rhumb_direct(c(0, 45), azi = 90, s = 1000000)
Compute accurate polygon area and perimeter on the ellipsoid:
# 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:
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))
All functions handle polar regions automatically using UPS (Universal Polar Stereographic) when appropriate:
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)
All functions are implemented in C++ and fully vectorized:
# 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,]))
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.