Overview

This library provides vectorized direct and inverse geodesic methods.

The motivation was to have the accurate and fast vectorized geodesic routines for sailboat routing project. There are few python libraries, with geographiclib being the most accurate and reliable. Haversine method, which is widely used as an example of fast inverse method can be vectorized rather easily, however the errors are expected to be at least 0.5%. There are no vectorized AND accurate options.

This library is based on numpy and uses Vincenty's formulae. It is heavily based on the Movable Type Scripts blog and Javascript Geodesy code.

Vincenty's inverse algorithm is accurate, but sensitive to nearly antipodal points. One approach would be to return NaN for such points, with the assumption that they are not frequently observed in practical applications, however as this discussion nicely pointed out the package cannot be complete if it cannot handle these situations. I found that the issue can be solved by relaxing one of convergence criteria, but it results in errors up to 0.25% vs geographiclib for these points.

So, instead, this library uses the vectorized Vincenty's formulae with geographiclib as a fallback for unconverged points.

See notebook for execution time comparisons vs geographiclib.

Direct method for 100,000 points

94.9 ms ± 25 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
vs
9.79 s ± 1.4 s per loop (mean ± std. dev. of 7 runs, 1 loop each)

Inverse method for 100,000 points

1.5 s ± 504 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
vs
24.2 s ± 3.91 s per loop (mean ± std. dev. of 7 runs, 1 loop each)

Setup

Requirements

• Python 3.7+
• Numpy
• Geographiclib

Installation

Install it directly into an activated virtual environment:

pip install geovectorslib

Usage

After installation, the package can imported:

$ python
>>> from geovectorslib import direct, inverse
>>> direct(lats1, lon1, bearings, distances, ellipsoid='WGS84')
>>> inverse(lats1, lons1, lats2, lons2, ellipsoid='WGS84')

Latitudes in decimal degrees [-90; +90].
Longitudes in decimal degrees [-180; +180].
Bearings in decimal degrees [0; 360].
Distances in meters.

You can specify a different ellipsoid or leave blank to use WGS84. Available ellipsoid keywords are:

• ETRS89
• WGS84
• GRS80
• GRS67
• Intl-1924
• Clarke-1880
• Airy-1830

References

Movable Type Scripts

Geodesy

Geopy

Geographiclib

Stackoverflow discussion