Introduction to geographiclib

Source: vignettes/geographiclib-overview.Rmd

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)))
#> [1] "55GEN2654152348"

# Reverse returns rich metadata
mgrs_rev(code)
#>       lon     lat        x       y zone northp precision convergence     scale
#> 1 147.325 -42.881 526541.5 5252348   55  FALSE         5  -0.2211584 0.9996087
#>   grid_zone square_100km        crs
#> 1       55G           EN EPSG:32755

Variable precision from 100km (0) to 1m (5):

mgrs_fwd(c(147.325, -42.881), precision = 0:5)
#> [1] "55GEN"

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)
#>           x       y zone northp convergence     scale     lon      lat
#> 1  526541.3 5252349   55  FALSE  -0.2211566 0.9996087 147.325 -42.8810
#> 2  583959.4 4507351   18   TRUE   0.6483920 0.9996868 -74.006  40.7128
#> 3 2000000.0 1777931    0   TRUE   0.0000000 0.9943028   0.000  88.0000
#>          crs
#> 1 EPSG:32755
#> 2 EPSG:32618
#> 3 EPSG:32661

Reverse conversion requires zone and hemisphere:

utm <- utmups_fwd(c(147.325, -42.881))
utmups_rev(utm$x, utm$y, utm$zone, utm$northp)
#>       lon     lat        x       y zone northp convergence     scale        crs
#> 1 147.325 -42.881 526541.3 5252349   55  FALSE  -0.2211566 0.9996087 EPSG:32755

Geohash

Encode locations as compact strings with hierarchical precision:

# Default length 12 (~19mm precision)
(gh <- geohash_fwd(c(147.325, -42.881)))
#> [1] "r22u03yb164p"

# Truncation preserves containment
substr(gh, 1, c(4, 6, 8))
#> [1] "r22u"

# Reverse shows resolution
geohash_rev(gh)
#>       lon     lat len lat_resolution lon_resolution
#> 1 147.325 -42.881  12   1.676381e-07   3.352761e-07

Resolution table for different lengths:

geohash_resolution(c(4, 6, 8, 10, 12))
#>   len lat_resolution lon_resolution
#> 1   4   1.757812e-01   3.515625e-01
#> 2   6   5.493164e-03   1.098633e-02
#> 3   8   1.716614e-04   3.433228e-04
#> 4  10   5.364418e-06   1.072884e-05
#> 5  12   1.676381e-07   3.352761e-07

GARS - Global Area Reference System

Military grid system with 3 precision levels:

gars_fwd(c(-74, 40.7), precision = 0)
#> [1] "213LX"
gars_fwd(c(-74, 40.7), precision = 1)
#> [1] "213LX3"
gars_fwd(c(-74, 40.7), precision = 2)
#> [1] "213LX31"

gars_rev("213LR29")
#>         lon      lat precision lat_resolution lon_resolution
#> 1 -73.54167 37.79167         2     0.08333333     0.08333333

Georef - World Geographic Reference System

Aviation-oriented grid system:

georef_fwd(c(-0.1, 51.5), precision = 0)
#> [1] "MKQG"
georef_fwd(c(-0.1, 51.5), precision = 1)
#> [1] "MKQG5430"
georef_fwd(c(-0.1, 51.5), precision = 2)
#> [1] "MKQG5430"

georef_rev("GJPJ3230")
#>         lon      lat precision lat_resolution lon_resolution
#> 1 -76.45833 38.50833         2    0.001666667    0.001666667

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)
#>           x        y convergence    scale  lon lat
#> 1      0.00      0.0   0.0000000 1.000000 -100  40
#> 2  84146.25 111521.9   0.6427876 1.000152  -99  41
#> 3 165789.24 224013.2   1.2855752 1.000613  -98  42

# Two standard parallels (secant cone)
lcc_fwd(pts, lon0 = -96, stdlat1 = 33, stdlat2 = 45)
#>           x        y convergence     scale  lon lat
#> 1 -339643.8 108321.8   -2.521985 0.9946660 -100  40
#> 2 -251122.5 215464.1   -1.891489 0.9950973  -99  41
#> 3 -164998.9 323691.8   -1.260993 0.9958400  -98  42

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
#>          x       y       azi     scale   lon   lat lon0 lat0
#> 1 -5302905 1761511 -128.7635 0.8770643 -74.0  40.7 -0.1 51.5
#> 2  5029432 8155050  156.2494 0.6652729 139.7  35.7 -0.1 51.5
#> 3 14783025 8374075  139.2137 0.1711044 151.2 -33.9 -0.1 51.5

# Distance from London = sqrt(x^2 + y^2)
sqrt(result$x^2 + result$y^2) / 1000
#> [1]  5587.820  9581.233 16990.084

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)
#>           x             y      azi        rk  lon lat
#> 1      0.00  7.069289e-10 90.00000 1.0000000 -100  40
#> 2  84133.35  1.115260e+05 90.65610 0.9999129  -99  41
#> 3 -86624.66 -1.105493e+05 89.37064 0.9999076 -101  39

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)
#>          x         y        azi        rk   lon  lat
#> 1 -3276316  127487.1 -113.67126 0.8893703 -74.0 40.7
#> 2  2972069 1123467.3   98.77886 0.8951499   2.3 48.9

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)
#>    easting northing convergence     scale    lon    lat
#> 1 529972.5 180390.9    1.466171 0.9998087 -0.127 51.507

