observingClouds commented 4 days ago

Hi :wave:

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.

https://github.com/SciTools/cartopy/issues/153 https://github.com/SciTools/cartopy/issues/1911 https://github.com/SciTools/cartopy/issues/2099

Additional information

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

Traceback 1

```bash 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 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

```bash ----> 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

```bash /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 commented 1 day ago

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 commented 15 hours ago

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.

SciTools / cartopy

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