Delete GeoTransform in assign_index

Closes #28

Author: dcherian

docs/raster_index/design_choices.md

@@ -6,3 +6,14 @@
 
 In designing {py:class}`RasterIndex`, we faced a few thorny questions. Below we discuss these considerations, and the approach we've taken.
 Ultimately, there are no easy answers and tradeoffs to be made.
+
+## Handling the `GeoTransform` attribute
+
+GDAL _chooses_ to track the affine transform using a `GeoTransform` attribute on a `spatial_ref` variable. The `"spatial_ref"` is a
+"grid mapping" variable (as termed by the CF-conventions). It also records CRS information. Currently, our design is that
+{py:class}`xproj.CRSIndex` controls the CRS information and handles the creation of the `"spatial_ref"` variable, or more generally,
+the grid mapping variable. Thus, {py:class}`RasterIndex` _cannot_ keep the `"GeoTransform"` attribute on `"spatial_ref"` up-to-date
+because it does not control it.
+
+Thus, {py:func}`assign_index` will delete the `"GeoTransform"` attribute on the grid mapping variable if it is detected, after using it
+to construct the affine matrix.

src/rasterix/utils.py

@@ -24,7 +24,9 @@ def get_grid_mapping_var(obj: xr.Dataset | xr.DataArray) -> xr.DataArray | None:
     return None
 
 
-def get_affine(obj: xr.Dataset | xr.DataArray, *, x_dim="x", y_dim="y") -> Affine:
+def get_affine(
+    obj: xr.Dataset | xr.DataArray, *, x_dim="x", y_dim="y", clear_transform: bool = False
+) -> Affine:
     """
     Grabs an affine transform from an Xarray object.
 
@@ -39,13 +41,17 @@ def get_affine(obj: xr.Dataset | xr.DataArray, *, x_dim="x", y_dim="y") -> Affin
         Name of the X dimension coordinate variable.
     y_dim: str, optional
         Name of the Y dimension coordinate variable.
+    clear_transform: bool
+       Whether to delete the ``GeoTransform`` attribute if detected.
 
     Returns
     -------
     affine: Affine
     """
     grid_mapping_var = get_grid_mapping_var(obj)
     if grid_mapping_var is not None and (transform := grid_mapping_var.attrs.get("GeoTransform")):
+        if clear_transform:
+            del grid_mapping_var.attrs["GeoTransform"]
         return Affine.from_gdal(*map(float, transform.split(" ")))
     else:
         x = obj.coords[x_dim]

tests/test_raster_index.py

@@ -85,6 +85,7 @@ def test_sel_slice():
     transform = Affine.identity()
     ds.coords["spatial_ref"] = ((), 0, {"GeoTransform": " ".join(map(str, transform.to_gdal()))})
     ds = assign_index(ds)
+    assert "GeoTransform" not in ds.spatial_ref.attrs
     assert ds.xindexes["x"].transform() == transform
 
     actual = ds.sel(x=slice(4), y=slice(3, 5))