diff --git a/docs/history.rst b/docs/history.rst index 343ad776..21dc29a2 100644 --- a/docs/history.rst +++ b/docs/history.rst @@ -3,6 +3,9 @@ History Latest ------ +- ENH: Add `convention` option to `set_options()` for future multi-convention support (pull #899) +- ENH: Add read support for Zarr spatial and proj conventions (pull #XXX) +- REF: Extract CF convention logic to `_convention/cf.py` module (pull #899) 0.21.0 diff --git a/docs/rioxarray.rst b/docs/rioxarray.rst index b32f8658..8d07a038 100644 --- a/docs/rioxarray.rst +++ b/docs/rioxarray.rst @@ -27,6 +27,15 @@ rioxarray.show_versions .. autofunction:: rioxarray.show_versions +rioxarray.enum module +--------------------- + +.. automodule:: rioxarray.enum + :members: + :undoc-members: + :show-inheritance: + + rioxarray `rio` accessors -------------------------- diff --git a/rioxarray/_convention/__init__.py b/rioxarray/_convention/__init__.py new file mode 100644 index 00000000..24e1d58a --- /dev/null +++ b/rioxarray/_convention/__init__.py @@ -0,0 +1 @@ +"""Convention modules for rioxarray.""" diff --git a/rioxarray/_convention/cf.py b/rioxarray/_convention/cf.py new file mode 100644 index 00000000..7b39327a --- /dev/null +++ b/rioxarray/_convention/cf.py @@ -0,0 +1,298 @@ +""" +CF (Climate and Forecasts) convention support for rioxarray. + +This module provides functions for reading and writing geospatial metadata according to +the CF conventions: https://github.com/cf-convention/cf-conventions +""" +from typing import Optional, Tuple, Union + +import numpy +import pyproj +import rasterio.crs +import xarray +from affine import Affine + +from rioxarray._options import EXPORT_GRID_MAPPING, get_option +from rioxarray.crs import crs_from_user_input + + +def _find_grid_mapping( + obj: Union[xarray.Dataset, xarray.DataArray], + grid_mapping: Optional[str] = None, +) -> Optional[str]: + """ + Find the grid_mapping coordinate name. + + Parameters + ---------- + obj : xarray.Dataset or xarray.DataArray + Object to search for grid_mapping + grid_mapping : str, optional + Explicit grid_mapping name to use + + Returns + ------- + str or None + The grid_mapping name, or None if not found + """ + if grid_mapping is not None: + return grid_mapping + + # Try to find grid_mapping attribute on data variables + if hasattr(obj, "data_vars"): + for data_var in obj.data_vars.values(): + if "grid_mapping" in data_var.attrs: + return data_var.attrs["grid_mapping"] + if "grid_mapping" in data_var.encoding: + return data_var.encoding["grid_mapping"] + + if hasattr(obj, "attrs") and "grid_mapping" in obj.attrs: + return obj.attrs["grid_mapping"] + + if hasattr(obj, "encoding") and "grid_mapping" in obj.encoding: + return obj.encoding["grid_mapping"] + + return None + + +def read_crs( + obj: Union[xarray.Dataset, xarray.DataArray], grid_mapping: Optional[str] = None +) -> Optional[rasterio.crs.CRS]: + """ + Read CRS from CF conventions. + + Parameters + ---------- + obj : xarray.Dataset or xarray.DataArray + Object to read CRS from + grid_mapping : str, optional + Name of the grid_mapping coordinate variable + + Returns + ------- + rasterio.crs.CRS or None + CRS object, or None if not found + """ + grid_mapping = _find_grid_mapping(obj, grid_mapping) + + if grid_mapping is not None: + try: + grid_mapping_coord = obj.coords[grid_mapping] + + # Look in wkt attributes first for performance + for crs_attr in ("spatial_ref", "crs_wkt"): + try: + return crs_from_user_input(grid_mapping_coord.attrs[crs_attr]) + except KeyError: + pass + + # Look in grid_mapping CF attributes + try: + return pyproj.CRS.from_cf(grid_mapping_coord.attrs) + except (KeyError, pyproj.exceptions.CRSError): + pass + except KeyError: + # grid_mapping coordinate doesn't exist, fall through to attrs check + pass + + # Fallback: look in attrs for 'crs' + try: + return crs_from_user_input(obj.attrs["crs"]) + except KeyError: + return None + + +def read_transform( + obj: Union[xarray.Dataset, xarray.DataArray], grid_mapping: Optional[str] = None +) -> Optional[Affine]: + """ + Read transform from CF conventions (GeoTransform attribute). + + Parameters + ---------- + obj : xarray.Dataset or xarray.DataArray + Object to read transform from + grid_mapping : str, optional + Name of the grid_mapping coordinate variable + + Returns + ------- + affine.Affine or None + Transform object, or None if not found + """ + grid_mapping = _find_grid_mapping(obj, grid_mapping) + + if grid_mapping is not None: + try: + transform = numpy.fromstring( + obj.coords[grid_mapping].attrs["GeoTransform"], sep=" " + ) + # Calling .tolist() to assure the arguments are Python float and JSON serializable + return Affine.from_gdal(*transform.tolist()) + except KeyError: + pass + + # Fallback: look in attrs for 'transform' + try: + return Affine(*obj.attrs["transform"][:6]) + except KeyError: + pass + + return None + + +def read_spatial_dimensions( + obj: Union[xarray.Dataset, xarray.DataArray], +) -> Optional[Tuple[str, str]]: + """ + Read spatial dimensions from CF conventions. + + This function detects spatial dimensions based on: + 1. Standard dimension names ('x'/'y', 'longitude'/'latitude') + 2. CF coordinate attributes ('axis', 'standard_name') + + Parameters + ---------- + obj : xarray.Dataset or xarray.DataArray + Object to read spatial dimensions from + + Returns + ------- + tuple of (y_dim, x_dim) or None + Tuple of dimension names, or None if not found + """ + x_dim = None + y_dim = None + + # Check standard dimension names + if "x" in obj.dims and "y" in obj.dims: + return "y", "x" + elif "longitude" in obj.dims and "latitude" in obj.dims: + return "latitude", "longitude" + + # Look for coordinates with CF attributes + for coord in obj.coords: + # Make sure to only look in 1D coordinates + # that has the same dimension name as the coordinate + if obj.coords[coord].dims != (coord,): + continue + if (obj.coords[coord].attrs.get("axis", "").upper() == "X") or ( + obj.coords[coord].attrs.get("standard_name", "").lower() + in ("longitude", "projection_x_coordinate") + ): + x_dim = coord + elif (obj.coords[coord].attrs.get("axis", "").upper() == "Y") or ( + obj.coords[coord].attrs.get("standard_name", "").lower() + in ("latitude", "projection_y_coordinate") + ): + y_dim = coord + + if x_dim is not None and y_dim is not None: + return y_dim, x_dim + + return None + + +def write_crs( + obj: Union[xarray.Dataset, xarray.DataArray], + crs: rasterio.crs.CRS, + grid_mapping_name: str, + inplace: bool = True, +) -> Union[xarray.Dataset, xarray.DataArray]: + """ + Write CRS using CF conventions. + + Parameters + ---------- + obj : xarray.Dataset or xarray.DataArray + Object to write CRS to + crs : rasterio.crs.CRS + CRS to write + grid_mapping_name : str + Name of the grid_mapping coordinate + inplace : bool, default True + If True, modify object in place + + Returns + ------- + xarray.Dataset or xarray.DataArray + Object with CRS written + """ + obj_out = obj if inplace else obj.copy(deep=True) + + # Get original transform before modifying + transform = read_transform(obj) + + # Remove old grid mapping coordinate if exists + try: + del obj_out.coords[grid_mapping_name] + except KeyError: + pass + + # Add grid mapping coordinate + obj_out.coords[grid_mapping_name] = xarray.Variable((), 0) + grid_map_attrs = {} + if get_option(EXPORT_GRID_MAPPING): + try: + grid_map_attrs = pyproj.CRS.from_user_input(crs).to_cf() + except KeyError: + pass + + # spatial_ref is for compatibility with GDAL + crs_wkt = crs.to_wkt() + grid_map_attrs["spatial_ref"] = crs_wkt + grid_map_attrs["crs_wkt"] = crs_wkt + if transform is not None: + grid_map_attrs["GeoTransform"] = " ".join( + [str(item) for item in transform.to_gdal()] + ) + obj_out.coords[grid_mapping_name].attrs = grid_map_attrs + + # Remove old crs if exists + obj_out.attrs.pop("crs", None) + + return obj_out + + +def write_transform( + obj: Union[xarray.Dataset, xarray.DataArray], + transform: Affine, + grid_mapping_name: str, + inplace: bool = True, +) -> Union[xarray.Dataset, xarray.DataArray]: + """ + Write transform using CF conventions (GeoTransform attribute). + + Parameters + ---------- + obj : xarray.Dataset or xarray.DataArray + Object to write transform to + transform : affine.Affine + Transform to write + grid_mapping_name : str + Name of the grid_mapping coordinate + inplace : bool, default True + If True, modify object in place + + Returns + ------- + xarray.Dataset or xarray.DataArray + Object with transform written + """ + obj_out = obj if inplace else obj.copy(deep=True) + + # Delete the old attribute to prevent confusion + obj_out.attrs.pop("transform", None) + + try: + grid_map_attrs = obj_out.coords[grid_mapping_name].attrs.copy() + except KeyError: + obj_out.coords[grid_mapping_name] = xarray.Variable((), 0) + grid_map_attrs = obj_out.coords[grid_mapping_name].attrs.copy() + + grid_map_attrs["GeoTransform"] = " ".join( + [str(item) for item in transform.to_gdal()] + ) + obj_out.coords[grid_mapping_name].attrs = grid_map_attrs + + return obj_out diff --git a/rioxarray/_convention/zarr.py b/rioxarray/_convention/zarr.py new file mode 100644 index 00000000..35c62fcf --- /dev/null +++ b/rioxarray/_convention/zarr.py @@ -0,0 +1,326 @@ +""" +Zarr spatial and proj convention support for rioxarray. + +This module provides functions for reading geospatial metadata according to: +- Zarr spatial convention: https://github.com/zarr-conventions/spatial +- Zarr geo-proj convention: https://github.com/zarr-experimental/geo-proj +""" +import json +from typing import Optional, Union + +import rasterio.crs +import xarray +from affine import Affine + +from rioxarray.crs import crs_from_user_input + +# Convention identifiers +PROJ_CONVENTION = { + "schema_url": "https://raw.githubusercontent.com/zarr-experimental/geo-proj/refs/tags/v1/schema.json", + "spec_url": "https://github.com/zarr-experimental/geo-proj/blob/v1/README.md", + "uuid": "f17cb550-5864-4468-aeb7-f3180cfb622f", + "name": "proj:", + "description": "Coordinate reference system information for geospatial data", +} + +SPATIAL_CONVENTION = { + "schema_url": "https://raw.githubusercontent.com/zarr-conventions/spatial/refs/tags/v1/schema.json", + "spec_url": "https://github.com/zarr-conventions/spatial/blob/v1/README.md", + "uuid": "689b58e2-cf7b-45e0-9fff-9cfc0883d6b4", + "name": "spatial:", + "description": "Spatial coordinate information", +} + + +def has_convention_declared(attrs: dict, convention_name: str) -> bool: + """ + Check if a specific convention is declared in zarr_conventions. + + Parameters + ---------- + attrs : dict + Attributes dictionary to check + convention_name : str + Name of convention to check for (e.g., "proj:" or "spatial:") + + Returns + ------- + bool + True if convention is declared + """ + zarr_conventions = attrs.get("zarr_conventions", []) + if not isinstance(zarr_conventions, list): + return False + + for convention in zarr_conventions: + if isinstance(convention, dict) and convention.get("name") == convention_name: + return True + return False + + +def get_declared_conventions(attrs: dict) -> set: + """ + Get set of declared convention names from attrs. + + Parameters + ---------- + attrs : dict + Attributes dictionary to check + + Returns + ------- + set + Set of declared convention names (e.g., {"proj:", "spatial:"}) + """ + zarr_conventions = attrs.get("zarr_conventions", []) + if not isinstance(zarr_conventions, list): + return set() + + declared = set() + for convention in zarr_conventions: + if isinstance(convention, dict) and "name" in convention: + declared.add(convention["name"]) + + return declared + + +# ============================================================================ +# Parsing utilities +# ============================================================================ + + +def parse_spatial_transform( + spatial_transform: Union[list, tuple], +) -> Optional[Affine]: + """ + Convert spatial:transform array to Affine object. + + Parameters + ---------- + spatial_transform : list or tuple + Transform as [a, b, c, d, e, f] array + + Returns + ------- + affine.Affine or None + Affine transform object, or None if invalid + """ + if not isinstance(spatial_transform, (list, tuple)): + return None + if len(spatial_transform) != 6: + return None + try: + return Affine(*spatial_transform) + except (TypeError, ValueError): + return None + + +def parse_proj_code(proj_code: str) -> Optional[rasterio.crs.CRS]: + """ + Parse proj:code to CRS. + + Parameters + ---------- + proj_code : str + Authority code string (e.g., "EPSG:4326") + + Returns + ------- + rasterio.crs.CRS or None + CRS object, or None if invalid + """ + if not isinstance(proj_code, str): + return None + return crs_from_user_input(proj_code) + + +def parse_proj_wkt2(proj_wkt2: str) -> Optional[rasterio.crs.CRS]: + """ + Parse proj:wkt2 to CRS. + + Parameters + ---------- + proj_wkt2 : str + WKT2 string representation of CRS + + Returns + ------- + rasterio.crs.CRS or None + CRS object, or None if invalid + """ + if not isinstance(proj_wkt2, str): + return None + return rasterio.crs.CRS.from_wkt(proj_wkt2) + + +def parse_proj_projjson( + proj_projjson: Union[dict, str], +) -> Optional[rasterio.crs.CRS]: + """ + Parse proj:projjson to CRS. + + Parameters + ---------- + proj_projjson : dict or str + PROJJSON object or JSON string + + Returns + ------- + rasterio.crs.CRS or None + CRS object, or None if invalid + """ + if isinstance(proj_projjson, str): + proj_projjson = json.loads(proj_projjson) + + if not isinstance(proj_projjson, dict): + return None + + return crs_from_user_input(json.dumps(proj_projjson)) + + +# ============================================================================ +# Internal parsing helpers +# ============================================================================ + + +def _parse_crs_from_attrs( + attrs: dict, convention_check: bool = True +) -> Optional[rasterio.crs.CRS]: + """ + Parse CRS from proj: attributes with fallback priority. + + Parameters + ---------- + attrs : dict + Attributes dictionary to parse from + convention_check : bool, default True + Whether to check for convention declaration + + Returns + ------- + rasterio.crs.CRS or None + Parsed CRS object, or None if not found + """ + if convention_check and not has_convention_declared(attrs, "proj:"): + return None + + for proj_attr, parser in [ + ("proj:wkt2", parse_proj_wkt2), + ("proj:code", parse_proj_code), + ("proj:projjson", parse_proj_projjson), + ]: + try: + proj_value = attrs.get(proj_attr) + if proj_value is not None: + parsed_crs = parser(proj_value) + if parsed_crs is not None: + return parsed_crs + except (KeyError, Exception): + pass + return None + + +def _parse_transform_from_attrs( + attrs: dict, convention_check: bool = True +) -> Optional[Affine]: + """ + Parse transform from spatial: attributes. + + Parameters + ---------- + attrs : dict + Attributes dictionary to parse from + convention_check : bool, default True + Whether to check for convention declaration + + Returns + ------- + affine.Affine or None + Parsed transform object, or None if not found + """ + if convention_check and not has_convention_declared(attrs, "spatial:"): + return None + + try: + spatial_transform = attrs.get("spatial:transform") + if spatial_transform is not None: + return parse_spatial_transform(spatial_transform) + except (KeyError, Exception): + pass + return None + + +# ============================================================================ +# Public read functions +# ============================================================================ + + +def read_crs( + obj: Union[xarray.Dataset, xarray.DataArray], +) -> Optional[rasterio.crs.CRS]: + """ + Read CRS from Zarr proj: convention. + + Parameters + ---------- + obj : xarray.Dataset or xarray.DataArray + Object to read CRS from + + Returns + ------- + rasterio.crs.CRS or None + CRS object, or None if not found + """ + return _parse_crs_from_attrs(obj.attrs) + + +def read_transform( + obj: Union[xarray.Dataset, xarray.DataArray], +) -> Optional[Affine]: + """ + Read transform from Zarr spatial: convention. + + Parameters + ---------- + obj : xarray.Dataset or xarray.DataArray + Object to read transform from + + Returns + ------- + affine.Affine or None + Transform object, or None if not found + """ + return _parse_transform_from_attrs(obj.attrs) + + +def read_spatial_dimensions( + obj: Union[xarray.Dataset, xarray.DataArray], +) -> Optional[tuple[str, str]]: + """ + Read spatial dimensions from Zarr spatial: convention. + + Parameters + ---------- + obj : xarray.Dataset or xarray.DataArray + Object to read spatial dimensions from + + Returns + ------- + tuple of (y_dim, x_dim) or None + Tuple of dimension names, or None if not found + """ + # Only interpret spatial:* attributes if convention is declared + if not has_convention_declared(obj.attrs, "spatial:"): + return None + + try: + spatial_dims = obj.attrs.get("spatial:dimensions") + if spatial_dims is not None and len(spatial_dims) >= 2: + # spatial:dimensions format is ["y", "x"] or similar + y_dim_name, x_dim_name = spatial_dims[-2:] # Take last two + if y_dim_name in obj.dims and x_dim_name in obj.dims: + return y_dim_name, x_dim_name + except (KeyError, Exception): + pass + + return None diff --git a/rioxarray/_options.py b/rioxarray/_options.py index 1dc55ffa..61a60395 100644 --- a/rioxarray/_options.py +++ b/rioxarray/_options.py @@ -6,19 +6,30 @@ This file was adopted from: https://github.com/pydata/xarray # noqa Source file: https://github.com/pydata/xarray/blob/2ab0666c1fcc493b1e0ebc7db14500c427f8804e/xarray/core/options.py # noqa """ -from typing import Any +from typing import Any, Optional + +from rioxarray.enum import Convention EXPORT_GRID_MAPPING = "export_grid_mapping" SKIP_MISSING_SPATIAL_DIMS = "skip_missing_spatial_dims" +CONVENTION = "convention" -OPTIONS = { +OPTIONS: dict[str, Any] = { EXPORT_GRID_MAPPING: True, SKIP_MISSING_SPATIAL_DIMS: False, + CONVENTION: None, } OPTION_NAMES = set(OPTIONS) + +def _validate_convention(value: Optional[Convention]) -> bool: + """Validate the convention option.""" + return value is None or isinstance(value, Convention) + + VALIDATORS = { EXPORT_GRID_MAPPING: lambda choice: isinstance(choice, bool), + CONVENTION: _validate_convention, } @@ -46,6 +57,7 @@ class set_options: # pylint: disable=invalid-name .. versionadded:: 0.3.0 .. versionadded:: 0.7.0 skip_missing_spatial_dims + .. versionadded:: 0.18.0 convention Parameters ---------- @@ -60,6 +72,11 @@ class set_options: # pylint: disable=invalid-name If True, it will not perform spatial operations on variables within a :class:`xarray.Dataset` if the spatial dimensions are not found. + convention: Convention, default=None + The convention to use for reading and writing geospatial metadata. + If None, CF convention is used as the default with Zarr as fallback + when Zarr conventions are explicitly declared in the data. + See :class:`rioxarray.enum.Convention` for available options. Usage as a context manager:: diff --git a/rioxarray/enum.py b/rioxarray/enum.py new file mode 100644 index 00000000..60fd0669 --- /dev/null +++ b/rioxarray/enum.py @@ -0,0 +1,49 @@ +"""Enums for rioxarray.""" +from enum import Enum + + +class Convention(Enum): + """ + Supported geospatial metadata conventions. + + rioxarray supports conventions for storing geospatial metadata. + Currently supported: + + - CF: Climate and Forecasts convention using grid_mapping coordinates + - Zarr: Zarr spatial and proj conventions for cloud-native formats + + The convention can be set globally using set_options() or per-method + using the convention parameter. + + Examples + -------- + Set global convention: + + >>> import rioxarray + >>> from rioxarray.enum import Convention + >>> rioxarray.set_options(convention=Convention.CF) + + Use specific convention for a method: + + >>> from rioxarray.enum import Convention + >>> data.rio.write_crs("EPSG:4326", convention=Convention.CF) + + See Also + -------- + rioxarray.set_options : Set global options including convention + + References + ---------- + .. [1] CF Conventions: https://github.com/cf-convention/cf-conventions + .. [2] Zarr spatial convention: https://github.com/zarr-conventions/spatial + .. [3] Zarr geo-proj convention: https://github.com/zarr-experimental/geo-proj + """ + + #: Climate and Forecasts convention (default) + #: https://github.com/cf-convention/cf-conventions + CF = "CF" + + #: Zarr spatial and proj conventions + #: https://github.com/zarr-conventions/spatial + #: https://github.com/zarr-experimental/geo-proj + Zarr = "Zarr" diff --git a/rioxarray/rioxarray.py b/rioxarray/rioxarray.py index 17d61e14..f5ca0e9e 100644 --- a/rioxarray/rioxarray.py +++ b/rioxarray/rioxarray.py @@ -11,7 +11,6 @@ from typing import Any, Literal, Optional, Union import numpy -import pyproj import rasterio.warp import rasterio.windows import xarray @@ -22,7 +21,8 @@ from rasterio.crs import CRS from rasterio.rpc import RPC -from rioxarray._options import EXPORT_GRID_MAPPING, get_option +from rioxarray._convention import cf, zarr +from rioxarray._options import CONVENTION, get_option from rioxarray._spatial_utils import ( # noqa: F401, pylint: disable=unused-import DEFAULT_GRID_MAP, _affine_has_rotation, @@ -35,6 +35,7 @@ affine_to_coords, ) from rioxarray.crs import crs_from_user_input +from rioxarray.enum import Convention from rioxarray.exceptions import ( DimensionError, DimensionMissingCoordinateError, @@ -57,30 +58,25 @@ def __init__(self, xarray_obj: Union[xarray.DataArray, xarray.Dataset]): self._x_dim: Optional[Hashable] = None self._y_dim: Optional[Hashable] = None - # Determine the spatial dimensions of the `xarray.DataArray` - if "x" in self._obj.dims and "y" in self._obj.dims: - self._x_dim = "x" - self._y_dim = "y" - elif "longitude" in self._obj.dims and "latitude" in self._obj.dims: - self._x_dim = "longitude" - self._y_dim = "latitude" + + # Read spatial dimensions using the global convention setting + # Both conventions are tried regardless of setting (priority changes only) + convention = get_option(CONVENTION) + + spatial_dims = None + if convention is Convention.Zarr: + # Try Zarr first, then CF + spatial_dims = zarr.read_spatial_dimensions(self._obj) + if spatial_dims is None: + spatial_dims = cf.read_spatial_dimensions(self._obj) else: - # look for coordinates with CF attributes - for coord in self._obj.coords: - # make sure to only look in 1D coordinates - # that has the same dimension name as the coordinate - if self._obj.coords[coord].dims != (coord,): - continue - if (self._obj.coords[coord].attrs.get("axis", "").upper() == "X") or ( - self._obj.coords[coord].attrs.get("standard_name", "").lower() - in ("longitude", "projection_x_coordinate") - ): - self._x_dim = coord - elif (self._obj.coords[coord].attrs.get("axis", "").upper() == "Y") or ( - self._obj.coords[coord].attrs.get("standard_name", "").lower() - in ("latitude", "projection_y_coordinate") - ): - self._y_dim = coord + # Default: Try CF first, then Zarr + spatial_dims = cf.read_spatial_dimensions(self._obj) + if spatial_dims is None: + spatial_dims = zarr.read_spatial_dimensions(self._obj) + + if spatial_dims is not None: + self._y_dim, self._x_dim = spatial_dims # properties self._count: Optional[int] = None @@ -98,32 +94,29 @@ def crs(self) -> Optional[rasterio.crs.CRS]: if self._crs is not None: return None if self._crs is False else self._crs - # look in wkt attributes to avoid using - # pyproj CRS if possible for performance - for crs_attr in ("spatial_ref", "crs_wkt"): - try: - self._set_crs( - self._obj.coords[self.grid_mapping].attrs[crs_attr], - inplace=True, - ) - return self._crs - except KeyError: - pass + # Read using global convention setting + # Both conventions are tried regardless of setting (priority changes only) + convention = get_option(CONVENTION) - # look in grid_mapping - try: - self._set_crs( - pyproj.CRS.from_cf(self._obj.coords[self.grid_mapping].attrs), - inplace=True, - ) - except (KeyError, pyproj.exceptions.CRSError): - try: - # look in attrs for 'crs' - self._set_crs(self._obj.attrs["crs"], inplace=True) - except KeyError: - self._crs = False - return None - return self._crs + parsed_crs = None + if convention is Convention.Zarr: + # Try Zarr first, then CF + parsed_crs = zarr.read_crs(self._obj) + if parsed_crs is None: + parsed_crs = cf.read_crs(self._obj, self.grid_mapping) + else: + # Default: Try CF first, then Zarr + parsed_crs = cf.read_crs(self._obj, self.grid_mapping) + if parsed_crs is None: + parsed_crs = zarr.read_crs(self._obj) + + if parsed_crs is not None: + self._set_crs(parsed_crs, inplace=True) + return self._crs + + # No CRS found + self._crs = False + return None def _get_obj(self, inplace: bool) -> Union[xarray.Dataset, xarray.DataArray]: """ @@ -280,12 +273,13 @@ def write_crs( self, input_crs: Optional[Any] = None, grid_mapping_name: Optional[str] = None, + convention: Optional[Convention] = None, inplace: bool = False, ) -> Union[xarray.Dataset, xarray.DataArray]: """ - Write the CRS to the dataset in a CF compliant manner. + Write the CRS to the dataset using the specified convention. - .. warning:: The grid_mapping attribute is written to the encoding. + .. warning:: When using CF convention, the grid_mapping attribute is written to the encoding. Parameters ---------- @@ -293,14 +287,17 @@ def write_crs( Anything accepted by `rasterio.crs.CRS.from_user_input`. grid_mapping_name: str, optional Name of the grid_mapping coordinate to store the CRS information in. - Default is the grid_mapping name of the dataset. + Only used with CF convention. Default is the grid_mapping name of the dataset. + convention: Convention, optional + Convention to use for writing CRS. If None, uses the global default + from set_options(). Currently only CF convention is supported. inplace: bool, optional If True, it will write to the existing dataset. Default is False. Returns ------- :obj:`xarray.Dataset` | :obj:`xarray.DataArray`: - Modified dataset with CF compliant CRS information. + Modified dataset with CRS information. Examples -------- @@ -317,45 +314,30 @@ def write_crs( else: data_obj = self._get_obj(inplace=inplace) - # get original transform - transform = self._cached_transform() - # remove old grid maping coordinate if exists - grid_mapping_name = ( - self.grid_mapping if grid_mapping_name is None else grid_mapping_name - ) - try: - del data_obj.coords[grid_mapping_name] - except KeyError: - pass - if data_obj.rio.crs is None: raise MissingCRS( "CRS not found. Please set the CRS with 'rio.write_crs()'." ) - # add grid mapping coordinate - data_obj.coords[grid_mapping_name] = xarray.Variable((), 0) - grid_map_attrs = {} - if get_option(EXPORT_GRID_MAPPING): - try: - grid_map_attrs = pyproj.CRS.from_user_input(data_obj.rio.crs).to_cf() - except KeyError: - pass - # spatial_ref is for compatibility with GDAL - crs_wkt = data_obj.rio.crs.to_wkt() - grid_map_attrs["spatial_ref"] = crs_wkt - grid_map_attrs["crs_wkt"] = crs_wkt - if transform is not None: - grid_map_attrs["GeoTransform"] = " ".join( - [str(item) for item in transform.to_gdal()] - ) - data_obj.coords[grid_mapping_name].rio.set_attrs(grid_map_attrs, inplace=True) - # remove old crs if exists - data_obj.attrs.pop("crs", None) + # Determine which convention to use + if convention is None: + convention = get_option(CONVENTION) or Convention.CF - return data_obj.rio.write_grid_mapping( - grid_mapping_name=grid_mapping_name, inplace=True - ) + if convention is Convention.CF: + grid_mapping_name = ( + self.grid_mapping if grid_mapping_name is None else grid_mapping_name + ) + data_obj = cf.write_crs( + data_obj, + data_obj.rio.crs, + grid_mapping_name, + inplace=True, + ) + return data_obj.rio.write_grid_mapping( + grid_mapping_name=grid_mapping_name, inplace=True + ) + else: + raise ValueError(f"Unsupported convention: {convention}") def estimate_utm_crs(self, datum_name: str = "WGS 84") -> rasterio.crs.CRS: """Returns the estimated UTM CRS based on the bounds of the dataset. @@ -401,36 +383,37 @@ def estimate_utm_crs(self, datum_name: str = "WGS 84") -> rasterio.crs.CRS: def _cached_transform(self) -> Optional[Affine]: """ - Get the transform from: - 1. The GeoTransform metatada property in the grid mapping - 2. The transform attribute. + Get the transform using the global convention setting. + Both conventions are tried regardless of setting (priority changes only). """ - try: - # look in grid_mapping - transform = numpy.fromstring( - self._obj.coords[self.grid_mapping].attrs["GeoTransform"], sep=" " - ) - # Calling .tolist() to assure the arguments are Python float and JSON serializable - return Affine.from_gdal(*transform.tolist()) + convention = get_option(CONVENTION) - except KeyError: - try: - return Affine(*self._obj.attrs["transform"][:6]) - except KeyError: - pass - return None + if convention is Convention.Zarr: + # Try Zarr first, then CF + transform = zarr.read_transform(self._obj) + if transform is None: + transform = cf.read_transform(self._obj, self.grid_mapping) + return transform + else: + # Default: Try CF first, then Zarr + transform = cf.read_transform(self._obj, self.grid_mapping) + if transform is None: + transform = zarr.read_transform(self._obj) + return transform def write_transform( self, transform: Optional[Affine] = None, grid_mapping_name: Optional[str] = None, + convention: Optional[Convention] = None, inplace: bool = False, ) -> Union[xarray.Dataset, xarray.DataArray]: """ .. versionadded:: 0.0.30 - Write the GeoTransform to the dataset where GDAL can read it in. + Write the transform to the dataset using the specified convention. + For CF convention, this writes the GeoTransform to the dataset where GDAL can read it in. https://gdal.org/drivers/raster/netcdf.html#georeference Parameters @@ -439,34 +422,40 @@ def write_transform( The transform of the dataset. If not provided, it will be calculated. grid_mapping_name: str, optional Name of the grid_mapping coordinate to store the transform information in. - Default is the grid_mapping name of the dataset. + Only used with CF convention. Default is the grid_mapping name of the dataset. + convention: Convention, optional + Convention to use for writing transform. If None, uses the global default + from set_options(). Currently only CF convention is supported. inplace: bool, optional If True, it will write to the existing dataset. Default is False. Returns ------- :obj:`xarray.Dataset` | :obj:`xarray.DataArray`: - Modified dataset with Geo Transform written. + Modified dataset with transform written. """ transform = transform or self.transform(recalc=True) data_obj = self._get_obj(inplace=inplace) - # delete the old attribute to prevent confusion - data_obj.attrs.pop("transform", None) - grid_mapping_name = ( - self.grid_mapping if grid_mapping_name is None else grid_mapping_name - ) - try: - grid_map_attrs = data_obj.coords[grid_mapping_name].attrs.copy() - except KeyError: - data_obj.coords[grid_mapping_name] = xarray.Variable((), 0) - grid_map_attrs = data_obj.coords[grid_mapping_name].attrs.copy() - grid_map_attrs["GeoTransform"] = " ".join( - [str(item) for item in transform.to_gdal()] - ) - data_obj.coords[grid_mapping_name].rio.set_attrs(grid_map_attrs, inplace=True) - return data_obj.rio.write_grid_mapping( - grid_mapping_name=grid_mapping_name, inplace=True - ) + + # Determine which convention to use + if convention is None: + convention = get_option(CONVENTION) or Convention.CF + + if convention is Convention.CF: + grid_mapping_name = ( + self.grid_mapping if grid_mapping_name is None else grid_mapping_name + ) + data_obj = cf.write_transform( + data_obj, + transform, + grid_mapping_name, + inplace=True, + ) + return data_obj.rio.write_grid_mapping( + grid_mapping_name=grid_mapping_name, inplace=True + ) + else: + raise ValueError(f"Unsupported convention: {convention}") def transform(self, recalc: bool = False) -> Affine: """ diff --git a/test/integration/test_integration_zarr_conventions.py b/test/integration/test_integration_zarr_conventions.py new file mode 100644 index 00000000..b1a45a42 --- /dev/null +++ b/test/integration/test_integration_zarr_conventions.py @@ -0,0 +1,214 @@ +"""Integration tests for reading Zarr conventions.""" +import numpy as np +import pyproj +import pytest +import xarray as xr +from affine import Affine +from rasterio.crs import CRS + +import rioxarray +from rioxarray import set_options +from rioxarray._convention import zarr +from rioxarray.enum import Convention + + +def _create_zarr_array_with_proj(): + """Create a DataArray with Zarr proj: convention attributes.""" + data = xr.DataArray( + np.random.rand(10, 20), + dims=["y", "x"], + coords={ + "y": np.arange(10), + "x": np.arange(20), + }, + ) + data.attrs["zarr_conventions"] = [zarr.PROJ_CONVENTION] + data.attrs["proj:wkt2"] = CRS.from_epsg(4326).to_wkt() + return data + + +def _create_zarr_array_with_spatial(): + """Create a DataArray with Zarr spatial: convention attributes.""" + data = xr.DataArray( + np.random.rand(10, 20), + dims=["lat", "lon"], + coords={ + "lat": np.arange(10), + "lon": np.arange(20), + }, + ) + data.attrs["zarr_conventions"] = [zarr.SPATIAL_CONVENTION] + data.attrs["spatial:transform"] = [1.0, 0.0, 100.0, 0.0, -1.0, 200.0] + data.attrs["spatial:dimensions"] = ["lat", "lon"] + return data + + +def _create_zarr_array_with_both(): + """Create a DataArray with both Zarr conventions.""" + data = xr.DataArray( + np.random.rand(10, 20), + dims=["lat", "lon"], + coords={ + "lat": np.arange(10), + "lon": np.arange(20), + }, + ) + data.attrs["zarr_conventions"] = [zarr.PROJ_CONVENTION, zarr.SPATIAL_CONVENTION] + data.attrs["proj:wkt2"] = CRS.from_epsg(32618).to_wkt() + data.attrs["spatial:transform"] = [10.0, 0.0, 500000.0, 0.0, -10.0, 4500000.0] + data.attrs["spatial:dimensions"] = ["lat", "lon"] + return data + + +def test_read_crs_from_zarr_convention(): + """Test reading CRS from DataArray with Zarr proj: convention.""" + data = _create_zarr_array_with_proj() + + crs = data.rio.crs + assert crs is not None + assert crs == CRS.from_epsg(4326) + + +def test_read_crs_from_zarr_convention__with_setting(): + """Test reading CRS with Convention.Zarr setting.""" + data = _create_zarr_array_with_proj() + + with set_options(convention=Convention.Zarr): + crs = data.rio.crs + assert crs is not None + assert crs == CRS.from_epsg(4326) + + +def test_read_transform_from_zarr_convention(): + """Test reading transform from DataArray with Zarr spatial: convention.""" + data = _create_zarr_array_with_spatial() + + # Access transform via rio accessor + transform = data.rio.transform() + # Transform should be calculated from coordinates since we didn't + # set up the full spatial metadata, but let's check the cached version + cached = data.rio._cached_transform() + assert cached is not None + assert cached == Affine(1.0, 0.0, 100.0, 0.0, -1.0, 200.0) + + +def test_read_spatial_dimensions_from_zarr_convention(): + """Test reading spatial dimensions from Zarr spatial: convention.""" + data = _create_zarr_array_with_spatial() + + assert data.rio.x_dim == "lon" + assert data.rio.y_dim == "lat" + + +def test_read_both_conventions(): + """Test reading from DataArray with both Zarr conventions.""" + data = _create_zarr_array_with_both() + + # CRS from proj: + crs = data.rio.crs + assert crs is not None + assert crs == CRS.from_epsg(32618) + + # Transform from spatial: + cached = data.rio._cached_transform() + assert cached is not None + assert cached == Affine(10.0, 0.0, 500000.0, 0.0, -10.0, 4500000.0) + + # Dimensions from spatial: + assert data.rio.x_dim == "lon" + assert data.rio.y_dim == "lat" + + +def test_fallback_zarr_to_cf(): + """Test that CF convention is tried as fallback when Zarr not found.""" + # Create data with CF convention + data = xr.DataArray( + np.random.rand(10, 20), + dims=["y", "x"], + coords={ + "y": np.arange(10), + "x": np.arange(20), + }, + ) + data.coords["spatial_ref"] = xr.Variable((), 0) + data.coords["spatial_ref"].attrs["spatial_ref"] = "EPSG:4326" + + # Even with Zarr preference, should fall back to CF + with set_options(convention=Convention.Zarr): + crs = data.rio.crs + assert crs is not None + assert crs == CRS.from_epsg(4326) + + +def test_fallback_cf_to_zarr(): + """Test that Zarr convention is tried as fallback when CF not found.""" + # Create data with Zarr convention only + data = _create_zarr_array_with_proj() + + # With CF preference (default), should fall back to Zarr + crs = data.rio.crs + assert crs is not None + assert crs == CRS.from_epsg(4326) + + +def test_priority_zarr_over_cf(): + """Test that Zarr convention takes priority when setting is Zarr.""" + # Create data with both conventions (different CRS values) + data = xr.DataArray( + np.random.rand(10, 20), + dims=["y", "x"], + coords={ + "y": np.arange(10), + "x": np.arange(20), + }, + ) + # CF convention + data.coords["spatial_ref"] = xr.Variable((), 0) + data.coords["spatial_ref"].attrs["spatial_ref"] = "EPSG:4326" + + # Zarr convention (different CRS) + data.attrs["zarr_conventions"] = [zarr.PROJ_CONVENTION] + data.attrs["proj:wkt2"] = CRS.from_epsg(32618).to_wkt() + + # With Zarr setting, should prefer Zarr CRS + with set_options(convention=Convention.Zarr): + crs = data.rio.crs + assert crs is not None + assert crs == CRS.from_epsg(32618) + + # Reset to check default + data2 = data.copy(deep=True) + data2.rio._crs = None # Reset cached CRS + + # With default setting (CF priority), should prefer CF CRS + crs = data2.rio.crs + assert crs is not None + assert crs == CRS.from_epsg(4326) + + +def test_read_proj_code(): + """Test reading CRS from proj:code attribute.""" + data = xr.DataArray( + np.random.rand(10, 20), + dims=["y", "x"], + ) + data.attrs["zarr_conventions"] = [zarr.PROJ_CONVENTION] + data.attrs["proj:code"] = "EPSG:32618" + + crs = data.rio.crs + assert crs is not None + assert crs == CRS.from_epsg(32618) + + +def test_read_proj_projjson(): + """Test reading CRS from proj:projjson attribute.""" + data = xr.DataArray( + np.random.rand(10, 20), + dims=["y", "x"], + ) + data.attrs["zarr_conventions"] = [zarr.PROJ_CONVENTION] + data.attrs["proj:projjson"] = pyproj.CRS.from_epsg(4326).to_json_dict() + + crs = data.rio.crs + assert crs is not None + assert crs == CRS.from_epsg(4326) diff --git a/test/unit/test_convention_cf.py b/test/unit/test_convention_cf.py new file mode 100644 index 00000000..99a79660 --- /dev/null +++ b/test/unit/test_convention_cf.py @@ -0,0 +1,179 @@ +"""Unit tests for the CF convention module.""" +import numpy as np +import xarray as xr +from affine import Affine +from rasterio.crs import CRS + +from rioxarray._convention import cf + + +def test_read_crs__from_grid_mapping_spatial_ref(): + """Test reading CRS from grid_mapping coordinate's spatial_ref attribute.""" + data = xr.DataArray(np.random.rand(10, 10), dims=["y", "x"]) + data.coords["spatial_ref"] = xr.Variable((), 0) + data.coords["spatial_ref"].attrs["spatial_ref"] = "EPSG:4326" + + crs = cf.read_crs(data, "spatial_ref") + assert crs is not None + assert crs == CRS.from_epsg(4326) + + +def test_read_crs__from_grid_mapping_crs_wkt(): + """Test reading CRS from grid_mapping coordinate's crs_wkt attribute.""" + data = xr.DataArray(np.random.rand(10, 10), dims=["y", "x"]) + data.coords["spatial_ref"] = xr.Variable((), 0) + data.coords["spatial_ref"].attrs["crs_wkt"] = CRS.from_epsg(4326).to_wkt() + + crs = cf.read_crs(data, "spatial_ref") + assert crs is not None + assert crs == CRS.from_epsg(4326) + + +def test_read_crs__from_attrs(): + """Test reading CRS from object's attrs.""" + data = xr.DataArray(np.random.rand(10, 10), dims=["y", "x"]) + data.attrs["crs"] = "EPSG:4326" + + crs = cf.read_crs(data) + assert crs is not None + assert crs == CRS.from_epsg(4326) + + +def test_read_crs__from_attrs_with_missing_grid_mapping(): + """Test reading CRS from attrs when grid_mapping is specified but doesn't exist. + + This tests a common case where rioxarray's grid_mapping property returns + "spatial_ref" as a default, but the coordinate doesn't actually exist. + The CRS should still be found in the object's attrs. + """ + data = xr.DataArray(np.random.rand(10, 10), dims=["y", "x"]) + data.attrs["crs"] = "EPSG:4326" + + # Pass grid_mapping that doesn't exist as a coordinate + crs = cf.read_crs(data, "spatial_ref") + assert crs is not None + assert crs == CRS.from_epsg(4326) + + +def test_read_crs__not_found(): + """Test that None is returned when no CRS is found.""" + data = xr.DataArray(np.random.rand(10, 10), dims=["y", "x"]) + + crs = cf.read_crs(data) + assert crs is None + + +def test_read_transform__from_geotransform(): + """Test reading transform from GeoTransform attribute.""" + data = xr.DataArray(np.random.rand(10, 10), dims=["y", "x"]) + data.coords["spatial_ref"] = xr.Variable((), 0) + # GeoTransform format: [c, a, b, f, d, e] (GDAL format) + data.coords["spatial_ref"].attrs["GeoTransform"] = "0.0 1.0 0.0 10.0 0.0 -1.0" + + transform = cf.read_transform(data, "spatial_ref") + assert transform is not None + assert transform == Affine(1.0, 0.0, 0.0, 0.0, -1.0, 10.0) + + +def test_read_transform__from_attrs(): + """Test reading transform from object's attrs.""" + data = xr.DataArray(np.random.rand(10, 10), dims=["y", "x"]) + # Transform stored as list in attrs + data.attrs["transform"] = [1.0, 0.0, 0.0, 0.0, -1.0, 10.0] + + transform = cf.read_transform(data) + assert transform is not None + assert transform == Affine(1.0, 0.0, 0.0, 0.0, -1.0, 10.0) + + +def test_read_transform__not_found(): + """Test that None is returned when no transform is found.""" + data = xr.DataArray(np.random.rand(10, 10), dims=["y", "x"]) + + transform = cf.read_transform(data) + assert transform is None + + +def test_read_spatial_dimensions__xy(): + """Test detecting x/y dimension names.""" + data = xr.DataArray(np.random.rand(10, 10), dims=["y", "x"]) + + dims = cf.read_spatial_dimensions(data) + assert dims == ("y", "x") + + +def test_read_spatial_dimensions__lonlat(): + """Test detecting longitude/latitude dimension names.""" + data = xr.DataArray(np.random.rand(10, 10), dims=["latitude", "longitude"]) + + dims = cf.read_spatial_dimensions(data) + assert dims == ("latitude", "longitude") + + +def test_read_spatial_dimensions__cf_axis(): + """Test detecting dimensions from CF axis attributes.""" + data = xr.DataArray( + np.random.rand(10, 10), + dims=["row", "col"], + coords={ + "row": ("row", np.arange(10)), + "col": ("col", np.arange(10)), + }, + ) + data.coords["row"].attrs["axis"] = "Y" + data.coords["col"].attrs["axis"] = "X" + + dims = cf.read_spatial_dimensions(data) + assert dims == ("row", "col") + + +def test_read_spatial_dimensions__cf_standard_name(): + """Test detecting dimensions from CF standard_name attributes.""" + data = xr.DataArray( + np.random.rand(10, 10), + dims=["lat", "lon"], + coords={ + "lat": ("lat", np.arange(10)), + "lon": ("lon", np.arange(10)), + }, + ) + data.coords["lat"].attrs["standard_name"] = "latitude" + data.coords["lon"].attrs["standard_name"] = "longitude" + + dims = cf.read_spatial_dimensions(data) + assert dims == ("lat", "lon") + + +def test_read_spatial_dimensions__not_found(): + """Test that None is returned when spatial dimensions are not found.""" + data = xr.DataArray(np.random.rand(10, 10), dims=["a", "b"]) + + dims = cf.read_spatial_dimensions(data) + assert dims is None + + +def test_write_crs(): + """Test writing CRS to a DataArray.""" + data = xr.DataArray(np.random.rand(10, 10), dims=["y", "x"]) + crs = CRS.from_epsg(4326) + + result = cf.write_crs(data, crs, "spatial_ref", inplace=False) + + assert "spatial_ref" in result.coords + assert result.coords["spatial_ref"].attrs["spatial_ref"] == crs.to_wkt() + assert result.coords["spatial_ref"].attrs["crs_wkt"] == crs.to_wkt() + + +def test_write_transform(): + """Test writing transform to a DataArray.""" + data = xr.DataArray(np.random.rand(10, 10), dims=["y", "x"]) + transform = Affine(1.0, 0.0, 0.0, 0.0, -1.0, 10.0) + + result = cf.write_transform(data, transform, "spatial_ref", inplace=False) + + assert "spatial_ref" in result.coords + assert "GeoTransform" in result.coords["spatial_ref"].attrs + assert ( + result.coords["spatial_ref"].attrs["GeoTransform"] + == "0.0 1.0 0.0 10.0 0.0 -1.0" + ) diff --git a/test/unit/test_convention_zarr.py b/test/unit/test_convention_zarr.py new file mode 100644 index 00000000..65b5a8e9 --- /dev/null +++ b/test/unit/test_convention_zarr.py @@ -0,0 +1,203 @@ +"""Unit tests for the Zarr convention module.""" +import json + +import numpy as np +import pytest +import xarray as xr +from affine import Affine +from rasterio.crs import CRS + +from rioxarray._convention import zarr + + +def test_has_convention_declared__proj(): + """Test checking for proj: convention declaration.""" + attrs = { + "zarr_conventions": [ + { + "name": "proj:", + "uuid": "f17cb550-5864-4468-aeb7-f3180cfb622f", + } + ] + } + assert zarr.has_convention_declared(attrs, "proj:") is True + assert zarr.has_convention_declared(attrs, "spatial:") is False + + +def test_has_convention_declared__spatial(): + """Test checking for spatial: convention declaration.""" + attrs = { + "zarr_conventions": [ + { + "name": "spatial:", + "uuid": "689b58e2-cf7b-45e0-9fff-9cfc0883d6b4", + } + ] + } + assert zarr.has_convention_declared(attrs, "spatial:") is True + assert zarr.has_convention_declared(attrs, "proj:") is False + + +def test_has_convention_declared__not_declared(): + """Test when no convention is declared.""" + attrs = {} + assert zarr.has_convention_declared(attrs, "proj:") is False + assert zarr.has_convention_declared(attrs, "spatial:") is False + + +def test_get_declared_conventions(): + """Test getting all declared conventions.""" + attrs = { + "zarr_conventions": [ + {"name": "proj:", "uuid": "test-uuid-1"}, + {"name": "spatial:", "uuid": "test-uuid-2"}, + ] + } + declared = zarr.get_declared_conventions(attrs) + assert declared == {"proj:", "spatial:"} + + +def test_parse_spatial_transform(): + """Test parsing spatial:transform array.""" + transform_array = [1.0, 0.0, 100.0, 0.0, -1.0, 200.0] + result = zarr.parse_spatial_transform(transform_array) + assert result == Affine(1.0, 0.0, 100.0, 0.0, -1.0, 200.0) + + +def test_parse_spatial_transform__invalid(): + """Test parsing invalid spatial:transform.""" + assert zarr.parse_spatial_transform([1, 2, 3]) is None + assert zarr.parse_spatial_transform("invalid") is None + + +def test_parse_proj_code(): + """Test parsing proj:code.""" + result = zarr.parse_proj_code("EPSG:4326") + assert result is not None + assert result == CRS.from_epsg(4326) + + +def test_parse_proj_wkt2(): + """Test parsing proj:wkt2.""" + wkt = CRS.from_epsg(4326).to_wkt() + result = zarr.parse_proj_wkt2(wkt) + assert result is not None + assert result == CRS.from_epsg(4326) + + +def test_parse_proj_projjson__dict(): + """Test parsing proj:projjson from dict.""" + # Use pyproj to get PROJJSON since rasterio CRS doesn't have to_json + import pyproj + + projjson = pyproj.CRS.from_epsg(4326).to_json_dict() + result = zarr.parse_proj_projjson(projjson) + assert result is not None + assert result == CRS.from_epsg(4326) + + +def test_parse_proj_projjson__string(): + """Test parsing proj:projjson from JSON string.""" + import pyproj + + projjson_str = pyproj.CRS.from_epsg(4326).to_json() + result = zarr.parse_proj_projjson(projjson_str) + assert result is not None + assert result == CRS.from_epsg(4326) + + +def test_read_crs__from_wkt2(): + """Test reading CRS from proj:wkt2 attribute.""" + data = xr.DataArray(np.random.rand(10, 10), dims=["y", "x"]) + data.attrs["zarr_conventions"] = [zarr.PROJ_CONVENTION] + data.attrs["proj:wkt2"] = CRS.from_epsg(4326).to_wkt() + + crs = zarr.read_crs(data) + assert crs is not None + assert crs == CRS.from_epsg(4326) + + +def test_read_crs__from_code(): + """Test reading CRS from proj:code attribute.""" + data = xr.DataArray(np.random.rand(10, 10), dims=["y", "x"]) + data.attrs["zarr_conventions"] = [zarr.PROJ_CONVENTION] + data.attrs["proj:code"] = "EPSG:4326" + + crs = zarr.read_crs(data) + assert crs is not None + assert crs == CRS.from_epsg(4326) + + +def test_read_crs__not_found(): + """Test that None is returned when no CRS is found.""" + data = xr.DataArray(np.random.rand(10, 10), dims=["y", "x"]) + + crs = zarr.read_crs(data) + assert crs is None + + +def test_read_crs__no_convention_declared(): + """Test that CRS is not read when convention is not declared.""" + data = xr.DataArray(np.random.rand(10, 10), dims=["y", "x"]) + # Add proj attributes but no convention declaration + data.attrs["proj:wkt2"] = CRS.from_epsg(4326).to_wkt() + + crs = zarr.read_crs(data) + assert crs is None + + +def test_read_transform__from_spatial_transform(): + """Test reading transform from spatial:transform attribute.""" + data = xr.DataArray(np.random.rand(10, 10), dims=["y", "x"]) + data.attrs["zarr_conventions"] = [zarr.SPATIAL_CONVENTION] + data.attrs["spatial:transform"] = [1.0, 0.0, 100.0, 0.0, -1.0, 200.0] + + transform = zarr.read_transform(data) + assert transform is not None + assert transform == Affine(1.0, 0.0, 100.0, 0.0, -1.0, 200.0) + + +def test_read_transform__not_found(): + """Test that None is returned when no transform is found.""" + data = xr.DataArray(np.random.rand(10, 10), dims=["y", "x"]) + + transform = zarr.read_transform(data) + assert transform is None + + +def test_read_transform__no_convention_declared(): + """Test that transform is not read when convention is not declared.""" + data = xr.DataArray(np.random.rand(10, 10), dims=["y", "x"]) + # Add spatial attributes but no convention declaration + data.attrs["spatial:transform"] = [1.0, 0.0, 100.0, 0.0, -1.0, 200.0] + + transform = zarr.read_transform(data) + assert transform is None + + +def test_read_spatial_dimensions(): + """Test reading spatial dimensions from spatial:dimensions attribute.""" + data = xr.DataArray(np.random.rand(10, 20), dims=["lat", "lon"]) + data.attrs["zarr_conventions"] = [zarr.SPATIAL_CONVENTION] + data.attrs["spatial:dimensions"] = ["lat", "lon"] + + dims = zarr.read_spatial_dimensions(data) + assert dims == ("lat", "lon") + + +def test_read_spatial_dimensions__not_found(): + """Test that None is returned when no spatial dimensions are found.""" + data = xr.DataArray(np.random.rand(10, 10), dims=["y", "x"]) + + dims = zarr.read_spatial_dimensions(data) + assert dims is None + + +def test_read_spatial_dimensions__no_convention_declared(): + """Test that spatial dims are not read when convention is not declared.""" + data = xr.DataArray(np.random.rand(10, 10), dims=["y", "x"]) + # Add spatial attributes but no convention declaration + data.attrs["spatial:dimensions"] = ["y", "x"] + + dims = zarr.read_spatial_dimensions(data) + assert dims is None diff --git a/test/unit/test_options.py b/test/unit/test_options.py index 4b7388f5..3de7e653 100644 --- a/test/unit/test_options.py +++ b/test/unit/test_options.py @@ -1,7 +1,8 @@ import pytest from rioxarray import set_options -from rioxarray._options import EXPORT_GRID_MAPPING, get_option +from rioxarray._options import CONVENTION, EXPORT_GRID_MAPPING, get_option +from rioxarray.enum import Convention def test_set_options__contextmanager(): @@ -37,3 +38,43 @@ def test_set_options__invalid_value(): ): with set_options(export_grid_mapping=12345): pass + + +def test_set_options__convention_default(): + """Test that convention defaults to None.""" + assert get_option(CONVENTION) is None + + +def test_set_options__convention_cf(): + """Test setting convention to CF.""" + assert get_option(CONVENTION) is None + with set_options(convention=Convention.CF): + assert get_option(CONVENTION) is Convention.CF + assert get_option(CONVENTION) is None + + +def test_set_options__convention_zarr(): + """Test setting convention to Zarr.""" + assert get_option(CONVENTION) is None + with set_options(convention=Convention.Zarr): + assert get_option(CONVENTION) is Convention.Zarr + assert get_option(CONVENTION) is None + + +def test_set_options__convention_none(): + """Test setting convention back to None.""" + with set_options(convention=Convention.CF): + assert get_option(CONVENTION) is Convention.CF + with set_options(convention=None): + assert get_option(CONVENTION) is None + assert get_option(CONVENTION) is Convention.CF + + +def test_set_options__convention_invalid(): + """Test that invalid convention values raise error.""" + with pytest.raises( + ValueError, + match="option 'convention' gave an invalid value: 'invalid'.", + ): + with set_options(convention="invalid"): + pass