Round-Trip Serialization Issue: Plotting capabilities depend on projection subclass.

[observingClouds]

Hi 👋

Thank you very much for maintaining this library! I really appreciate all the work that is put into this package. Projections are tough!

Setting the scene

I like to save a projection (pyproj, cartopy) following the CF-conventions, read it back in and visualize.
WKT2 seems to be one of the best formats to store the complete CRS without loosing information such that one needs to add the optional crs_wkt attribute to the single-property attributes in a CF compliant file. Using cartopy which supports WKT strings seems to be the natural choice to visualize such projections.

The problem

Doing a roundtrip from a Cartopy projection to a WKT string and back is expected to work as the WKT format is lossless (opposed to e.g. proj4). The following code snippet confirms this:

>>> import cartopy.crs as ccrs
>>> from pyproj import CRS

>>> forward = lambda c: CRS.from_user_input(c).to_wkt()
>>> backward = lambda c: ccrs.CRS(CRS.from_wkt(c))
>>> backward(forward(ccrs.PlateCarree())) == ccrs.PlateCarree()
True

This works as expected. However, during this conversion, some properties of the PlateCarree() object are lost. Some of these properties are necessary for e.g. plotting:

>>> import matplotlib.pyplot as plt

>>> fig, ax = plt.subplots(subplot_kw=dict(projection=ccrs.PlateCarree()))  # works
>>>  fig, ax = plt.subplots(subplot_kw=dict(projection=backward(forward(PlateCarree()))))  # fails as `_as_mpl_axes` is missing
# See traceback 1

Converting the CRS to a Projection does not succeed either:

>>> fig, ax = plt.subplots(subplot_kw=dict(projection=ccrs.Projection(backward(forward(PlateCarree())))))  # fails with self.bounds not implemented
# See traceback 2

Partly solution

Digging a bit deeper, providing the BBOX argument to a WKT string can fix part of the issue:

