Skip to content

Add CRS accessor for cartopy #577

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 7 commits into
base: main
Choose a base branch
from
Draft

Add CRS accessor for cartopy #577

wants to merge 7 commits into from

Conversation

aulemahal
Copy link
Contributor

@aulemahal aulemahal commented Jul 2, 2025

This fixes #395.

  • Add a crs accessor to both the Dataset and DataArray accessors.
    • For Datasets, we can only return something if only one grid mapping is defined.
    • For DataArrays, the grid mapping must be a coordinate
    • For both, the absence of a grid mapping and the presence of a "longitude" entry in the dataset is assumed to mean that the grid mapping is "latitude_longitude". CF allows an absent grid mapping for this CRS.
  • Include this in the plotting function (?)
  • Add tests

@dcherian
Copy link
Contributor

dcherian commented Jul 3, 2025

Also, +1 to automatic use in the plotting function, that was my initial intent!

@aulemahal
Copy link
Contributor Author

Ugh... So it's not as easy as I thought. ccrs.Projection(pyproj.CRS.from_cf(CF)) is not equivalent to the native cartopy projection because it seems to be missing information about the limits of the projection domain. If I try to set this PR's cartopy_crs as the "projection" of a matplotlib figure and also use it as the transform of the data, I get this error:

---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
Cell In[23], line 1
----> 1 ds.Temperature_surface.isel(time=0).cf.plot()

File ~/Projets/cf-xarray/cf_xarray/accessor.py:1154, in __call__(self, *args, **kwargs)
   1145 def __call__(self, *args, **kwargs):
   1146     """
   1147     Allows .plot()
   1148     """
   1149     plot = _getattr(
   1150         obj=self._obj,
   1151         attr="plot",
   1152         accessor=self.accessor,
   1153         key_mappers=dict.fromkeys(self._keys, (_single(_get_all),)),
-> 1154     )
   1155     return self._plot_decorator(plot)(*args, **kwargs)

File <string>:40, in _plot_wrapper(*args, **kwargs)

File ~/Projets/cf-xarray/cf_xarray/accessor.py:763, in _getattr.<locals>.wrapper(*args, **kwargs)
    759 posargs, arguments = accessor._process_signature(
    760     func, args, kwargs, key_mappers
    761 )
    762 final_func = extra_decorator(func) if extra_decorator else func
--> 763 result = final_func(*posargs, **arguments)
    764 if wrap_classes and isinstance(result, _WRAPPED_CLASSES):
    765     result = _CFWrappedClass(result, accessor)

File ~/miniforge3/envs/cfxr-dev/lib/python3.13/site-packages/xarray/plot/accessor.py:48, in DataArrayPlotAccessor.__call__(self, **kwargs)
     46 @functools.wraps(dataarray_plot.plot, assigned=("__doc__", "__annotations__"))
     47 def __call__(self, **kwargs) -> Any:
---> 48     return dataarray_plot.plot(self._da, **kwargs)

File ~/miniforge3/envs/cfxr-dev/lib/python3.13/site-packages/xarray/plot/dataarray_plot.py:310, in plot(darray, row, col, col_wrap, ax, hue, subplot_kws, **kwargs)
    306     plotfunc = hist
    308 kwargs["ax"] = ax
--> 310 return plotfunc(darray, **kwargs)

File ~/miniforge3/envs/cfxr-dev/lib/python3.13/site-packages/xarray/plot/dataarray_plot.py:1604, in _plot2d.<locals>.newplotfunc(***failed resolving arguments***)
   1600 if "imshow" == plotfunc.__name__ and isinstance(aspect, str):
   1601     # forbid usage of mpl strings
   1602     raise ValueError("plt.imshow's `aspect` kwarg is not available in xarray")
-> 1604 ax = get_axis(figsize, size, aspect, ax, **subplot_kws)
   1606 primitive = plotfunc(
   1607     xplt,
   1608     yplt,
   (...)   1615     **kwargs,
   1616 )
   1618 # Label the plot with metadata

File ~/miniforge3/envs/cfxr-dev/lib/python3.13/site-packages/xarray/plot/utils.py:497, in get_axis(figsize, size, aspect, ax, **subplot_kws)
    494     raise ValueError("cannot use subplot_kws with existing ax")
    496 if ax is None:
