Implement proper geodesic buffering

[astrojuanlu]

Expected behavior and actual behavior

When doing GeoVector.buffer the user would expect for proper geodesic buffering to happen, as described in this article:

https://www.esri.com/news/arcuser/0111/geodesic.html

Instead, the wrong result is displayed.

Steps to reproduce the problem

tl.GeoVector(
    Point(125, 39.36827914916014), WGS84_CRS
).reproject(WEB_MERCATOR_CRS).buffer(10000e3)

Operating system and installation details

Linux, latest telluric.

[astrojuanlu]

(Notice that there is the combination of two problems here, see #23)

[astrojuanlu]

A good starting point for geodesic buffering is this exquisite paper "Algorithms for geodesics" by Charles F. F. Karney, where modern algorithms to compute the direct and inverse problems of the geodesic distance are developed:

https://arxiv.org/abs/1109.4448

(discovered from https://www.wikiwand.com/en/Vincenty%27s_formulae)

An undocumented version of the direct problem is implemented in GeoPy:

https://github.com/geopy/geopy/blob/355c8f98f43c5e3fa987f072e57ca419bdb3d132/geopy/distance.py#L393-L415

and also in Pyproj:

https://jswhit.github.io/pyproj/pyproj.Geod-class.html#fwd

[astrojuanlu]

Hit by this today, I might try and fix it.

[astrojuanlu]

There is a broader question here, that is: what to do with an operation that produces a result that is not representable in the original CRS?

This is exactly the case of the geodesic buffering: imagine that we have a point in EPSG:3857 close to 85 deg latitude. If we do a big buffer, the result can contain the pole, which is not representable in EPSG:3857 anymore. Naturally we would "cut-off" this part for visualization, but what about for storing it? Options:

• Cut the data - problematic, since (I think) requires specific knowledge of the singularities of the projection at hand, and also breaks the bidirectionality.
• Return invalid data - this will fail when any reprojection is requested, and also probably will produce some unexpected results.
• Return something in a different CRS where it is always valid and let the user decide - slightly inconsistent, but might be the only viable solution (for this particular case, returning the buffering in an azimuthal equidistant projection centered in the specific point or the centroid of the polygon).

More ideas?

[astrojuanlu]

Cartopy already visualizes it correctly:

But the core problem of telluric is still that a reproject might turn a valid geometry into an invalid one:

In [6]: geod_buffer = p0.reproject(azeq_p0).buffer(5000_000)

In [7]: geod_buffer.is_valid
Out [7]: True

In [8]: geod_buffer.reproject(WEB_MERCATOR_CRS).is_valid
Out [8]: False

Complete notebook: https://nbviewer.jupyter.org/gist/Juanlu001/a01e396e2364274a4a50af61eb15f6ce

[astrojuanlu]

See also #154.