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 rectagular 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 geodist from github with:
# install.packages("devtools")
devtools::install_github("hypertidy/geodist")
library(geodist)
# current verison
packageVersion("geodist")
#> [1] '0.0.1'
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))
#> [1] 10 10
y <- tibble::tibble (x = -180 + 360 * runif (2 * n),
y = -90 + 180 * runif (2 * n))
dim (geodist (x, y))
#> [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))
#> [1] 10 10
Distances 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)
#> haversine vincenty cheap
#> absolute 0.821979685 0.821979685 0.573589772
#> relative 0.002206036 0.002206036 0.001613667
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%.
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.158 1.000
#> 2 d2 <- geodist(x, measure = "haversine") 10 0.243 1.538
#> 3 d3 <- geodist(x, measure = "vincenty") 10 0.394 2.494
#> 4 d4 <- geodist(x, measure = "geodesic") 10 4.559 28.854
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)
#> Loading required package: 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]
#> Linking to GEOS 3.5.0, GDAL 2.1.3, proj.4 4.9.2
#> test replications elapsed relative
#> 2 geo_dist(x) 10 0.080 1.000
#> 1 sf_dist(xsf) 10 0.353 4.412
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] "7.4506e-09"
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]
#> test replications elapsed relative
#> 1 fgeodist() 10 0.024 1.000
#> 2 fgeosph() 10 0.100 4.167
geodist
is thus around 4 times faster than both sf
for highly accurate
geodesic distance calculations, and geosphere
for calculation of sequential
distances.
require (devtools)
require (testthat)
date()
#> [1] "Mon Jul 2 11:21:57 2018"
devtools::test("tests/")
#> Loading geodist
#> Testing geodist
#> ✔ | OK F W S | Context
✔ | 17 | misc tests [0.1 s]
✔ | 12 | geodist input formats
✔ | 18 | geodist measures
#>
#> ══ Results ════════════════════════════════════════════════════════════════
#> Duration: 0.3 s
#>
#> OK: 47
#> Failed: 0
#> Warnings: 0
#> Skipped: 0