Awesome
<!-- README.md is generated from README.Rmd. Please edit that file -->
geodist
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
Installation
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")
#> [1] '0.0.8.6'
Detailed Usage
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))
#> Maximum distance is > 100km. The 'cheap' measure is inaccurate over such
#> large distances, you'd likely be better using a different 'measure'.
#> [1] 10 10
y <- tibble::tibble (
x = -180 + 360 * runif (2 * n),
y = -90 + 180 * runif (2 * n)
)
dim (geodist (x, y))
#> Maximum distance is > 100km. The 'cheap' measure is inaccurate over such
#> large distances, you'd likely be better using a different 'measure'.
#> [1] 10 20
x <- cbind (
-180 + 360 * runif (n),
-90 + 100 * runif (n),
seq (n), runif (n)
)
colnames (x) <- c ("lon", "lat", "a", "b")
dim (geodist (x))
#> Maximum distance is > 100km. The 'cheap' measure is inaccurate over such
#> large distances, you'd likely be better using a different 'measure'.
#> [1] 10 10
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).
Benchmarks of geodesic accuracy
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)
#> haversine vincenty cheap
#> absolute 0.733504364 0.733504364 0.546134467
#> relative 0.002072017 0.002072017 0.001601954
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%.
Performance comparison
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]
#> test replications elapsed relative
#> 1 d1 <- geodist(x, measure = "cheap") 10 0.068 1.000
#> 2 d2 <- geodist(x, measure = "haversine") 10 0.139 2.044
#> 3 d3 <- geodist(x, measure = "vincenty") 10 0.229 3.368
#> 4 d4 <- geodist(x, measure = "geodesic") 10 3.315 48.750
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)
#> Warning in CPL_gdal_init(): GDAL Error 1: libpodofo.so.0.9.8: cannot open
#> shared object file: No such file or directory
#> Warning in CPL_gdal_init(): GDAL Error 1: libpodofo.so.0.9.8: cannot open
#> shared object file: No such file or directory
#> Warning in CPL_gdal_init(): GDAL Error 1: libpodofo.so.0.9.8: cannot open
#> shared object file: No such file or directory
#> Warning in CPL_gdal_init(): GDAL Error 1: libpodofo.so.0.9.8: cannot open
#> shared object file: No such file or directory
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]
#> test replications elapsed relative
#> 2 geo_dist(x) 10 0.070 1.000
#> 1 sf_dist(xsf) 10 0.144 2.057
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")
#> [1] "3.7676e+04"
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:
fgeodist <- function () geodist (x, measure = "vincenty", sequential = TRUE)
fgeosph <- function () geosphere::distVincentySphere (x)
rbenchmark::benchmark (
replications = 10, order = "test",
fgeodist (),
fgeosph ()
) [, 1:4]
#> The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
#> which was just loaded, will retire in October 2023.
#> Please refer to R-spatial evolution reports for details, especially
#> https://r-spatial.org/r/2023/05/15/evolution4.html.
#> It may be desirable to make the sf package available;
#> package maintainers should consider adding sf to Suggests:.
#> The sp package is now running under evolution status 2
#> (status 2 uses the sf package in place of rgdal)
#> test replications elapsed relative
#> 1 fgeodist() 10 0.017 1.000
#> 2 fgeosph() 10 0.032 1.882
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.
Test Results
require (devtools)
require (testthat)
date ()
#> [1] "Mon Aug 21 11:34:48 2023"
devtools::test ("tests/")
#> ℹ Testing geodist
#> Starting 2 test processes
#> ✔ | F W S OK | Context
#> ⠋ [ FAIL 0 | WARN 0 | SKIP 0 | PASS 0 ] Starting up... ⠙ [ FAIL 0 | WARN 0 | SKIP 0 | PASS 0 ] Starting up... ⠹ [ FAIL 0 | WARN 0 | SKIP 0 | PASS 0 ] Starting up... ⠸ [ FAIL 0 | WARN 0 | SKIP 0 | PASS 0 ] Starting up... ⠼ [ FAIL 0 | WARN 0 | SKIP 0 | PASS 0 ] Starting up... ⠴ [ FAIL 0 | WARN 0 | SKIP 0 | PASS 0 ] Starting up... ⠦ [ FAIL 0 | WARN 0 | SKIP 0 | PASS 0 ] Starting up... ✔ | 4 | geodist-min
#> ⠧ [ FAIL 0 | WARN 0 | SKIP 0 | PASS 16 ] @ geodist ⠇ [ FAIL 0 | WARN 0 | SKIP 0 | PASS 16 ] @ geodist ✔ | 52 | geodist
#> ⠏ [ FAIL 0 | WARN 0 | SKIP 0 | PASS 87 ] @ georange ⠋ [ FAIL 0 | WARN 0 | SKIP 0 | PASS 87 ] @ georange ✔ | 37 | georange
#> ⠙ [ FAIL 0 | WARN 0 | SKIP 0 | PASS 97 ] @ input-format ⠹ [ FAIL 0 | WARN 0 | SKIP 0 | PASS 98 ] @ input-format ✔ | 18 | input-format
#> ⠸ [ FAIL 0 | WARN 0 | SKIP 0 | PASS 111 ] @ measures ⠼ [ FAIL 0 | WARN 0 | SKIP 0 | PASS 111 ] @ measures ✔ | 18 | measures
#> ⠴ [ FAIL 0 | WARN 0 | SKIP 0 | PASS 129 ] Starting up...
#> ══ Results ════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════
#> Duration: 1.5 s
#>
#> [ FAIL 0 | WARN 0 | SKIP 0 | PASS 129 ]
Contributors
<!-- ALL-CONTRIBUTORS-LIST:START - Do not remove or modify this section --> <!-- prettier-ignore-start --> <!-- markdownlint-disable -->All contributions to this project are gratefully acknowledged using the
allcontributors
package following the
all-contributors specification.
Contributions of any kind are welcome!