contourf fails with lat/lon coordinates and transform for Stereographic projection but works correctly with projection coordinates

[guziy]

Description

I am trying to plot a field on a stereographic projection with contourf. The contourf fails if lat/lon are provided along with the transformation. It seems to be working if I project the coordinates first and pass the projected coordinates to contourf. If I convert the data array into a plain numpy array contourf does not fail but the plot is wrong.

Below I provide the input file along with the code to reproduce the problem.

Code to reproduce

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import numpy as np
import pickle

from pathlib import Path
import xarray

def get_data():
    """
    read data
    """

    nc_pth = Path("2024010400-dump.nc")

    if not nc_pth.exists():
        with open('2024010400-dump.dat','rb') as inFile:
            XY, dataPlot, levels = pickle.load(inFile)

        lons, lats = XY["native"]
        
        ds = xarray.Dataset()
        
        ds["dataPlot"] = xarray.DataArray(dataPlot,
                                          dims=("x", "y"),
                                          coords={"lon": (("x", "y"), lons),
                                                  "lat": (("x", "y"), lats)})
        
        ds["levels"] = xarray.DataArray(levels, dims=("levels", ))
        
        ds.to_netcdf(nc_pth)
    else:
        with xarray.open_dataset(nc_pth) as ds:
            XY = {"native": [ds[k].data for k in ["lon", "lat"]]}
            dataPlot = ds["dataPlot"].to_masked_array()
            levels = ds["levels"].data

    return dataPlot, XY, levels
    

dataPlot, XY, levels = get_data()

    
# plot projection
plt_prj = ccrs.Stereographic(central_longitude=301.0, central_latitude=50.5)

# source data projection
src_prj = ccrs.PlateCarree()


coord_prj = plt_prj.transform_points(src_prj, XY['native'][0], XY['native'][1])


# works if using pre-projected coordinates
fig1 = plt.figure()
ax = fig1.add_axes([0.1,0.1,0.8,0.8], projection=plt_prj)
cs = ax.contourf(coord_prj[..., 0], coord_prj[..., 1], 
            dataPlot,
            levels=levels,
            extend="both")

ax.set_extent([XY['native'][0].min(), 
               XY['native'][0].max(), 
               XY['native'][1].min(), 
               XY['native'][1].max()])
ax.coastlines()
cb = plt.colorbar(cs, extend="both")
fig1.savefig("works.png")

# wrong plot, when converting masked array to a plain numpy array
fig1 = plt.figure()
ax = fig1.add_axes([0.1,0.1,0.8,0.8], projection=plt_prj)
to_plot = np.where(~dataPlot.mask, dataPlot, 0.)
cs = ax.contourf(XY['native'][0], XY['native'][1], 
            to_plot,
            levels=levels,
            extend="both", transform=src_prj)

ax.set_extent([XY['native'][0].min(), 
               XY['native'][0].max(), 
               XY['native'][1].min(), 
               XY['native'][1].max()])

ax.coastlines()
cb = plt.colorbar(cs, extend="both")
fig1.savefig("wrong.png")

# fails
fig1 = plt.figure()
ax = fig1.add_axes([0.1,0.1,0.8,0.8], projection=plt_prj)
cs = ax.contourf(XY['native'][0], XY['native'][1], 
                dataPlot,
                levels=levels,
                extend="both", transform=src_prj)

input file 2024010400-dump.nc.zip

Traceback

