[](https://github.com/hypertidy/geodist/actions?query=workflow%3AR-CMD-check)
[](https://github.com/hypertidy/geodist/actions?query=workflow%3Apkgcheck)
[](https://www.repostatus.org/)
[](https://app.codecov.io/gh/hypertidy/geodist)
[](https://cran.r-project.org/package=geodist/)
# 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`).
``` r
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:
``` r
install.packages ("geodist") # current CRAN version
```
Alternatively, current development versions can be installed using any
of the following options:
``` r
# 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
``` r
library (geodist)
packageVersion ("geodist")
#> [1] '0.1.1'
```
## Detailed Usage
Input(s) to the `geodist()` function can be in arbitrary rectangular
format.
``` r
n <- 1e1
x <- tibble::tibble (
x = -180 + 360 * runif (n),
y = -90 + 180 * runif (n)
)
dim (geodist (x, measure = "haversine"))
#> [1] 10 10
y <- tibble::tibble (
x = -180 + 360 * runif (2 * n),
y = -90 + 180 * runif (2 * n)
)
dim (geodist (x, y, measure = "haversine"))
#> [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, measure = "haversine"))
#> [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](https://github.com/mapbox/cheap-ruler-cpp/blob/master/include/mapbox/cheap_ruler.hpp),
and the “reference” implementation of [Karney
(2013)](https://link.springer.com/content/pdf/10.1007/s00190-012-0578-z.pdf),
as implemented in the package
[`sf`](https://cran.r-project.org/package=sf). (Note that `geodist` does
not accept [`sf`](https://cran.r-project.org/package=sf)-format objects;
the [`sf`](https://cran.r-project.org/package=sf) package itself should
be used for that.) The [mapbox cheap ruler
algorithm](https://github.com/mapbox/cheap-ruler-cpp) 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)](https://link.springer.com/content/pdf/10.1007/s00190-012-0578-z.pdf).
``` r
geodist_benchmark (lat = 30, d = 1000)
#> haversine vincenty cheap
#> absolute 0.728786064 0.728786064 0.552930383
#> relative 0.002009509 0.002009509 0.001564216
```
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.
``` r
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",
cheap <- geodist (x, measure = "cheap"),
haversine <- geodist (x, measure = "haversine"),
vincenty <- geodist (x, measure = "vincenty"),
geodesic <- geodist (x, measure = "geodesic")
) [, 1:4]
#> test replications elapsed relative
#> 1 cheap <- geodist(x, measure = "cheap") 10 0.036 1.000
#> 4 geodesic <- geodist(x, measure = "geodesic") 10 1.397 38.806
#> 2 haversine <- geodist(x, measure = "haversine") 10 0.048 1.333
#> 3 vincenty <- geodist(x, measure = "vincenty") 10 0.085 2.361
```
Geodesic distance calculation is available in the [`sf`
package](https://cran.r-project.org/package=sf). Comparing computation
speeds requires conversion of sets of numeric lon-lat points to `sf`
form with the following code:
``` r
x_to_sf <- function (x) {
sapply (seq_len (nrow (x)), function (i) {
sf::st_point (x [i, ]) |>
sf::st_sfc ()
}) |>
sf::st_sfc (crs = 4326)
}
```
Distances in `sf` are by default calculated via
[`s2::s2_distance()`](https://r-spatial.github.io/s2/reference/s2_is_collection.html),
with alternative calculations using the same geodesic algorithm as here.
``` r
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, s2 = TRUE) {
if (s2) {
sf::sf_use_s2 (TRUE) # s2::s2_distance()
} else {
sf::sf_use_s2 (FALSE) # Karney's geodesic algorithm
}
sf::st_distance (xsf, xsf)
}
geo_dist <- function (x) geodist (x, measure = "geodesic")
rbenchmark::benchmark (
replications = 10, order = "test",
sf_dist (xsf, s2 = TRUE),
sf_dist (xsf, s2 = FALSE),
geo_dist (x)
) [, 1:4]
#> Spherical geometry (s2) switched off
#> Linking to GEOS 3.14.1, GDAL 3.13.0, PROJ 9.8.1; sf_use_s2() is FALSE
#> Spherical geometry (s2) switched on
#> Spherical geometry (s2) switched off
#> test replications elapsed relative
#> 3 geo_dist(x) 10 0.028 1.000
#> 2 sf_dist(xsf, s2 = FALSE) 10 0.096 3.429
#> 1 sf_dist(xsf, s2 = TRUE) 10 0.069 2.464
```
The [`geosphere` package](https://cran.r-project.org/package=geosphere)
also offers sequential calculation which is benchmarked with the
following code:
``` r
fgeodist <- function () geodist (x, measure = "vincenty", sequential = TRUE)
fgeosph <- function () geosphere::distVincentySphere (x)
rbenchmark::benchmark (
replications = 10, order = "test",
fgeodist (),
fgeosph ()
) [, 1:4]
#> test replications elapsed relative
#> 1 fgeodist() 10 0.007 1.000
#> 2 fgeosph() 10 0.017 2.429
```
`geodist` is thus at least twice as fast as both `sf` for highly
accurate geodesic distance calculations, and `geosphere` for calculation
of sequential distances.
## Contributors
All contributions to this project are gratefully acknowledged using the
[`allcontributors` package](https://github.com/ropensci/allcontributors)
following the [allcontributors](https://allcontributors.org)
specification. Contributions of any kind are welcome!
### Code
### Issue Authors
### Issue Contributors
