From 81d28d77d6b0d1b40d3abb9574d59b961d80cc1f Mon Sep 17 00:00:00 2001 From: Pascal Bourgault Date: Wed, 2 Jul 2025 17:51:36 -0400 Subject: [PATCH 1/6] Add basic crs accessors --- cf_xarray/accessor.py | 63 ++++++++++++++++++++++++++++++++++--------- 1 file changed, 50 insertions(+), 13 deletions(-) diff --git a/cf_xarray/accessor.py b/cf_xarray/accessor.py index 52d525b5..7d93a209 100644 --- a/cf_xarray/accessor.py +++ b/cf_xarray/accessor.py @@ -48,6 +48,13 @@ ) +try: + import pyproj + import cartopy.crs +except ImportError: + pyproj = None + + from . import parametric, sgrid from .criteria import ( _DSG_ROLES, @@ -2745,6 +2752,20 @@ def grid_mapping_names(self) -> dict[str, list[str]]: results[v].append(k) return results + @property + def crs(self): + """Cartopy CRS of the dataset's grid mapping.""" + if pyproj is None: + raise ImportError('`crs` accessor requires optional packages `pyproj` and `cartopy`.') + gmaps = list(itertools.chain(*self.grid_mapping_names.values())) + if len(gmaps) > 1: + raise ValueError("Multiple grid mappings found.") + if len(gmaps) == 0: + if 'longitude' in self: + return cartopy.crs.PlateCarree() + raise ValueError('No grid mapping nor longitude found in dataset.') + return cartopy.crs.Projection(pyproj.CRS.from_cf(self._obj[gmaps[0]].attrs)) + def decode_vertical_coords( self, *, outnames: dict[str, str] | None = None, prefix: str | None = None ) -> None: @@ -2899,6 +2920,21 @@ def formula_terms(self) -> dict[str, str]: # numpydoc ignore=SS06 terms[key] = value return terms + def _get_grid_mapping(self, ignore_missing=False) -> DataArray: + da = self._obj + + attrs_or_encoding = ChainMap(da.attrs, da.encoding) + grid_mapping = attrs_or_encoding.get("grid_mapping", None) + if not grid_mapping: + if ignore_missing: + return None + raise ValueError("No 'grid_mapping' attribute present.") + + if grid_mapping not in da._coords: + raise ValueError(f"Grid Mapping variable {grid_mapping} not present.") + + return da[grid_mapping] + @property def grid_mapping_name(self) -> str: """ @@ -2919,21 +2955,22 @@ def grid_mapping_name(self) -> str: >>> rotds.cf["temp"].cf.grid_mapping_name 'rotated_latitude_longitude' """ - - da = self._obj - - attrs_or_encoding = ChainMap(da.attrs, da.encoding) - grid_mapping = attrs_or_encoding.get("grid_mapping", None) - if not grid_mapping: - raise ValueError("No 'grid_mapping' attribute present.") - - if grid_mapping not in da._coords: - raise ValueError(f"Grid Mapping variable {grid_mapping} not present.") - - grid_mapping_var = da[grid_mapping] - + grid_mapping_var = self._get_grid_mapping() return grid_mapping_var.attrs["grid_mapping_name"] + @property + def crs(self): + """Cartopy CRS of the dataset's grid mapping.""" + if pyproj is None: + raise ImportError('`crs` accessor requires optional packages `pyproj` and `cartopy`.') + + grid_mapping_var = self._get_grid_mapping(ignore_missing=True) + if grid_mapping_var is None: + if 'longitude' in self: + return cartopy.crs.PlateCarree() + raise ValueError('No grid mapping nor longitude found.') + return cartopy.crs.Projection(pyproj.CRS.from_cf(grid_mapping_var.attrs)) + def __getitem__(self, key: Hashable | Iterable[Hashable]) -> DataArray: """ Index into a DataArray making use of CF attributes. From b8a797caa86404c395ff6e81c9004ea9e5ad5593 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 2 Jul 2025 21:56:59 +0000 Subject: [PATCH 2/6] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- cf_xarray/accessor.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/cf_xarray/accessor.py b/cf_xarray/accessor.py index 7d93a209..4a903f59 100644 --- a/cf_xarray/accessor.py +++ b/cf_xarray/accessor.py @@ -49,8 +49,8 @@ try: - import pyproj import cartopy.crs + import pyproj except ImportError: pyproj = None @@ -2756,14 +2756,16 @@ def grid_mapping_names(self) -> dict[str, list[str]]: def crs(self): """Cartopy CRS of the dataset's grid mapping.""" if pyproj is None: - raise ImportError('`crs` accessor requires optional packages `pyproj` and `cartopy`.') + raise ImportError( + "`crs` accessor requires optional packages `pyproj` and `cartopy`." + ) gmaps = list(itertools.chain(*self.grid_mapping_names.values())) if len(gmaps) > 1: raise ValueError("Multiple grid mappings found.") if len(gmaps) == 0: - if 'longitude' in self: + if "longitude" in self: return cartopy.crs.PlateCarree() - raise ValueError('No grid mapping nor longitude found in dataset.') + raise ValueError("No grid mapping nor longitude found in dataset.") return cartopy.crs.Projection(pyproj.CRS.from_cf(self._obj[gmaps[0]].attrs)) def decode_vertical_coords( @@ -2962,13 +2964,15 @@ def grid_mapping_name(self) -> str: def crs(self): """Cartopy CRS of the dataset's grid mapping.""" if pyproj is None: - raise ImportError('`crs` accessor requires optional packages `pyproj` and `cartopy`.') + raise ImportError( + "`crs` accessor requires optional packages `pyproj` and `cartopy`." + ) grid_mapping_var = self._get_grid_mapping(ignore_missing=True) if grid_mapping_var is None: - if 'longitude' in self: + if "longitude" in self: return cartopy.crs.PlateCarree() - raise ValueError('No grid mapping nor longitude found.') + raise ValueError("No grid mapping nor longitude found.") return cartopy.crs.Projection(pyproj.CRS.from_cf(grid_mapping_var.attrs)) def __getitem__(self, key: Hashable | Iterable[Hashable]) -> DataArray: From 67a03ad8b7f6a7fba1f8c106ec4d7f824d8d644c Mon Sep 17 00:00:00 2001 From: Pascal Bourgault Date: Thu, 3 Jul 2025 09:53:27 -0400 Subject: [PATCH 3/6] Better errors - add tests - better latlon guess --- cf_xarray/accessor.py | 9 +++++---- cf_xarray/tests/__init__.py | 1 + cf_xarray/tests/test_accessor.py | 27 +++++++++++++++++++++++++++ cf_xarray/utils.py | 13 +++++++++++++ ci/environment.yml | 2 ++ 5 files changed, 48 insertions(+), 4 deletions(-) diff --git a/cf_xarray/accessor.py b/cf_xarray/accessor.py index 7d93a209..94515179 100644 --- a/cf_xarray/accessor.py +++ b/cf_xarray/accessor.py @@ -83,6 +83,7 @@ always_iterable, emit_user_level_warning, invert_mappings, + is_latitude_longitude, parse_cell_methods_attr, parse_cf_standard_name_table, ) @@ -2761,9 +2762,9 @@ def crs(self): if len(gmaps) > 1: raise ValueError("Multiple grid mappings found.") if len(gmaps) == 0: - if 'longitude' in self: + if is_latitude_longitude(self._obj): return cartopy.crs.PlateCarree() - raise ValueError('No grid mapping nor longitude found in dataset.') + raise ValueError('No grid mapping found and dataset guessed as not latitude_longitude.') return cartopy.crs.Projection(pyproj.CRS.from_cf(self._obj[gmaps[0]].attrs)) def decode_vertical_coords( @@ -2966,9 +2967,9 @@ def crs(self): grid_mapping_var = self._get_grid_mapping(ignore_missing=True) if grid_mapping_var is None: - if 'longitude' in self: + if is_latitude_longitude(self._obj): return cartopy.crs.PlateCarree() - raise ValueError('No grid mapping nor longitude found.') + raise ValueError('No grid mapping found and dataset guesses as not latitude_longitude.') return cartopy.crs.Projection(pyproj.CRS.from_cf(grid_mapping_var.attrs)) def __getitem__(self, key: Hashable | Iterable[Hashable]) -> DataArray: diff --git a/cf_xarray/tests/__init__.py b/cf_xarray/tests/__init__.py index 8c83df3a..8665233e 100644 --- a/cf_xarray/tests/__init__.py +++ b/cf_xarray/tests/__init__.py @@ -69,3 +69,4 @@ def LooseVersion(vstring): has_pooch, requires_pooch = _importorskip("pooch") _, requires_rich = _importorskip("rich") has_regex, requires_regex = _importorskip("regex") +has_cartopy, requires_cartopy = _importorskip("cartopy") diff --git a/cf_xarray/tests/test_accessor.py b/cf_xarray/tests/test_accessor.py index c3f7005a..381d2340 100644 --- a/cf_xarray/tests/test_accessor.py +++ b/cf_xarray/tests/test_accessor.py @@ -41,6 +41,7 @@ ) from . import ( raise_if_dask_computes, + requires_cartopy, requires_cftime, requires_pint, requires_pooch, @@ -1084,6 +1085,32 @@ def test_bad_grid_mapping_attribute(): ds.cf.get_associated_variable_names("temp", error=False) +@requires_cartopy +def test_crs() -> None: + from pyproj import CRS + import cartopy.crs as ccrs + + # Dataset with explicit grid mapping + # ccrs.RotatedPole is not the same as CRS.from_cf(rotated_pole)... + # They are equivalent though, but specified differently + exp = ccrs.Projection(CRS.from_cf(rotds.rotated_pole.attrs)) + assert rotds.cf.crs == exp + with pytest.raises(ValueError, match='Grid Mapping variable rotated_pole not present'): + rotds.temp.cf.crs + assert rotds.set_coords('rotated_pole').temp.cf.crs == exp + + # Dataset with regular latlon (no grid mapping ) + exp = ccrs.PlateCarree() + assert forecast.cf.crs == exp + assert forecast.sst.cf.crs == exp + + # Dataset with no grid mapping specified but not on latlon (error) + with pytest.raises(ValueError, match="No grid mapping found"): + mollwds.cf.crs + with pytest.raises(ValueError, match="No grid mapping found"): + mollwds.lon_bounds.cf.crs + + def test_docstring() -> None: assert "One of ('X'" in airds.cf.groupby.__doc__ assert "Time variable accessor e.g. 'T.month'" in airds.cf.groupby.__doc__ diff --git a/cf_xarray/utils.py b/cf_xarray/utils.py index bdc2605e..3a867539 100644 --- a/cf_xarray/utils.py +++ b/cf_xarray/utils.py @@ -193,3 +193,16 @@ def emit_user_level_warning(message, category=None): """Emit a warning at the user level by inspecting the stack trace.""" stacklevel = find_stack_level() warnings.warn(message, category=category, stacklevel=stacklevel) + + +def is_latitude_longitude(ds): + """ + A dataset is probably using the latitude_longitude grid mapping implicitly if + - it has both longitude and latitude variables + - they are 1D (so either a list of points or a regular grid) + """ + return ( + 'longitude' in ds.cf and 'latitude' in ds.cf + and ds.cf['longitude'].ndim == 1 + and ds.cf['latitude'].ndim == 1 + ) diff --git a/ci/environment.yml b/ci/environment.yml index bcfb0780..83f13c93 100644 --- a/ci/environment.yml +++ b/ci/environment.yml @@ -6,6 +6,7 @@ dependencies: - pytest - pytest-xdist - dask + - cartopy - flox - lxml - matplotlib-base @@ -13,6 +14,7 @@ dependencies: - pandas - pint - pooch + - pyproj - regex - rich - pooch From 97e8fd68ea43ed24dd5606569f2bcbfbebee308b Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 3 Jul 2025 13:54:31 +0000 Subject: [PATCH 4/6] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- cf_xarray/accessor.py | 8 ++++++-- cf_xarray/tests/test_accessor.py | 8 +++++--- cf_xarray/utils.py | 7 ++++--- 3 files changed, 15 insertions(+), 8 deletions(-) diff --git a/cf_xarray/accessor.py b/cf_xarray/accessor.py index 4311ec6b..9621da9d 100644 --- a/cf_xarray/accessor.py +++ b/cf_xarray/accessor.py @@ -2766,7 +2766,9 @@ def crs(self): if len(gmaps) == 0: if is_latitude_longitude(self._obj): return cartopy.crs.PlateCarree() - raise ValueError('No grid mapping found and dataset guessed as not latitude_longitude.') + raise ValueError( + "No grid mapping found and dataset guessed as not latitude_longitude." + ) return cartopy.crs.Projection(pyproj.CRS.from_cf(self._obj[gmaps[0]].attrs)) def decode_vertical_coords( @@ -2973,7 +2975,9 @@ def crs(self): if grid_mapping_var is None: if is_latitude_longitude(self._obj): return cartopy.crs.PlateCarree() - raise ValueError('No grid mapping found and dataset guesses as not latitude_longitude.') + raise ValueError( + "No grid mapping found and dataset guesses as not latitude_longitude." + ) return cartopy.crs.Projection(pyproj.CRS.from_cf(grid_mapping_var.attrs)) def __getitem__(self, key: Hashable | Iterable[Hashable]) -> DataArray: diff --git a/cf_xarray/tests/test_accessor.py b/cf_xarray/tests/test_accessor.py index 381d2340..3cf02cb0 100644 --- a/cf_xarray/tests/test_accessor.py +++ b/cf_xarray/tests/test_accessor.py @@ -1087,17 +1087,19 @@ def test_bad_grid_mapping_attribute(): @requires_cartopy def test_crs() -> None: - from pyproj import CRS import cartopy.crs as ccrs + from pyproj import CRS # Dataset with explicit grid mapping # ccrs.RotatedPole is not the same as CRS.from_cf(rotated_pole)... # They are equivalent though, but specified differently exp = ccrs.Projection(CRS.from_cf(rotds.rotated_pole.attrs)) assert rotds.cf.crs == exp - with pytest.raises(ValueError, match='Grid Mapping variable rotated_pole not present'): + with pytest.raises( + ValueError, match="Grid Mapping variable rotated_pole not present" + ): rotds.temp.cf.crs - assert rotds.set_coords('rotated_pole').temp.cf.crs == exp + assert rotds.set_coords("rotated_pole").temp.cf.crs == exp # Dataset with regular latlon (no grid mapping ) exp = ccrs.PlateCarree() diff --git a/cf_xarray/utils.py b/cf_xarray/utils.py index 3a867539..213c219f 100644 --- a/cf_xarray/utils.py +++ b/cf_xarray/utils.py @@ -202,7 +202,8 @@ def is_latitude_longitude(ds): - they are 1D (so either a list of points or a regular grid) """ return ( - 'longitude' in ds.cf and 'latitude' in ds.cf - and ds.cf['longitude'].ndim == 1 - and ds.cf['latitude'].ndim == 1 + "longitude" in ds.cf + and "latitude" in ds.cf + and ds.cf["longitude"].ndim == 1 + and ds.cf["latitude"].ndim == 1 ) From 019f03520d6aa5d0056b59006539f56fb06f9130 Mon Sep 17 00:00:00 2001 From: Pascal Bourgault Date: Fri, 4 Jul 2025 09:56:26 -0400 Subject: [PATCH 5/6] changes after review --- cf_xarray/accessor.py | 12 +++++++----- cf_xarray/utils.py | 11 ++++++----- pyproject.toml | 1 + 3 files changed, 14 insertions(+), 10 deletions(-) diff --git a/cf_xarray/accessor.py b/cf_xarray/accessor.py index 9621da9d..b5e027f5 100644 --- a/cf_xarray/accessor.py +++ b/cf_xarray/accessor.py @@ -1,6 +1,7 @@ from __future__ import annotations import functools +import importlib import inspect import itertools import re @@ -48,10 +49,11 @@ ) -try: +if importlib.util.find_spec("cartopy"): + # pyproj is a dep of cartopy import cartopy.crs import pyproj -except ImportError: +else: pyproj = None @@ -2754,7 +2756,7 @@ def grid_mapping_names(self) -> dict[str, list[str]]: return results @property - def crs(self): + def cartopy_crs(self): """Cartopy CRS of the dataset's grid mapping.""" if pyproj is None: raise ImportError( @@ -2925,7 +2927,7 @@ def formula_terms(self) -> dict[str, str]: # numpydoc ignore=SS06 terms[key] = value return terms - def _get_grid_mapping(self, ignore_missing=False) -> DataArray: + def _get_grid_mapping(self, ignore_missing=False) -> DataArray | None: da = self._obj attrs_or_encoding = ChainMap(da.attrs, da.encoding) @@ -2964,7 +2966,7 @@ def grid_mapping_name(self) -> str: return grid_mapping_var.attrs["grid_mapping_name"] @property - def crs(self): + def cartopy_crs(self): """Cartopy CRS of the dataset's grid mapping.""" if pyproj is None: raise ImportError( diff --git a/cf_xarray/utils.py b/cf_xarray/utils.py index 213c219f..e4292b24 100644 --- a/cf_xarray/utils.py +++ b/cf_xarray/utils.py @@ -198,12 +198,13 @@ def emit_user_level_warning(message, category=None): def is_latitude_longitude(ds): """ A dataset is probably using the latitude_longitude grid mapping implicitly if - - it has both longitude and latitude variables + - it has both longitude and latitude coordinates - they are 1D (so either a list of points or a regular grid) """ + coords = ds.cf.coordinates return ( - "longitude" in ds.cf - and "latitude" in ds.cf - and ds.cf["longitude"].ndim == 1 - and ds.cf["latitude"].ndim == 1 + "longitude" in coords + and "latitude" in coords + and ds[coords["longitude"][0]].ndim == 1 + and ds[coords["latitude"][0]].ndim == 1 ) diff --git a/pyproject.toml b/pyproject.toml index 2269b133..93a470d1 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -137,6 +137,7 @@ enable_error_code = ["ignore-without-code", "redundant-expr", "truthy-bool"] [[tool.mypy.overrides]] module=[ + "cartopy", "cftime", "pandas", "pooch", From 3d7d9c7a019161ae0ccbb60ac872cb72e5e07894 Mon Sep 17 00:00:00 2001 From: Pascal Bourgault Date: Fri, 4 Jul 2025 10:14:41 -0400 Subject: [PATCH 6/6] Plotting draft - Broken maybe from cartopy issue --- cf_xarray/accessor.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/cf_xarray/accessor.py b/cf_xarray/accessor.py index b5e027f5..8611c322 100644 --- a/cf_xarray/accessor.py +++ b/cf_xarray/accessor.py @@ -1102,6 +1102,10 @@ def _plot_wrapper(*args, **kwargs): func.__name__ == "wrapper" and (kwargs.get("hue") or self._obj.ndim == 1) ) + is_grid_plot = (func.__name__ in ["contour", "countourf", "pcolormsh"]) or ( + func.__name__ == "wrapper" + and (self._obj.ndim - sum(["col" in kwargs, "row" in kwargs])) == 2 + ) if is_line_plot: hue = kwargs.get("hue") if "x" not in kwargs and "y" not in kwargs: @@ -1111,6 +1115,20 @@ def _plot_wrapper(*args, **kwargs): else: kwargs = self._process_x_or_y(kwargs, "x", skip=kwargs.get("y")) kwargs = self._process_x_or_y(kwargs, "y", skip=kwargs.get("x")) + if is_grid_plot and pyproj is not None: + from cartopy.mpl.geoaxes import GeoAxes + + ax = kwargs.get("ax") + if ax is None or isinstance(ax, GeoAxes): + try: + kwargs["transform"] = self._obj.cf.cartopy_crs + except ValueError: + pass + else: + if ax is None: + kwargs.setdefault("subplot_kws", {}).setdefault( + "projection", kwargs["transform"] + ) # Now set some nice properties kwargs = self._set_axis_props(kwargs, "x")