Introduction to geographiclib

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
library(geographiclib)

Contents

MGRS - Military Grid Reference System package provides R bindings to the

GeographicLib 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:

# 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)

UTM/UPS - Universal Transverse Mercator

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)

Geohash

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))

GARS - Global Area Reference System

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")

Georef - World Geographic Reference System

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")

Lambert Conformal Conic Projection

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)

Azimuthal Equidistant Projection

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

Cassini-Soldner Projection

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)

Gnomonic Projection

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)

OSGB - Ordnance Survey National Grid

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")

Local Cartesian (ENU) Coordinates

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)

Geocentric (ECEF) Coordinates

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)

WGS84 Ellipsoid

Access ellipsoid parameters and derived quantities:

ellipsoid_params()

ellipsoid_curvature(c(0, 45, 90))

Geodesic Calculations

Solve geodesic problems on the WGS84 ellipsoid.

Inverse problem: distance and bearing between points

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

Direct problem: destination from start, bearing, and distance

# 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:

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

Distance 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.

# 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)

Polygon Area

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))

Polar Regions

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)

Performance

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,]))

Further Reading



Try the geographiclib package in your browser

Any scripts or data that you put into this service are public.

geographiclib documentation built on March 4, 2026, 9:07 a.m.