0.26.0: set_extent over a full 360° longitude range collapses the x-limits on polar/azimuthal projections for some central longitudes (map disappears, LinAlgError at some dpi)

[aneveningglow1115]

Description

With cartopy 0.26.0, GeoAxes.set_extent() with a longitude range of exactly 360°
returns an essentially zero-width x range on azimuthal projections
(NorthPolarStereo, SouthPolarStereo, Orthographic, LambertAzimuthalEqualArea),
depending on the central longitude. cartopy 0.25.0 returns the expected extent
(the full disk) with the same code.

import cartopy.crs as ccrs
import matplotlib.pyplot as plt

ax = plt.figure().add_subplot(projection=ccrs.NorthPolarStereo(central_longitude=180))
ax.set_extent([-180, 180, 0, 89.5], crs=ccrs.PlateCarree(central_longitude=180))
print(ax.get_xlim(), ax.get_ylim())

version	x range	y range
0.25.0	(-1.2687e7, 1.2693e7), width 2.54e7 m	(-1.2620e7, 1.2714e7)
0.26.0	(-1.557e-9, -6.8e-12), width 1.5e-9 m	(5.58e4, 1.2714e7)

Observations:

• The result depends on central_longitude in an irregular way (same value passed
  to the projection and to the PlateCarree used for crs): 0, 140, 179, 181 and
  200 give the correct extent; 90, 179.9, 180, 180.1, -180, 270 and 360 collapse.
• Same collapse with crs=ccrs.PlateCarree() and [0, 360, 0, 89.5], for any
  lat_min between 0 and 60, and for SouthPolarStereo(central_longitude=180),
  Orthographic(180, 90) and LambertAzimuthalEqualArea(180, 90).
• A longitude span slightly smaller than 360° works: [-179.999, 179.999, 0, 89.5],
  [0, 359.999, 0, 89.5] or [-180, 179.9, 0, 89.5] all give width 2.54e7 m.

Consequences: since GeoAxes uses adjustable="box", the axes box also becomes
zero-width and the map is drawn as a vertical line, i.e. it disappears from the
figure. If gridlines are added, Gridliner._draw_gridliner inverts a singular
transform and savefig raises numpy.linalg.LinAlgError: Singular matrix for
some dpi values (150, 300 and 600 fail; 100, 200, 250 and 400 happen to pass,
but the map is still missing).

This looks like a side effect of the PlateCarree rework for PROJ 9.8 (#2653,
eqc + lon_0 replaced by latlong + pm): the two edges of the extent
rectangle at ±180° end up on the same meridian and the projected rectangle
degenerates to a line.

Code to reproduce

import io
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

lon = np.arange(0.0, 360.0, 2.5)
lat = np.arange(-90.0, 90.1, 2.5)
LON, LAT = np.meshgrid(lon, lat)
Z = 250.0 + 40.0 * np.cos(np.radians(LAT)) + 5.0 * np.cos(np.radians(LON))

fig = plt.figure(figsize=(6.4, 4.8))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.NorthPolarStereo(central_longitude=180.0))
ax.set_extent([-180.0, 180.0, 0.0, 89.5], crs=ccrs.PlateCarree(central_longitude=180.0))
print("xlim:", ax.get_xlim(), "ylim:", ax.get_ylim())   # x width ~1e-9 on 0.26.0
cf = ax.contourf(lon, lat, Z, levels=21, transform=ccrs.PlateCarree())
fig.colorbar(cf, ax=ax)
ax.coastlines()
ax.gridlines()
fig.savefig(io.BytesIO(), format="png", dpi=300)         # LinAlgError on 0.26.0

Traceback

