ortho projection with lat_0>45 fails

[akubaryk]

I don't know if this issue belongs upstream, so apologies in advance if this isn't basemap's problem.

I compiled geos 3.9.0 from source and then installed Basemap from https://github.com/matplotlib/basemap/archive/master.zip -- my pip freeze shows basemap 1.2.1 and matplotlib 3.1.1. The base install is Anaconda 2019.10 (Python 3.7.4).

Creating any Basemap with projection='ortho' and lat_0>45 produces an error in the GEOS lib before segfaulting.

Example:

$ python -c "from mpl_toolkits.basemap import Basemap; m = Basemap(projection='ortho',lon_0=0,lat_0=45.1)"
GEOS_ERROR: b'IllegalArgumentException: CGAlgorithmsDD::orientationIndex encountered NaN/Inf numbers'
Segmentation fault (core dumped)

Strangely, this is a non-issue with lat_0 < -45.

I am more than happy to compile other versions of geos (I note 3.3.3 is included in this repo), or do whatever to help debug this.

Edit: Tried to reinstall basemap with the repo's geos-3.3.3 and the error is different, but the result is the same.

$ python -c "from mpl_toolkits.basemap import Basemap; m = Basemap(projection='ortho',lon_0=0,lat_0=44.9)"
$ python -c "from mpl_toolkits.basemap import Basemap; m = Basemap(projection='ortho',lon_0=0,lat_0=45.1)"
GEOS_ERROR: b'IllegalArgumentException: RobustDeterminant encountered non-finite numbers '
Segmentation fault (core dumped)

Another edit... pyproj 1.9.6 works with both of the above commands while 2.4.2-post1 does not.

[akubaryk]

pyproj 2.4.2 (without the post1) didn't work, but 2.4.1 and below seemed to from some limited testing.

Seems as though pyproj4/pyproj#494 is very relevant to this issue, and a pyproj dev claims it's intentional behavior in that issue.

After a lot of painful debugging, I found the error to be here in __init__.py:

# this is a workaround to avoid
# "GEOS_ERROR: TopologyException:
# found non-noded intersection between ..."
if not poly.is_valid():
    poly=poly.fix()

so up where b is defined, adding the last line here fixes the issue:

if name == 'gshhs' and as_polygons and self.projection in tostere:
    b[:,0], b[:,1] = maptran(b[:,0], b[:,1])
else:
    b[:,0], b[:,1] = self(b[:,0], b[:,1])
b = np.where(np.isposinf(np.where(np.isneginf(b),-1.e20,b)),1.e20,b)

And huzzah, I can do NH polar orthographic again.

I don't know if this is the fix y'all want in terms of implementation, but I'll open a pull if it is.

Just now seeing that you're deprecating basemap lol. Will simply close this.