Traceback (most recent call last):
  File "/fs/homeu2/eccc/cmd/cmde/olh001/Python/test/debug-cartopy/test.py", line 58, in <module>
    cs = ax.contourf(XY['native'][0], XY['native'][1],
  File "/fs/ssm/eccc/cmd/cmds/env/python/py310_2023.07.28_all/py310/lib/python3.10/site-packages/cartopy/mpl/geoaxes.py", line 318, in wrapper
    return func(self, *args, **kwargs)
  File "/fs/ssm/eccc/cmd/cmds/env/python/py310_2023.07.28_all/py310/lib/python3.10/site-packages/cartopy/mpl/geoaxes.py", line 362, in wrapper
    return func(self, *args, **kwargs)
  File "/fs/ssm/eccc/cmd/cmds/env/python/py310_2023.07.28_all/py310/lib/python3.10/site-packages/cartopy/mpl/geoaxes.py", line 1665, in contourf
    bboxes = [col.get_datalim(self.transData)
  File "/fs/ssm/eccc/cmd/cmds/env/python/py310_2023.07.28_all/py310/lib/python3.10/site-packages/cartopy/mpl/geoaxes.py", line 1665, in <listcomp>
    bboxes = [col.get_datalim(self.transData)
  File "/fs/ssm/eccc/cmd/cmds/env/python/py310_2023.07.28_all/py310/lib/python3.10/site-packages/matplotlib/collections.py", line 267, in get_datalim
    paths = [transform.transform_path_non_affine(p) for p in paths]
  File "/fs/ssm/eccc/cmd/cmds/env/python/py310_2023.07.28_all/py310/lib/python3.10/site-packages/matplotlib/collections.py", line 267, in <listcomp>
    paths = [transform.transform_path_non_affine(p) for p in paths]
  File "/fs/ssm/eccc/cmd/cmds/env/python/py310_2023.07.28_all/py310/lib/python3.10/site-packages/matplotlib/transforms.py", line 2436, in transform_path_non_affine
    return self._a.transform_path_non_affine(path)
  File "/fs/ssm/eccc/cmd/cmds/env/python/py310_2023.07.28_all/py310/lib/python3.10/site-packages/cartopy/mpl/geoaxes.py", line 190, in transform_path_non_affine
    proj_geom = self.target_projection.project_geometry(
  File "/fs/ssm/eccc/cmd/cmds/env/python/py310_2023.07.28_all/py310/lib/python3.10/site-packages/cartopy/crs.py", line 808, in project_geometry
    return getattr(self, method_name)(geometry, src_crs)
  File "/fs/ssm/eccc/cmd/cmds/env/python/py310_2023.07.28_all/py310/lib/python3.10/site-packages/cartopy/crs.py", line 963, in _project_polygon
    return self._rings_to_multi_polygon(rings, is_ccw)
  File "/fs/ssm/eccc/cmd/cmds/env/python/py310_2023.07.28_all/py310/lib/python3.10/site-packages/cartopy/crs.py", line 1224, in _rings_to_multi_polygon
    multi_poly = sgeom.MultiPolygon(polygon_bits)
  File "/fs/ssm/eccc/cmd/cmds/env/python/py310_2023.07.28_all/py310/lib/python3.10/site-packages/shapely/geometry/multipolygon.py", line 73, in __new__
    raise ValueError("Sequences of multi-polygons are not valid arguments")
ValueError: Sequences of multi-polygons are not valid arguments

Full environment definition

Operating system

Linux

Cartopy version

0.21.1

conda list

pip list

[rcomer]

Hello @guziy, thanks for the report and the reproducing code.

Unfortunately pickle files are not portable across python environments, so it is unlikely that we will be able to run the reproducing example as-is. Could you save in some other format? NetCDF should be fine if you are used to using that.

It would also be useful to know what versions of other packages you are using, particularly shapely and matplotlib. I note that the error is being raised in shapely.

[guziy]

Here are the versions of matplotlib and shapely:

Python 3.10.12 | packaged by conda-forge | (main, Jun 23 2023, 22:40:32) [GCC 12.3.0]
Type 'copyright', 'credits' or 'license' for more information
IPython 8.14.0 -- An enhanced Interactive Python. Type '?' for help.

In [1]: import matplotlib, shapely

In [2]: matplotlib.__version__
Out[2]: '3.6.3'

In [3]: shapely.__version__
Out[3]: '2.0.1'

I've updated the code and uploaded the input netcdf file.

[rcomer]

Thanks @guziy I confirm I have reproduced the error with our main development branch.

[rcomer]

The specific error here is caused by this line

cartopy/lib/cartopy/crs.py

Line 1210 in 97c6a64

polygon = sgeom.Polygon(ring).buffer(0)

It seems that adding a zero buffer to a Polygon can create a MultiPolygon.

I tried this

diff --git a/lib/cartopy/crs.py b/lib/cartopy/crs.py
index ba505dab..1c6c7f2d 100644
--- a/lib/cartopy/crs.py
+++ b/lib/cartopy/crs.py
@@ -1226,7 +1226,11 @@ class Projection(CRS, metaclass=ABCMeta):
                     polygon = boundary_poly.intersection(polygon)
 
                     if not polygon.is_empty:
-                        polygon_bits.append(polygon)
+                        if isinstance(polygon, sgeom.MultiPolygon):
+                            polygon_bits.extend(polygon.geoms)
+                            print('multipoly')
+                        else:
+                            polygon_bits.append(polygon)
 
         if polygon_bits:
             multi_poly = sgeom.MultiPolygon(polygon_bits)

which does get us past the error. However the plot that comes out looks like

when it should look more like

So my patch appears to treat one symptom rather than the cause of the problem.