xlim: (np.float64(-1.5569669666696745e-09), np.float64(-6.8393253762197455e-12)) ylim: (np.float64(55847.329866713815), np.float64(12713600.098850505))
Traceback (most recent call last):
  File "repro.py", line 15, in <module>
    fig.savefig(io.BytesIO(), format="png", dpi=300)
  File "matplotlib/figure.py", line 3490, in savefig
    self.canvas.print_figure(fname, **kwargs)
  File "matplotlib/backend_bases.py", line 2186, in print_figure
    result = print_method(
  File "matplotlib/backends/backend_agg.py", line 481, in print_png
    self._print_pil(filename_or_obj, "png", pil_kwargs, metadata)
  File "matplotlib/backends/backend_agg.py", line 429, in _print_pil
    FigureCanvasAgg.draw(self)
  File "matplotlib/backends/backend_agg.py", line 382, in draw
    self.figure.draw(self.renderer)
  File "matplotlib/figure.py", line 3257, in draw
    mimage._draw_list_compositing_images(
  File "matplotlib/image.py", line 134, in _draw_list_compositing_images
    a.draw(renderer)
  File "cartopy/mpl/geoaxes.py", line 520, in draw
    return super().draw(renderer=renderer, **kwargs)
  File "matplotlib/axes/_base.py", line 3226, in draw
    mimage._draw_list_compositing_images(
  File "matplotlib/image.py", line 134, in _draw_list_compositing_images
    a.draw(renderer)
  File "cartopy/mpl/gridliner.py", line 1227, in draw
    self._draw_gridliner(renderer=renderer)
  File "cartopy/mpl/gridliner.py", line 575, in _draw_gridliner
    lon_lim, lat_lim = self._axes_domain(nx=nx, ny=ny)
  File "cartopy/mpl/gridliner.py", line 1146, in _axes_domain
    ax_to_bkg_patch = self.axes.transAxes - self.axes.patch.get_transform()
  File "matplotlib/transforms.py", line 1463, in __sub__
    return self + other.inverted()
  File "matplotlib/transforms.py", line 2443, in inverted
    self._b.inverted(), self._a.inverted())
  File "matplotlib/transforms.py", line 2443, in inverted
    self._b.inverted(), self._a.inverted())
  File "matplotlib/transforms.py", line 1877, in inverted
    self._inverted = Affine2D(inv(mtx), shorthand_name=shorthand_name)
  File "numpy/linalg/_linalg.py", line 648, in inv
    ainv = _umath_linalg.inv(a, signature=signature)
  File "numpy/linalg/_linalg.py", line 145, in _raise_linalgerror_singular
    raise LinAlgError("Singular matrix")
numpy.linalg.LinAlgError: Singular matrix

With cartopy 0.25.0 the same script prints
xlim: (-12686968.2, 12693088.5) ylim: (-12619811.3, 12713600.1) and saves the
figure at every dpi.

Full environment definition

Operating system

macOS 15.7.5 (arm64)

Cartopy version

0.26.0 (conda-forge build np2py312h3775662_0); 0.25.0 works

conda list (relevant packages, both from conda-forge)

Failing environment:

cartopy                   0.26.0          np2py312h3775662_0
geos                      3.14.1               ha834cd8_0
matplotlib-base           3.10.8          py312h605b88b_0
numpy                     2.4.3           py312h84a4f5f_0
proj                      9.9.0                h51dee2d_0
pyproj                    3.8.0           py312h16d7083_2
pyshp                     3.1.6              pyhcf101f3_0
python                    3.12.14         hd05a0c4_3_cpython
shapely                   2.1.2           py312h0fc4a54_3

Working environment (only cartopy, proj and pyproj differ):

cartopy                   0.25.0          py312h5978115_1
geos                      3.14.1               h5afe852_0
matplotlib-base           3.10.8          py312h605b88b_0
numpy                     2.4.3           py312h84a4f5f_0
proj                      9.7.1                hfb14a63_3
pyproj                    3.7.2           py312hce17ad6_3
python                    3.12.13         h8561d8f_0_cpython
shapely                   2.1.2           py312h35cd81b_2

pip list

Not used (conda-forge environment).

[QuLogic]

Can you try with the older proj as well?

[aneveningglow1115]

Yes, I downgraded PROJ in the same conda-forge environment (macOS arm64) and re-ran exactly the same checks. The PROJ version makes no difference:

cartopy	pyproj	PROJ	full-longitude set_extent on NorthPolarStereo(180)	savefig with gridline labels
0.25.0	3.7.2	9.7.1	x width 2.54e7 m (OK)	OK at every dpi
0.26.0	3.8.0	9.9.0	x width 1.55e-9 m (collapsed)	LinAlgError at dpi 150 / 300 / 600
0.26.0	3.7.2	9.7.1	x width 1.55e-9 m (collapsed)	LinAlgError at dpi 150 / 300 / 600

The full sweep (12 central longitudes, [0, 360] via PlateCarree(), SouthPolarStereo / Orthographic / LambertAzimuthalEqualArea, savefig at 5 dpi values) produced byte-identical output with PROJ 9.9.0 and 9.7.1 apart from the version line. The set of central longitudes that collapse (90, 179.9, 180, 180.1, -180, 270, 360) is the same for both.

One more observation that may help locate it: projecting the extent rectangle itself, which is what set_extent does internally, returns an empty geometry in 0.26.0 for exactly those central longitudes:

import shapely.geometry as sgeom
import cartopy.crs as ccrs

for clon in [0, 140, 179, 180, 181, 270]:
    src = ccrs.PlateCarree(central_longitude=clon)
    dst = ccrs.NorthPolarStereo(central_longitude=clon)
    g = dst.project_geometry(sgeom.box(-180, 0, 180, 89.5), src)
    print(clon, g.is_empty, g.area)

# cartopy 0.25.0 (PROJ 9.7.1): every line is  False 5.025e+14  (the full disk)
# cartopy 0.26.0 (PROJ 9.7.1):
0   False 5.025e+14
140 False 5.025e+14
179 False 5.025e+14
180 True  0.0
181 False 5.025e+14
270 True  0.0

So the rectangle whose left and right edges sit at ±180° (in the source CRS) is dropped entirely during projection for some central longitudes, and set_extent then fall
s back to something degenerate for the x limits. The individual corner points transform fine (transform_point gives the expected values, with x ≈ 0 for both ±180° corne
rs as in 0.25.0), so it is the geometry/boundary handling rather than the point transform.