# Get alphanumeric grid reference
osgb_gridref(london, precision = 3)  # 100m precision
#> [1] "TQ299803"

# Parse a grid reference
osgb_gridref_rev("TQ308080")
#>          lon      lat easting northing precision
#> 1 -0.1407134 50.85658  530850   108050         3

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)
#>           x         y         z  lon  lat h
#> 1     0.000      0.00   0.00000 -0.1 51.5 0
#> 2 -6928.841  11130.61 -13.47326 -0.2 51.6 0
#> 3  6959.234 -11120.93 -13.48954  0.0 51.4 0

Geocentric (ECEF) Coordinates

Convert between geodetic and Earth-Centered Earth-Fixed coordinates:

geocentric_fwd(c(-0.1, 51.5))
#>         X         Y       Z  lon  lat h
#> 1 3978642 -6944.048 4968362 -0.1 51.5 0

# With height
geocentric_fwd(c(-0.1, 51.5), h = 1000)
#>         X         Y       Z  lon  lat    h
#> 1 3979265 -6945.135 4969145 -0.1 51.5 1000

WGS84 Ellipsoid

Access ellipsoid parameters and derived quantities:

ellipsoid_params()
#> $a
#> [1] 6378137
#> 
#> $f
#> [1] 0.003352811
#> 
#> $b
#> [1] 6356752
#> 
#> $e2
#> [1] 0.00669438
#> 
#> $ep2
#> [1] 0.006739497
#> 
#> $n
#> [1] 0.00167922
#> 
#> $area
#> [1] 5.100656e+14
#> 
#> $volume
#> [1] 1.083207e+21

ellipsoid_curvature(c(0, 45, 90))
#>   lat meridional transverse
#> 1   0    6335439    6378137
#> 2  45    6367382    6388838
#> 3  90    6399594    6399594

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))
#>   lon1 lat1 lon2 lat2     s12      azi1      azi2     m12       M12       M21
#> 1 -0.1 51.5  -74 40.7 5587820 -71.62462 -128.7635 4900877 0.6407216 0.6404073
#>             S12
#> 1 -4.040644e+13

Direct problem: destination from start, bearing, and distance 

# 1000km east from London
geodesic_direct(c(-0.1, 51.5), azi = 90, s = 1000000)
#>   lon1 lat1 azi1   s12     lon2     lat2     azi2    m12       M12       M21
#> 1 -0.1 51.5   90 1e+06 14.12014 50.62607 101.0838 995914 0.9877522 0.9877514
#>            S12
#> 1 7.838198e+12

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
#>         lon      lat        azi       s
#> 1  -0.10000 51.50000  -71.62462       0
#> 2 -20.47028 53.76472  -87.88020 1396955
#> 3 -41.26329 52.38360 -104.56966 2793910
#> 4 -59.44254 47.71610 -118.56584 4190865
#> 5 -74.00000 40.70000 -128.76352 5587820

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)
#>       [,1]  [,2]  [,3]  [,4]
#> [1,]     0  5588  9581 16990
#> [2,]  5588     0 10873 15991
#> [3,]  9581 10873     0  7797
#> [4,] 16990 15991  7797     0

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
#> [1] 5587.82
rhumb_inverse(c(-0.1, 51.5), c(-74, 40.7))$s12 / 1000     # km (longer!)
#> [1] 5812.568

# Generate rhumb line path
rhumb_path(c(-0.1, 51.5), c(-74, 40.7), n = 5)
#> Warning in format.data.frame(if (omit) x[seq_len(n0), , drop = FALSE] else x, :
#> corrupt data frame: columns will be truncated or padded with NAs
#>         lon      lat       s     azi12
#> 1  -0.10000 51.50000       0 -101.9189
#> 2 -20.00072 48.80191 1453142      <NA>
#> 3 -38.86061 46.10255 2906284      <NA>
#> 4 -56.82095 43.40192 4359426      <NA>
#> 5 -74.00000 40.70000 5812568      <NA>

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)
#>   lon1 lat1 azi12   s12     lon2 lat2          S12
#> 1    0   45    90 1e+06 12.68282   45 6.338984e+12

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
#> [1] -1.444896e+14
#> 
#> $perimeter
#> [1] 38568531
#> 
#> $n
#> [1] 3

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

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))
#>   id       area perimeter n
#> 1  1 6154854787  378793.4 3
#> 2  2 6154854787  378793.4 3

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)
#> [1] "ZAF0000088973" "BAL0000088973"

# UTM/UPS shows zone 0 for polar
utmups_fwd(polar)
#>       x       y zone northp convergence     scale lon lat        crs
#> 1 2e+06 1888973    0   TRUE           0 0.9940757   0  89 EPSG:32661
#> 2 2e+06 1888973    0  FALSE        -180 0.9940757 180 -89 EPSG:32761

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))
#>    user  system elapsed 
#>   0.035   0.000   0.035

system.time(geohash_fwd(x))
#>    user  system elapsed 
#>   0.013   0.000   0.012

system.time(geodesic_distance_matrix(x[1:100,]))
#>    user  system elapsed 
#>   0.059   0.000   0.060

Further Reading