--> 497     ax = _maybe_gca(**subplot_kws)
    499 return ax

File ~/miniforge3/envs/cfxr-dev/lib/python3.13/site-packages/xarray/plot/utils.py:513, in _maybe_gca(**subplot_kws)
    509 if f.axes:
    510     # can not pass kwargs to active axes
    511     return plt.gca()
--> 513 return plt.axes(**subplot_kws)

File ~/miniforge3/envs/cfxr-dev/lib/python3.13/site-packages/matplotlib/pyplot.py:1353, in axes(arg, **kwargs)
   1351 if arg is None:
   1352     if pos is None:
-> 1353         return fig.add_subplot(**kwargs)
   1354     else:
   1355         return fig.add_axes(pos, **kwargs)

File ~/miniforge3/envs/cfxr-dev/lib/python3.13/site-packages/matplotlib/figure.py:768, in FigureBase.add_subplot(self, *args, **kwargs)
    766         args = tuple(map(int, str(args[0])))
    767     projection_class, pkw = self._process_projection_requirements(**kwargs)
--> 768     ax = projection_class(self, *args, **pkw)
    769     key = (projection_class, pkw)
    770 return self._add_axes_internal(ax, key)

File ~/miniforge3/envs/cfxr-dev/lib/python3.13/site-packages/cartopy/mpl/geoaxes.py:407, in GeoAxes.__init__(self, *args, **kwargs)
    403     raise ValueError("A GeoAxes can only be created with a "
    404                      "projection of type cartopy.crs.Projection")
    405 self.projection = projection
--> 407 super().__init__(*args, **kwargs)
    408 self.img_factories = []
    409 self._done_img_factory = False

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

File ~/miniforge3/envs/cfxr-dev/lib/python3.13/site-packages/cartopy/mpl/geoaxes.py:587, in GeoAxes.clear(self)
    585 """Clear the current Axes and add boundary lines."""
    586 result = super().clear()
--> 587 self.__clear()
    588 return result

File ~/miniforge3/envs/cfxr-dev/lib/python3.13/site-packages/cartopy/mpl/geoaxes.py:575, in GeoAxes.__clear(self)
    572 self._tight = True
    573 self.set_aspect('equal')
--> 575 self._boundary()
    577 # XXX consider a margin - but only when the map is not global...
    578 # self._xmargin = 0.15
    579 # self._ymargin = 0.15
    581 self.dataLim.intervalx = self.projection.x_limits

File ~/miniforge3/envs/cfxr-dev/lib/python3.13/site-packages/cartopy/mpl/geoaxes.py:1538, in GeoAxes._boundary(self)
   1531 def _boundary(self):
   1532     """
   1533     Add the map's boundary to this GeoAxes.
   1534 
   1535     The :data:`.patch` and :data:`.spines['geo']` are updated to match.
   1536 
   1537     """
-> 1538     path, = cpatch.geos_to_path(self.projection.boundary)
   1540     # Get the outline path in terms of self.transData
   1541     proj_to_data = self.projection._as_mpl_transform(self) - self.transData

File ~/miniforge3/envs/cfxr-dev/lib/python3.13/site-packages/cartopy/crs.py:701, in Projection.boundary(self)
    698 @property
    699 def boundary(self):
    700     if self.bounds is None:
--> 701         raise NotImplementedError
    702     x0, x1, y0, y1 = self.bounds
    703     return sgeom.LineString([(x0, y0), (x0, y1), (x1, y1), (x1, y0),
    704                              (x0, y0)])

NotImplementedError: 

If instead I use "PlateCarree" as the projection for the figure, my rotated pole test data plots like this:
Capture d’écran du 2025-07-14 15-16-26

Here is the same plot using cartopy's RotatedPole projection instead of the "trick" implemented in this PR:
Capture d’écran du 2025-07-14 15-16-22

Thus, it seems my hopes are voided, this trick is not sufficient. Maybe it's only a problem for Rotated Pole, but that's my bread and butter. I have some code to re-implement from_cf directly in cartopy, I will try to push that instead. We'll see.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Set cartopy crs automatically
2 participants