geodist

Ultra lightweight, ultra fast calculation of geo distances

https://github.com/hypertidy/geodist

Repository

Ultra lightweight, ultra fast calculation of geo distances

README.Rmd

---
output: github_document
---



```{r, echo = FALSE}
knitr::opts_chunk$set (
    collapse = TRUE,
    comment = "#>",
    fig.path = "README-"
)
```

[![R build
status](https://github.com/hypertidy/geodist/workflows/R-CMD-check/badge.svg)](https://github.com/hypertidy/geodist/actions?query=workflow%3AR-CMD-check)
[![pkgcheck](https://github.com/hypertidy/geodist/workflows/pkgcheck/badge.svg)](https://github.com/hypertidy/geodist/actions?query=workflow%3Apkgcheck)
[![Project Status: Active – The project has reached a stable, usable state and is being actively developed.](http://www.repostatus.org/badges/latest/active.svg)](https://www.repostatus.org/)
[![codecov](https://codecov.io/gh/hypertidy/geodist/branch/master/graph/badge.svg)](https://app.codecov.io/gh/hypertidy/geodist)
[![CRAN_Status_Badge](http://www.r-pkg.org/badges/version/geodist)](http://cran.r-project.org/web/packages/geodist/) 
![downloads](http://cranlogs.r-pkg.org/badges/grand-total/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, eval = FALSE}
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 cran-installation, eval = FALSE}
install.packages ("geodist") # current CRAN version
```
Alternatively, current development versions can be installed using any of the
following options:
```{r gh-installation, eval = FALSE}
# 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}
library (geodist)
packageVersion ("geodist")
```

## 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"))
y <- tibble::tibble (
    x = -180 + 360 * runif (2 * n),
    y = -90 + 180 * runif (2 * n)
)
dim (geodist (x, y, measure = "haversine"))
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"))
```
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)
(see their [blog
post](https://blog.mapbox.com/fast-geodesic-approximations-with-cheap-ruler-106f229ad016)),
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)
```
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%.


```{r plot, eval = FALSE, echo = FALSE}
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]])
    )
}
```
![](vignettes/fig1.png)


### Performance comparison

The following code demonstrates the relative speed advantages of the different
distance measures implemented in the `geodist` package.
```{r benchmark2}
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](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, message = FALSE}
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 benchmark}
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]
```
The [`geosphere` package](https://cran.r-project.org/package=geosphere) also
offers sequential calculation which is benchmarked with the following code:
```{r, echo = FALSE}
n <- 1e4
x <- cbind (-180 + 360 * runif (n), -90 + 180 * runif (n))
colnames (x) <- c ("x", "y")
```
```{r sequential}
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 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









mpadge






jlacko






mdsumner






daniellemccool






kadyb






olivroy






### Issue Authors









edzer






marcosci






mem48






dcooley






Robinlovelace






espinielli






Maschette






### Issue Contributors









njtierney






mkuehn10






asardaes

Owner

CodeMeta (codemeta.json)

{
  "@context": "https://doi.org/10.5063/schema/codemeta-2.0",
  "@type": "SoftwareSourceCode",
  "identifier": "geodist",
  "description": "Dependency-free, ultra fast calculation of geodesic distances. Includes the reference nanometre-accuracy geodesic distances of Karney (2013) <doi:10.1007/s00190-012-0578-z>, as used by the 'sf' package, as well as Haversine and Vincenty distances. Default distance measure is the \"Mapbox cheap ruler\" which is generally more accurate than Haversine or Vincenty for distances out to a few hundred kilometres, and is considerably faster. The main function accepts one or two inputs in almost any generic rectangular form, and returns either matrices of pairwise distances, or vectors of sequential distances.",
  "name": "geodist: Fast, Dependency-Free Geodesic Distance Calculations",
  "relatedLink": [
    "https://hypertidy.github.io/geodist/",
    "https://CRAN.R-project.org/package=geodist"
  ],
  "codeRepository": "https://github.com/hypertidy/geodist",
  "issueTracker": "https://github.com/hypertidy/geodist/issues",
  "license": "https://spdx.org/licenses/MIT",
  "version": "0.1.1.001",
  "programmingLanguage": {
    "@type": "ComputerLanguage",
    "name": "R",
    "url": "https://r-project.org"
  },
  "runtimePlatform": "R version 4.4.3 (2025-02-28)",
  "provider": {
    "@id": "https://cran.r-project.org",
    "@type": "Organization",
    "name": "Comprehensive R Archive Network (CRAN)",
    "url": "https://cran.r-project.org"
  },
  "author": [
    {
      "@type": "Person",
      "givenName": "Mark",
      "familyName": "Padgham",
      "email": "mark.padgham@email.com"
    },
    {
      "@type": "Person",
      "givenName": "Michael D.",
      "familyName": "Sumner"
    }
  ],
  "copyrightHolder": [
    {
      "@type": "Person",
      "givenName": "Charles F.F",
      "familyName": "Karney"
    }
  ],
  "maintainer": [
    {
      "@type": "Person",
      "givenName": "Mark",
      "familyName": "Padgham",
      "email": "mark.padgham@email.com"
    }
  ],
  "softwareSuggestions": [
    {
      "@type": "SoftwareApplication",
      "identifier": "knitr",
      "name": "knitr",
      "provider": {
        "@id": "https://cran.r-project.org",
        "@type": "Organization",
        "name": "Comprehensive R Archive Network (CRAN)",
        "url": "https://cran.r-project.org"
      },
      "sameAs": "https://CRAN.R-project.org/package=knitr"
    },
    {
      "@type": "SoftwareApplication",
      "identifier": "rmarkdown",
      "name": "rmarkdown",
      "provider": {
        "@id": "https://cran.r-project.org",
        "@type": "Organization",
        "name": "Comprehensive R Archive Network (CRAN)",
        "url": "https://cran.r-project.org"
      },
      "sameAs": "https://CRAN.R-project.org/package=rmarkdown"
    },
    {
      "@type": "SoftwareApplication",
      "identifier": "testthat",
      "name": "testthat",
      "provider": {
        "@id": "https://cran.r-project.org",
        "@type": "Organization",
        "name": "Comprehensive R Archive Network (CRAN)",
        "url": "https://cran.r-project.org"
      },
      "sameAs": "https://CRAN.R-project.org/package=testthat"
    }
  ],
  "softwareRequirements": {
    "SystemRequirements": {}
  },
  "fileSize": "2217.814KB",
  "releaseNotes": "https://github.com/hypertidy/geodist/blob/master/NEWS.md",
  "readme": "https://github.com/hypertidy/geodist/blob/main/README.md",
  "contIntegration": [
    "https://github.com/hypertidy/geodist/actions?query=workflow%3AR-CMD-check",
    "https://github.com/hypertidy/geodist/actions?query=workflow%3Apkgcheck",
    "https://app.codecov.io/gh/hypertidy/geodist"
  ],
  "developmentStatus": "https://www.repostatus.org/"
}

GitHub Events

Committers

Last synced: 9 months ago

Issues and Pull Requests

Last synced: 6 months ago

Packages