From d5a4c6d36ed4602823ff6d6bf2d7fc699487540f Mon Sep 17 00:00:00 2001 From: Deepak Cherian Date: Wed, 26 Jun 2024 16:10:23 -0400 Subject: [PATCH] Detect grid_mapping in encode_geometries --- cf_xarray/geometry.py | 58 +++++++++++++++++++++++++++----- cf_xarray/tests/test_geometry.py | 14 ++++---- 2 files changed, 56 insertions(+), 16 deletions(-) diff --git a/cf_xarray/geometry.py b/cf_xarray/geometry.py index 760436f7..febac8bc 100644 --- a/cf_xarray/geometry.py +++ b/cf_xarray/geometry.py @@ -17,6 +17,19 @@ ] +# Useful convention language: +# 1. Whether linked to normal CF space-time coordinates with a nodes attribute or not, inclusion of such coordinates is +# recommended to maintain backward compatibility with software that has not implemented geometry capabilities. +# 2. The geometry node coordinate variables must each have an axis attribute whose allowable values are X, Y, and Z. +# 3. If a coordinates attribute is carried by the geometry container variable or its parent data variable, then those coordinate variables +# that have a meaningful correspondence with node coordinates are indicated as such by a nodes attribute that names the corresponding node +# coordinates, but only if the grid_mapping associated the geometry node variables is the same as that of the coordinate variables. +# If a different grid mapping is used, then the provided coordinates must not have the nodes attribute. +# +# Interpretation: +# 1. node coordinates are exact; the 'normal' coordinates are a reasonable value to use, if you do not know how to interpret the nodes. + + def decode_geometries(encoded: xr.Dataset) -> xr.Dataset: """ Decode CF encoded geometries to a numpy object array containing shapely geometries. @@ -282,6 +295,14 @@ def shapely_to_cf(geometries: xr.DataArray | Sequence, grid_mapping: str | None ---------- Please refer to the CF conventions document: http://cfconventions.org/Data/cf-conventions/cf-conventions-1.8/cf-conventions.html#geometries """ + + if isinstance(geometries, xr.DataArray) and grid_mapping is not None: + raise DeprecationWarning( + "Explicitly passing `grid_mapping` with DataArray of geometries is deprecated. " + "Please set a `grid_mapping` attribute on `geometries`, ", + "and set the grid mapping variable as a coordinate", + ) + # Get all types to call the appropriate translation function. types = { geom.item().geom_type if isinstance(geom, xr.DataArray) else geom.geom_type @@ -300,19 +321,38 @@ def shapely_to_cf(geometries: xr.DataArray | Sequence, grid_mapping: str | None ds[GEOMETRY_CONTAINER_NAME].attrs.update(coordinates="crd_x crd_y") + if ( + grid_mapping is None + and isinstance(geometries, xr.DataArray) + and (grid_mapping_varname := geometries.attrs.get("grid_mapping")) + ): + if grid_mapping_varname in geometries.coords: + grid_mapping = geometries.coords[grid_mapping_varname].attrs[ + "grid_mapping_name" + ] + for name_ in ["x", "y", "crd_x", "crd_y"]: + ds[name_].attrs["grid_mapping"] = grid_mapping_varname + # Special treatment of selected grid mappings - if grid_mapping == "longitude_latitude": - # Special case for longitude_latitude grid mapping + if grid_mapping in ["latitude_longitude", "rotated_latitude_longitude"]: + # Special case for longitude_latitude type grid mappings ds = ds.rename(crd_x="lon", crd_y="lat") - ds.lon.attrs.update(units="degrees_east", standard_name="longitude") - ds.lat.attrs.update(units="degrees_north", standard_name="latitude") + if grid_mapping == "latitude_longitude": + ds.lon.attrs.update(units="degrees_east", standard_name="longitude") + ds.x.attrs.update(units="degrees_east", standard_name="longitude") + ds.lat.attrs.update(units="degrees_north", standard_name="latitude") + ds.y.attrs.update(units="degrees_north", standard_name="latitude") + elif grid_mapping == "rotated_latitude_longitude": + ds.lon.attrs.update(units="degrees", standard_name="grid_longitude") + ds.x.attrs.update(units="degrees", standard_name="grid_longitude") + ds.lat.attrs.update(units="degrees", standard_name="grid_latitude") + ds.y.attrs.update(units="degrees", standard_name="grid_latitude") ds[GEOMETRY_CONTAINER_NAME].attrs.update(coordinates="lon lat") - ds.x.attrs.update(units="degrees_east", standard_name="longitude") - ds.y.attrs.update(units="degrees_north", standard_name="latitude") elif grid_mapping is not None: - raise NotImplementedError( - f"Only grid mapping longitude_latitude is implemented. Got {grid_mapping}." - ) + ds.crd_x.attrs.update(standard_name="projection_x_coordinate") + ds.x.attrs.update(standard_name="projection_x_coordinate") + ds.crd_y.attrs.update(standard_name="projection_y_coordinate") + ds.y.attrs.update(standard_name="projection_y_coordinate") return ds diff --git a/cf_xarray/tests/test_geometry.py b/cf_xarray/tests/test_geometry.py index b026f23e..c412f3c7 100644 --- a/cf_xarray/tests/test_geometry.py +++ b/cf_xarray/tests/test_geometry.py @@ -265,8 +265,8 @@ def test_shapely_to_cf(geometry_ds): [ in_ds.drop_vars("geometry").isel(index=slice(1, None)), cfxr.shapely_to_cf( - in_ds.geometry.isel(index=slice(1, None)), - grid_mapping="longitude_latitude", + in_ds.geometry.isel(index=slice(1, None)).data, + grid_mapping="latitude_longitude", ), ] ) @@ -355,10 +355,11 @@ def test_shapely_to_cf_errors(): with pytest.raises(ValueError, match="Mixed geometry types are not supported"): cfxr.shapely_to_cf(geoms) - with pytest.raises( - NotImplementedError, match="Only grid mapping longitude_latitude" - ): - cfxr.shapely_to_cf([Point(4, 5)], grid_mapping="albers_conical_equal_area") + encoded = cfxr.shapely_to_cf( + [Point(4, 5)], grid_mapping="albers_conical_equal_area" + ) + assert encoded["x"].attrs["standard_name"] == "projection_x_coordinate" + assert encoded["y"].attrs["standard_name"] == "projection_y_coordinate" @requires_shapely @@ -491,7 +492,6 @@ def test_reshape_unique_geometries(geometry_ds): @requires_shapely def test_encode_decode(geometry_ds, polygon_geometry): - geom_dim_ds = xr.Dataset() geom_dim_ds = geom_dim_ds.assign_coords( xr.Coordinates(