>>> ccrs.PlateCarree().to_wkt()
'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["Unknown based on WGS 84 ellipsoid",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",METHOD["Equidistant Cylindrical",ID["EPSG",1028]],PARAMETER["Latitude of 1st standard parallel",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["unknown",111319.490793274],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["unknown",111319.490793274],ID["EPSG",8807]]],CS[Cartesian,3],AXIS["(E)",east,ORDER[1],LENGTHUNIT["unknown",111319.490793274]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["unknown",111319.490793274]],AXIS["ellipsoidal height (h)",up,ORDER[3],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]'
>>> wkt = PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["Unknown based on WGS 84 ellipsoid",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",METHOD["Equidistant Cylindrical",ID["EPSG",1028]],PARAMETER["Latitude of 1st standard parallel",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["unknown",111319.490793274],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["unknown",111319.490793274],ID["EPSG",8807]]],CS[Cartesian,3],AXIS["(E)",east,ORDER[1],LENGTHUNIT["unknown",111319.490793274]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["unknown",111319.490793274]],AXIS["ellipsoidal height (h)",up,ORDER[3],LENGTHUNIT["metre",1,ID["EPSG",9001]]], USAGE[BBOX[-90,-180,90,180]]]
>>> proj = pyproj.crs.CRS.from_wkt(wkt)
<Projected CRS: PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["Unk ...>
Name: unknown
Axis Info [cartesian]:
- E[east]: Easting (unknown)
- N[north]: Northing (unknown)
- h[up]: Ellipsoidal height (metre)
Area of Use:
- name: undefined
- bounds: (-180.0, -90.0, 180.0, 90.0)
Coordinate Operation:
- name: unknown
- method: Equidistant Cylindrical
Datum: Unknown based on WGS 84 ellipsoid
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich
>>> fig, ax = plt.subplots(subplot_kw=dict(projection=Projection(proj)))

The problem here is that the following transformation returns NaNs for the points and therefore the bounds are set to NaN
https://github.com/SciTools/cartopy/blob/main/lib/cartopy/crs.py#L690-L692

For projections, that are defined by authorities (and include a BBOX) the workflow works directly:

import pyproj
import cartopy
import matplotlib.pyplot as plt
MercatorWKT = pyproj.crs.CRS.from_epsg("3395").to_wkt()
MercatorPyProj =  pyproj.crs.CRS.from_wkt(MercatorWKT)
MercatorCartopy = cartopy.crs.Projection(MercatorPyProj)
fig, ax = plt.subplots(subplot_kw=dict(projection=MercatorCartopy))
ax.coastlines()
fig.show()

UPDATE: the BBOX needs to be within the valid boundary of the projection (e.g. is not allowed to extend to the poles for some projection methods
This seems to be the case whenever, the projection is registered with an authority. So something (apart from e.g. the BBOX parameter) is written to WKT that actually enables cartopy to create a fully functional Projection object. What are these requirements?

Expectations

• Round-trip via lossless protocol (e.g. WKT2) creates object with same properties
• Plotting works for all CRS that are read via WKT2, especially those that are exported from a cartopy object
  • plotting works independently of bounding box defintions

Related issues

I think there are a lot of related issues that all touch on the cartopy backend being mostly proj4 opposed to WKT. I understand that this is a major refactoring undertaking, so I would like to particularly know how we could make the plotting more general, that it works with any/most Projection objects.

#153
#1911
#2099

Additional information

pyproj.__version__ == 3.7.0
cartopy.__version__ == main (0d4e39f)  # but also tested v0.24.0

Traceback 1

1 fig, ax = plt.subplots(subplot_kw=dict(projection=backward(forward(PlateCarree()))))

File ~/mambaforge/envs/pyproj_env/lib/python3.13/site-packages/matplotlib/pyplot.py:1760, in subplots(nrows, ncols, sharex, sharey, squeeze, width_ratios, height_ratios, subplot_kw, gridspec_kw, **fig_kw)
   1615 """
   1616 Create a figure and a set of subplots.
   1617 
   (...)
   1757 
   1758 """
   1759 fig = figure(**fig_kw)
-> 1760 axs = fig.subplots(nrows=nrows, ncols=ncols, sharex=sharex, sharey=sharey,
   1761                    squeeze=squeeze, subplot_kw=subplot_kw,
   1762                    gridspec_kw=gridspec_kw, height_ratios=height_ratios,
   1763                    width_ratios=width_ratios)
   1764 return fig, axs

File ~/mambaforge/envs/pyproj_env/lib/python3.13/site-packages/matplotlib/figure.py:861, in FigureBase.subplots(self, nrows, ncols, sharex, sharey, squeeze, width_ratios, height_ratios, subplot_kw, gridspec_kw)
    858     gridspec_kw['width_ratios'] = width_ratios
    860 gs = self.add_gridspec(nrows, ncols, figure=self, **gridspec_kw)
--> 861 axs = gs.subplots(sharex=sharex, sharey=sharey, squeeze=squeeze,
    862                   subplot_kw=subplot_kw)
    863 return axs

File ~/mambaforge/envs/pyproj_env/lib/python3.13/site-packages/matplotlib/gridspec.py:283, in GridSpecBase.subplots(self, sharex, sharey, squeeze, subplot_kw)
    281         subplot_kw["sharex"] = shared_with[sharex]
    282         subplot_kw["sharey"] = shared_with[sharey]
--> 283         axarr[row, col] = figure.add_subplot(
    284             self[row, col], **subplot_kw)
    286 # turn off redundant tick labeling
    287 if sharex in ["col", "all"]:

File ~/mambaforge/envs/pyproj_env/lib/python3.13/site-packages/matplotlib/figure.py:709, in FigureBase.add_subplot(self, *args, **kwargs)
    706 if (len(args) == 1 and isinstance(args[0], Integral)
    707         and 100 <= args[0] <= 999):
    708     args = tuple(map(int, str(args[0])))
--> 709 projection_class, pkw = self._process_projection_requirements(**kwargs)
    710 ax = projection_class(self, *args, **pkw)
    711 key = (projection_class, pkw)

File ~/mambaforge/envs/pyproj_env/lib/python3.13/site-packages/matplotlib/figure.py:1722, in FigureBase._process_projection_requirements(self, axes_class, polar, projection, **kwargs)
   1720         kwargs.update(**extra_kwargs)
   1721     else:
-> 1722         raise TypeError(
   1723             f"projection must be a string, None or implement a "
   1724             f"_as_mpl_axes method, not {projection!r}")
   1725 return projection_class, kwargs

TypeError: projection must be a string, None or implement a _as_mpl_axes method, not <Projected CRS: PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["Unk ...>
Name: unknown
Axis Info [cartesian]:
- E[east]: Easting (unknown)
- N[north]: Northing (unknown)
- h[up]: Ellipsoidal height (metre)
Area of Use:
- undefined
Coordinate Operation:
- name: unknown
- method: Equidistant Cylindrical
Datum: Unknown based on WGS 84 ellipsoid
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

Traceback 2

----> 1 fig, ax = plt.subplots(subplot_kw=dict(projection=ccrs.Projection(backward(forward(PlateCarree())))))

File ~/mambaforge/envs/pyproj_env/lib/python3.13/site-packages/matplotlib/pyplot.py:1760, in subplots(nrows, ncols, sharex, sharey, squeeze, width_ratios, height_ratios, subplot_kw, gridspec_kw, **fig_kw)
   1615 """
   1616 Create a figure and a set of subplots.
   1617 
   (...)
   1757 
   1758 """
   1759 fig = figure(**fig_kw)
-> 1760 axs = fig.subplots(nrows=nrows, ncols=ncols, sharex=sharex, sharey=sharey,
   1761                    squeeze=squeeze, subplot_kw=subplot_kw,
   1762                    gridspec_kw=gridspec_kw, height_ratios=height_ratios,
   1763                    width_ratios=width_ratios)
   1764 return fig, axs

File ~/mambaforge/envs/pyproj_env/lib/python3.13/site-packages/matplotlib/figure.py:861, in FigureBase.subplots(self, nrows, ncols, sharex, sharey, squeeze, width_ratios, height_ratios, subplot_kw, gridspec_kw)
    858     gridspec_kw['width_ratios'] = width_ratios
    860 gs = self.add_gridspec(nrows, ncols, figure=self, **gridspec_kw)
--> 861 axs = gs.subplots(sharex=sharex, sharey=sharey, squeeze=squeeze,
    862                   subplot_kw=subplot_kw)
    863 return axs

File ~/mambaforge/envs/pyproj_env/lib/python3.13/site-packages/matplotlib/gridspec.py:283, in GridSpecBase.subplots(self, sharex, sharey, squeeze, subplot_kw)
    281         subplot_kw["sharex"] = shared_with[sharex]
    282         subplot_kw["sharey"] = shared_with[sharey]
--> 283         axarr[row, col] = figure.add_subplot(
    284             self[row, col], **subplot_kw)
    286 # turn off redundant tick labeling
    287 if sharex in ["col", "all"]:

File ~/mambaforge/envs/pyproj_env/lib/python3.13/site-packages/matplotlib/figure.py:710, in FigureBase.add_subplot(self, *args, **kwargs)
    708         args = tuple(map(int, str(args[0])))
    709     projection_class, pkw = self._process_projection_requirements(**kwargs)
--> 710     ax = projection_class(self, *args, **pkw)
    711     key = (projection_class, pkw)
    712 return self._add_axes_internal(ax, key)

File ~/Documents/GitHub/cartopy/lib/cartopy/mpl/geoaxes.py:392, in GeoAxes.__init__(self, *args, **kwargs)
    388     raise ValueError("A GeoAxes can only be created with a "
    389                      "projection of type cartopy.crs.Projection")
    390 self.projection = projection
--> 392 super().__init__(*args, **kwargs)
    393 self.img_factories = []
    394 self._done_img_factory = False

File ~/mambaforge/envs/pyproj_env/lib/python3.13/site-packages/matplotlib/axes/_base.py:681, in _AxesBase.__init__(self, fig, facecolor, frameon, sharex, sharey, label, xscale, yscale, box_aspect, forward_navigation_events, *args, **kwargs)
    678 self.set_axisbelow(mpl.rcParams['axes.axisbelow'])
    680 self._rasterization_zorder = None
--> 681 self.clear()
    683 # funcs used to format x and y - fall back on major formatters
    684 self.fmt_xdata = None

File ~/Documents/GitHub/cartopy/lib/cartopy/mpl/geoaxes.py:572, in GeoAxes.clear(self)
    570 """Clear the current Axes and add boundary lines."""
    571 result = super().clear()
--> 572 self.__clear()
    573 return result

File ~/Documents/GitHub/cartopy/lib/cartopy/mpl/geoaxes.py:560, in GeoAxes.__clear(self)
    557 self._tight = True
    558 self.set_aspect('equal')
--> 560 self._boundary()
    562 # XXX consider a margin - but only when the map is not global...
    563 # self._xmargin = 0.15
    564 # self._ymargin = 0.15
    566 self.dataLim.intervalx = self.projection.x_limits

File ~/Documents/GitHub/cartopy/lib/cartopy/mpl/geoaxes.py:1523, in GeoAxes._boundary(self)
   1516 def _boundary(self):
   1517     """
   1518     Add the map's boundary to this GeoAxes.
   1519 
   1520     The :data:`.patch` and :data:`.spines['geo']` are updated to match.
   1521 
   1522     """
-> 1523     path = cpath.shapely_to_path(self.projection.boundary)
   1525     # Get the outline path in terms of self.transData
   1526     proj_to_data = self.projection._as_mpl_transform(self) - self.transData

File ~/Documents/GitHub/cartopy/lib/cartopy/crs.py:707, in Projection.boundary(self)
    704 @property
    705 def boundary(self):
    706     if self.bounds is None:
--> 707         raise NotImplementedError
    708     x0, x1, y0, y1 = self.bounds
    709     return sgeom.LineString([(x0, y0), (x0, y1), (x1, y1), (x1, y0),
    710                              (x0, y0)])

NotImplementedError:

Traceback 3

/Users/haukeschulz/mambaforge/envs/pyproj_env/lib/python3.13/site-packages/shapely/creation.py:119: RuntimeWarning: invalid value encountered in linestrings
  return lib.linestrings(coords, out=out, **kwargs)
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
File ~/Documents/GitHub/cartopy/lib/cartopy/crs.py:755, in Projection.domain(self)
    754 try:
--> 755     domain = self._domain
    756 except AttributeError:

AttributeError: 'Projection' object has no attribute '_domain'

During handling of the above exception, another exception occurred:

TopologicalError                          Traceback (most recent call last)
Cell In[89], line 1
----> 1 fig, ax = plt.subplots(subplot_kw=dict(projection=Projection(proj)))

File ~/mambaforge/envs/pyproj_env/lib/python3.13/site-packages/matplotlib/pyplot.py:1760, in subplots(nrows, ncols, sharex, sharey, squeeze, width_ratios, height_ratios, subplot_kw, gridspec_kw, **fig_kw)
   1615 """
   1616 Create a figure and a set of subplots.
   1617 
   (...)
   1757 
   1758 """
   1759 fig = figure(**fig_kw)
-> 1760 axs = fig.subplots(nrows=nrows, ncols=ncols, sharex=sharex, sharey=sharey,
   1761                    squeeze=squeeze, subplot_kw=subplot_kw,
   1762                    gridspec_kw=gridspec_kw, height_ratios=height_ratios,
   1763                    width_ratios=width_ratios)
   1764 return fig, axs

File ~/mambaforge/envs/pyproj_env/lib/python3.13/site-packages/matplotlib/figure.py:861, in FigureBase.subplots(self, nrows, ncols, sharex, sharey, squeeze, width_ratios, height_ratios, subplot_kw, gridspec_kw)
    858     gridspec_kw['width_ratios'] = width_ratios
    860 gs = self.add_gridspec(nrows, ncols, figure=self, **gridspec_kw)
--> 861 axs = gs.subplots(sharex=sharex, sharey=sharey, squeeze=squeeze,
    862                   subplot_kw=subplot_kw)
    863 return axs

File ~/mambaforge/envs/pyproj_env/lib/python3.13/site-packages/matplotlib/gridspec.py:283, in GridSpecBase.subplots(self, sharex, sharey, squeeze, subplot_kw)
    281         subplot_kw["sharex"] = shared_with[sharex]
    282         subplot_kw["sharey"] = shared_with[sharey]
--> 283         axarr[row, col] = figure.add_subplot(
    284             self[row, col], **subplot_kw)
    286 # turn off redundant tick labeling
    287 if sharex in ["col", "all"]:

File ~/mambaforge/envs/pyproj_env/lib/python3.13/site-packages/matplotlib/figure.py:710, in FigureBase.add_subplot(self, *args, **kwargs)
    708         args = tuple(map(int, str(args[0])))
    709     projection_class, pkw = self._process_projection_requirements(**kwargs)
--> 710     ax = projection_class(self, *args, **pkw)
    711     key = (projection_class, pkw)
    712 return self._add_axes_internal(ax, key)

File ~/Documents/GitHub/cartopy/lib/cartopy/mpl/geoaxes.py:392, in GeoAxes.__init__(self, *args, **kwargs)
    388     raise ValueError("A GeoAxes can only be created with a "
    389                      "projection of type cartopy.crs.Projection")
    390 self.projection = projection
--> 392 super().__init__(*args, **kwargs)
    393 self.img_factories = []
    394 self._done_img_factory = False

File ~/mambaforge/envs/pyproj_env/lib/python3.13/site-packages/matplotlib/axes/_base.py:681, in _AxesBase.__init__(self, fig, facecolor, frameon, sharex, sharey, label, xscale, yscale, box_aspect, forward_navigation_events, *args, **kwargs)
    678 self.set_axisbelow(mpl.rcParams['axes.axisbelow'])
    680 self._rasterization_zorder = None
--> 681 self.clear()
    683 # funcs used to format x and y - fall back on major formatters
    684 self.fmt_xdata = None

File ~/Documents/GitHub/cartopy/lib/cartopy/mpl/geoaxes.py:572, in GeoAxes.clear(self)
    570 """Clear the current Axes and add boundary lines."""
    571 result = super().clear()
--> 572 self.__clear()
    573 return result

File ~/Documents/GitHub/cartopy/lib/cartopy/mpl/geoaxes.py:560, in GeoAxes.__clear(self)
    557 self._tight = True
    558 self.set_aspect('equal')
--> 560 self._boundary()
    562 # XXX consider a margin - but only when the map is not global...
    563 # self._xmargin = 0.15
    564 # self._ymargin = 0.15
    566 self.dataLim.intervalx = self.projection.x_limits

File ~/Documents/GitHub/cartopy/lib/cartopy/mpl/geoaxes.py:1527, in GeoAxes._boundary(self)
   1525 # Get the outline path in terms of self.transData
   1526 proj_to_data = self.projection._as_mpl_transform(self) - self.transData
-> 1527 trans_path = proj_to_data.transform_path(path)
   1529 # Set the boundary - we can make use of the rectangular clipping.
   1530 self.set_boundary(trans_path)

File ~/mambaforge/envs/pyproj_env/lib/python3.13/site-packages/matplotlib/transforms.py:1610, in Transform.transform_path(self, path)
   1603 def transform_path(self, path):
   1604     """
   1605     Apply the transform to `.Path` *path*, returning a new `.Path`.
   1606 
   1607     In some cases, this transform may insert curves into the path
   1608     that began as line segments.
   1609     """
-> 1610     return self.transform_path_affine(self.transform_path_non_affine(path))

File ~/Documents/GitHub/cartopy/lib/cartopy/mpl/geoaxes.py:174, in InterProjectionTransform.transform_path_non_affine(self, src_path)
    171     return mpath.Path(self.transform(src_path.vertices))
    173 geom = cpath.path_to_shapely(src_path)
--> 174 transformed_geom = self.target_projection.project_geometry(
    175     geom, self.source_projection)
    177 result = cpath.shapely_to_path(transformed_geom)
    179 # store the result in the cache for future performance boosts

File ~/Documents/GitHub/cartopy/lib/cartopy/crs.py:833, in Projection.project_geometry(self, geometry, src_crs)
    831 if not method_name:
    832     raise ValueError(f'Unsupported geometry type {geom_type!r}')
--> 833 return getattr(self, method_name)(geometry, src_crs)

File ~/Documents/GitHub/cartopy/lib/cartopy/crs.py:839, in Projection._project_line_string(self, geometry, src_crs)
    838 def _project_line_string(self, geometry, src_crs):
--> 839     return cartopy.trace.project_linear(geometry, src_crs, self)

File lib/cartopy/trace.pyx:568, in cartopy.trace.project_linear()

File ~/Documents/GitHub/cartopy/lib/cartopy/crs.py:757, in Projection.domain(self)
    755     domain = self._domain
    756 except AttributeError:
--> 757     domain = self._domain = sgeom.Polygon(self.boundary)
    758 return domain

File ~/mambaforge/envs/pyproj_env/lib/python3.13/site-packages/shapely/geometry/polygon.py:230, in Polygon.__new__(self, shell, holes)
    228     return shell
    229 else:
--> 230     shell = LinearRing(shell)
    232 if holes is not None:
    233     if len(holes) == 0:
    234         # shapely constructor cannot handle holes=[]

File ~/mambaforge/envs/pyproj_env/lib/python3.13/site-packages/shapely/geometry/polygon.py:72, in LinearRing.__new__(self, coordinates)
     70     return coordinates
     71 elif not coordinates.is_valid:
---> 72     raise TopologicalError("An input LineString must be valid.")
     73 else:
     74     # LineString
     75     # TODO convert LineString to LinearRing more directly?
     76     coordinates = coordinates.coords

TopologicalError: An input LineString must be valid.

[greglucas]

I'm not sure I totally follow the request here. Cartopy is providing conveniences for projections and adding nice things to plot with like boundaries, but the ultimate source of projections is pyproj under the hood. So if you want to create projections and move them around, pyproj is really where the work should be happening. For Cartopy, I think the serialization would have to happen in pickling/unpickling classes, not a standard like WKT.

I might be missing something though, and we are definitely interested in making sharing projections across the various libraries as seamless as possible, so any PRs / updates on that are very welcome.

[observingClouds]

Thanks for your feedback @greglucas! The issue developed a bit after I experimented a bit more so it probably lost it's clarity there. I think one of the main issues here is that plotting arbitrary Projection objects is currently not possible because self.bounds of the Projection is needed. If this requirement could be relaxed somehow then this would open the plotting to projections that have no bounding box defined.