Folders and files

Name		Name	Last commit message	Last commit date
Latest commit History 258 Commits
.github/workflows		.github/workflows
.hooks		.hooks
R		R
man		man
src		src
tests		tests
vignettes		vignettes
.Rbuildignore		.Rbuildignore
.gitignore		.gitignore
.pre-commit-config.yaml		.pre-commit-config.yaml
CONDUCT.md		CONDUCT.md
CONTRIBUTING.md		CONTRIBUTING.md
DESCRIPTION		DESCRIPTION
LICENSE		LICENSE
NAMESPACE		NAMESPACE
NEWS.md		NEWS.md
README.Rmd		README.Rmd
README.md		README.md
codemeta.json		codemeta.json
cran-comments.md		cran-comments.md
makefile		makefile

Repository files navigation

geodist

An ultra-lightweight, zero-dependency package for very fast calculationof geodesic distances. Main eponymous function,geodist(), acceptsonly one or two primary arguments, which must be rectangular objectswith unambiguously labelled longitude and latitude columns (that is,some variant oflon/lat, orx/y).

n<-50x<- 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 matrixd1<- geodist (x,y)# A 50-by-100 matrixd2<- geodist (x,sequential=TRUE)# Vector of length 49d2<- geodist (x,sequential=TRUE,pad=TRUE)# Vector of length 50

Installation

You can install latest stable version ofgeodist from CRAN with:

install.packages ("geodist")# current CRAN version

Alternatively, current development versions can be installed using anyof 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.1.0.6'

Detailed Usage

Input(s) to thegeodist() function can be in arbitrary rectangularformat.

n<-1e1x<-tibble::tibble (x=-180+360* runif (n),y=-90+180* runif (n))dim (geodist (x,measure="haversine"))#> [1] 10 10y<-tibble::tibble (x=-180+360* runif (2*n),y=-90+180* runif (2*n))dim (geodist (x,y,measure="haversine"))#> [1] 10 20x<- 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 ofspherical and elliptical distance measures. Distance measures currentlyimplemented are Haversine, Vincenty (spherical and elliptical)), thevery fastmapbox cheapruler(see theirblogpost),and the “reference” implementation ofKarney(2013),as implemented in the packagesf. (Note thatgeodist doesnot acceptsf-format objects;thesf package itself shouldbe used for that.) Themapbox cheap ruleralgorithm is intended toprovide approximate yet very fast distance calculations within smallareas (tens to a few hundred kilometres across).

Benchmarks of geodesic accuracy

Thegeodist_benchmark() function - the only other function provided bythegeodist package - compares the accuracy of the different metricsto the nanometre-accuracy standard ofKarney(2013).

geodist_benchmark (lat=30,d=1000)#>            haversine    vincenty      cheap#> absolute 0.818681120 0.818681120 0.57910408#> relative 0.002173428 0.002173428 0.00160801

All distances (d) are in metres, and all measures are accurate towithin 1m over distances out to several km (at the chosen latitude of 30degrees). The following plots compare the absolute and relativeaccuracies of the different distance measures implemented here. Themapbox cheap ruler algorithm is the most accurate for distances out toaround 100km, beyond which it becomes extremely inaccurate. Averagerelative errors of Vincenty distances remain generally constant ataround 0.2%, while relative errors of cheap-ruler distances out to 100kmare around 0.16%.

Performance comparison

The following code demonstrates the relative speed advantages of thedifferent distance measures implemented in thegeodist package.

n<-1e3dx<-dy<-0.01x<- 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.070      1.0#> 2 d2 <- geodist(x, measure = "haversine")           10   0.133      1.9#> 3  d3 <- geodist(x, measure = "vincenty")           10   0.224      3.2#> 4  d4 <- geodist(x, measure = "geodesic")           10   2.856     40.8

Geodesic distance calculation is available in thesfpackage. Comparing computationspeeds requires conversion of sets of numeric lon-lat points tosfform with the following code:

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 insf are by default calculated vias2::s2_distance(),with alternative calculations using the same geodesic algorithm as here.

n<-1e2x<- 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.13.0, GDAL 3.10.0, PROJ 9.5.0; 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.062    1.000#> 2 sf_dist(xsf, s2 = FALSE)           10   0.207    3.339#> 1  sf_dist(xsf, s2 = TRUE)           10   0.146    2.355

Thegeosphere packagealso offers sequential calculation which is benchmarked with thefollowing 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.017    1.000#> 2  fgeosph()           10   0.037    2.176

geodist is thus at least twice as fast as bothsf for highlyaccurate geodesic distance calculations, andgeosphere for calculationof sequential distances.

Contributors

All contributions to this project are gratefully acknowledged using theallcontributors package following theallcontributors specification. Contributions of any kind are welcome!