Giving weird coordinates to Basemap segfaults

[ohnorobo]

This may only be a problem with llc projections

Code sample that recreates the bug: https://gist.github.com/ohnorobo/2718c0238c3dcc941dab

[molinav]

9 years later, I tried to reproduce your gist and the error is still there! I tested with Python 3.11, basemap 1.3.9.

Example 1:

# Works
from mpl_toolkits.basemap import Basemap
boundaries = (42.88679, 41.187053, -69.858861, -73.508142)
north, south, east, west = boundaries
b = Basemap(projection="lcc",
            llcrnrlon=west,
            llcrnrlat=south,
            urcrnrlon=east,
            urcrnrlat=north,
            lat_0=(south+north)/2,
            lon_0=(east+west)/2)
print(b)

<mpl_toolkits.basemap.Basemap object at 0x0000024000777290>

Example 2:

# Errors correctly
from mpl_toolkits.basemap import Basemap
'''
# RuntimeError: conic lat_1 = -lat_2
'''
boundaries = (90, -90, 180, -180)
north, south, east, west = boundaries
b = Basemap(projection="lcc",
            llcrnrlon=west,
            llcrnrlat=south,
            urcrnrlon=east,
            urcrnrlat=north,
            lat_0=(south+north)/2,
            lon_0=(east+west)/2)
print(b)

---------------------------------------------------------------------------
CRSError                                  Traceback (most recent call last)
Cell In[1], line 8
      6 boundaries = (90, -90, 180, -180)
      7 north, south, east, west = boundaries
----> 8 b = Basemap(projection="lcc",
      9             llcrnrlon=west,
     10             llcrnrlat=south,
     11             urcrnrlon=east,
     12             urcrnrlat=north,
     13             lat_0=(south+north)/2,
     14             lon_0=(east+west)/2)
     15 print(b)

File C:\ProgramData\pyenv-win\pyenv-win\versions\py39\lib\site-packages\mpl_toolkits\basemap\__init__.py:990, in Basemap.__init__(self, llcrnrlon, llcrnrlat, urcrnrlon, urcrnrlat, llcrnrx, llcrnry, urcrnrx, urcrnry, width, height, projection, resolution, area_thresh, rsphere, ellps, lat_ts, lat_1, lat_2, lat_0, lon_0, lon_1, lon_2, o_lon_p, o_lat_p, k_0, no_rot, suppress_ticks, satellite_height, boundinglat, fix_aspect, anchor, celestial, round, epsg, ax)
    987     raise ValueError(_unsupported_projection % projection)
    989 # initialize proj4
--> 990 proj = Proj(projparams,self.llcrnrlon,self.llcrnrlat,self.urcrnrlon,self.urcrnrlat)
    992 # make sure axis ticks are suppressed.
    993 self.noticks = suppress_ticks

File C:\ProgramData\pyenv-win\pyenv-win\versions\py39\lib\site-packages\mpl_toolkits\basemap\proj.py:206, in Proj.__init__(self, projparams, llcrnrlon, llcrnrlat, urcrnrlon, urcrnrlat, urcrnrislatlon)
    204         raise ValueError(msg)
    205 else:
--> 206     self._proj4 = pyproj.Proj(projparams)
    207     llcrnrx, llcrnry = self(llcrnrlon,llcrnrlat)
    208     if self.projection == 'aeqd': self._fulldisk=False

File C:\ProgramData\pyenv-win\pyenv-win\versions\py39\lib\site-packages\pyproj\proj.py:109, in Proj.__init__(self, projparams, preserve_units, **kwargs)
     47 def __init__(
     48     self, projparams: Optional[Any] = None, preserve_units: bool = True, **kwargs
     49 ) -> None:
     50     """
     51     A Proj class instance is initialized with proj map projection
     52     control parameter key/value pairs. The key/value pairs can
   (...)
    107     '116.366 39.867'
    108     """
--> 109     self.crs = CRS.from_user_input(projparams, **kwargs)
    110     # make sure units are meters if preserve_units is False.
    111     if not preserve_units and "foot" in self.crs.axis_info[0].unit_name:
    112         # ignore export to PROJ string deprecation warning

File C:\ProgramData\pyenv-win\pyenv-win\versions\py39\lib\site-packages\pyproj\crs\crs.py:501, in CRS.from_user_input(cls, value, **kwargs)
    499 if isinstance(value, cls):
    500     return value
--> 501 return cls(value, **kwargs)

File C:\ProgramData\pyenv-win\pyenv-win\versions\py39\lib\site-packages\pyproj\crs\crs.py:348, in CRS.__init__(self, projparams, **kwargs)
    346     self._local.crs = projparams
    347 else:
--> 348     self._local.crs = _CRS(self.srs)

File pyproj\_crs.pyx:2366, in pyproj._crs._CRS.__init__()

CRSError: Invalid projection: +proj=lcc +R=6370997.0 +units=m +lat_0=0.0 +lon_0=0.0 +lat_1=0.0 +lat_2=0.0 +type=crs: (Internal Proj Error: proj_create: Error 1027 (Invalid value for an argument): lcc: Invalid value for lat_1 and lat_2: |lat_1 + lat_2| should be > 0)

Example 3:

# Segfaults
from mpl_toolkits.basemap import Basemap
'''
$ python breakexample.py
GEOS_ERROR: TopologyException: Input geom 1 is invalid: Self-intersection at or near point -355 89.999638676664148 at -355 89.999638676664148
[1]    62852 segmentation fault  python breakexample.py
$ python3 breakexample.py
GEOS_ERROR: b'%s'
[1]    62912 segmentation fault  python3 breakexample.py
'''
boundaries = (90, -80, 180, -170)
north, south, east, west = boundaries
b = Basemap(projection="lcc",
            llcrnrlon=west,
            llcrnrlat=south,
            urcrnrlon=east,
            urcrnrlat=north,
            lat_0=(south+north)/2,
            lon_0=(east+west)/2)
print(b)

GEOS_ERROR: b'TopologyException: Input geom 1 is invalid: Self-intersection at or near point -355 89.480615067409872 at -355 89.480615067409872'

The third error is not a segfault anymore, but it is still aborting abruptly because of an error raised inside the _geoslib extension. This GEOS_ERROR could be handled definitely in a more Pythonic way.