SciTools
contourf fails with lat/lon coordinates and transform for Stereographic projection but works correctly with projection coordinates
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
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.
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.
Thanks @guziy I confirm I have reproduced the error with our main development branch.
The specific error here is caused by this line
https://github.com/SciTools/cartopy/blob/97c6a6488d2edfd4780ac2d8df7b0014560e4917/lib/cartopy/crs.py#L1210
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.