knitr::opts_chunk$set ( collapse = TRUE, comment = "#>", fig.path = "README-" )
An ultra-lightweight, zero-dependency package for very fast calculation of
geodesic distances. Main eponymous function, geodist()
, accepts only one or
two primary arguments, which must be rectangular objects with unambiguously
labelled longitude and latitude columns (that is, some variant of lon
/lat
,
or x
/y
).
n <- 50 x <- cbind (-10 + 20 * runif (n), -10 + 20 * runif (n)) y <- cbind (-10 + 20 * runif (2 * n), -10 + 20 * runif (2 * n)) colnames (x) <- colnames (y) <- c ("x", "y") d0 <- geodist (x) # A 50-by-50 matrix d1 <- geodist (x, y) # A 50-by-100 matrix d2 <- geodist (x, sequential = TRUE) # Vector of length 49 d2 <- geodist (x, sequential = TRUE, pad = TRUE) # Vector of length 50
You can install latest stable version of geodist
from CRAN with:
install.packages ("geodist") # current CRAN version
Alternatively, current development versions can be installed using any of the following options:
# install.packages("remotes") remotes::install_git ("https://git.sr.ht/~mpadge/geodist") remotes::install_git ("https://codeberg.org/hypertidy/geodist") remotes::install_bitbucket ("hypertidy/geodist") remotes::install_gitlab ("hypertidy/geodist") remotes::install_github ("hypertidy/geodist")
Then load with
library (geodist) packageVersion ("geodist")
Input(s) to the geodist()
function can be in arbitrary rectangular format.
n <- 1e1 x <- tibble::tibble ( x = -180 + 360 * runif (n), y = -90 + 180 * runif (n) ) dim (geodist (x)) y <- tibble::tibble ( x = -180 + 360 * runif (2 * n), y = -90 + 180 * runif (2 * n) ) dim (geodist (x, y)) x <- cbind ( -180 + 360 * runif (n), -90 + 100 * runif (n), seq (n), runif (n) ) colnames (x) <- c ("lon", "lat", "a", "b") dim (geodist (x))
All outputs are distances in metres, calculated with a variety of spherical and
elliptical distance measures. Distance measures currently implemented are
Haversine, Vincenty (spherical and elliptical)), the very fast mapbox cheap
ruler
(see their blog
post),
and the "reference" implementation of Karney
(2013), as
implemented in the package sf
. (Note
that geodist
does not accept
sf
-format objects; the
sf
package itself should be used for
that.) The mapbox cheap ruler
algorithm is intended to provide
approximate yet very fast distance calculations within small areas (tens to a few
hundred kilometres across).
The geodist_benchmark()
function - the only other function provided by the
geodist
package - compares the accuracy of the different metrics to the
nanometre-accuracy standard of Karney
(2013).
geodist_benchmark (lat = 30, d = 1000)
All distances (d)
are in metres, and all measures are accurate to within 1m
over distances out to several km (at the chosen latitude of 30 degrees). The
following plots compare the absolute and relative accuracies of the different
distance measures implemented here. The mapbox cheap ruler algorithm is the
most accurate for distances out to around 100km, beyond which it becomes
extremely inaccurate. Average relative errors of Vincenty distances remain
generally constant at around 0.2%, while relative errors of cheap-ruler
distances out to 100km are around 0.16%.
lat <- 30 d <- 10^(1:35 / 5) # 1m to 100 km y <- lapply (d, function (i) geodist_benchmark (lat = lat, d = i)) # yabs <- do.call (rbind, lapply (y, function (i) i [1, ])) [, c (1, 2, 4)] # yrel <- 100 * do.call (rbind, lapply (y, function (i) i [2, ])) [, c (1, 2, 4)] yabs <- do.call (rbind, lapply (y, function (i) i [1, ])) yrel <- 100 * do.call (rbind, lapply (y, function (i) i [2, ])) yvals <- list (yabs, yrel) cols <- c ("skyblue", "lawngreen", "tomato") par (mfrow = c (1, 2)) ylabs <- c ("Absolute error (m)", "Relative error (%)") ylims <- list (range (yvals [[1]]), c (min (yvals [[2]]), 1)) for (i in 1:2) { plot (NULL, NULL, xlim = range (d / 1000), ylim = ylims [[i]], bty = "l", log = "xy", xaxt = "n", yaxt = "n", xlab = "distance (km)", ylab = ylabs [i] ) axis (d / 1000, side = 1, at = c (0.001, 0.1, 10, 1e3, 1e4), labels = c ("0.001", "0.1", "10", "1000", "") ) if (i == 1) { yl <- 10^(-3:5) axis (yvals [[i]], side = 2, at = c (0.001, 0.1, 10, 100, 10000), labels = c ("0.001", "0.1", "10", "100", "1000") ) } else { yl <- c (0.1, 0.2, 0.3, 0.4, 0.5, 1, 2) axis (yvals [[i]], side = 2, at = yl, labels = c ("0.1", "0.2", "0.3", "0.4", "0.5", "1", "2") ) } junk <- sapply (yl, function (j) { lines (range (d / 1000), rep (j, 2), col = "grey", lty = 2 ) }) xl <- 10^(-3:6) junk <- sapply (xl, function (j) { lines (rep (j, 2), range (yvals [[i]]), col = "grey", lty = 2 ) }) for (j in (1:ncol (yabs)) [-1]) { lines (d / 1000, yvals [[i]] [, j], col = cols [j]) } legend ("topleft", lwd = 1, col = cols, bty = "n", legend = colnames (yvals [[i]]) ) }
The following code demonstrates the relative speed advantages of the different
distance measures implemented in the geodist
package.
n <- 1e3 dx <- dy <- 0.01 x <- cbind (-100 + dx * runif (n), 20 + dy * runif (n)) y <- cbind (-100 + dx * runif (2 * n), 20 + dy * runif (2 * n)) colnames (x) <- colnames (y) <- c ("x", "y") rbenchmark::benchmark ( replications = 10, order = "test", d1 <- geodist (x, measure = "cheap"), d2 <- geodist (x, measure = "haversine"), d3 <- geodist (x, measure = "vincenty"), d4 <- geodist (x, measure = "geodesic") ) [, 1:4]
Geodesic distance calculation is available in the sf
package. Comparing computation speeds
requires conversion of sets of numeric lon-lat points to sf
form with the
following code:
require (magrittr) x_to_sf <- function (x) { sapply (seq (nrow (x)), function (i) { sf::st_point (x [i, ]) %>% sf::st_sfc () }) %>% sf::st_sfc (crs = 4326) }
n <- 1e2 x <- cbind (-180 + 360 * runif (n), -90 + 180 * runif (n)) colnames (x) <- c ("x", "y") xsf <- x_to_sf (x) sf_dist <- function (xsf) sf::st_distance (xsf, xsf) geo_dist <- function (x) geodist (x, measure = "geodesic") rbenchmark::benchmark ( replications = 10, order = "test", sf_dist (xsf), geo_dist (x) ) [, 1:4]
Confirm that the two give almost identical results:
ds <- matrix (as.numeric (sf_dist (xsf)), nrow = length (xsf)) dg <- geodist (x, measure = "geodesic") formatC (max (abs (ds - dg)), format = "e")
All results are in metres, so the two differ by only around 10 nanometres.
The geosphere
package also
offers sequential calculation which is benchmarked with the following code:
n <- 1e4 x <- cbind (-180 + 360 * runif (n), -90 + 180 * runif (n)) colnames (x) <- c ("x", "y")
fgeodist <- function () geodist (x, measure = "vincenty", sequential = TRUE) fgeosph <- function () geosphere::distVincentySphere (x) rbenchmark::benchmark ( replications = 10, order = "test", fgeodist (), fgeosph () ) [, 1:4]
geodist
is thus around 3 times faster than sf
for highly accurate geodesic
distance calculations, and around twice as fast as geosphere
for calculation
of sequential distances.
require (devtools) require (testthat)
date () devtools::test ("tests/")
All contributions to this project are gratefully acknowledged using the allcontributors
package following the all-contributors specification. Contributions of any kind are welcome!
mpadge |
daniellemccool |
mdsumner |
edzer |
njtierney |
mkuehn10 |
asardaes |
marcosci |
mem48 |
dcooley |
Robinlovelace |
espinielli |
Maschette |
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.