The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

Introduction to geographiclib

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)))
#> [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.034   0.000   0.035

system.time(geohash_fwd(x))
#>    user  system elapsed 
#>   0.008   0.005   0.012

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

Further Reading

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.