fix: embed EPSG in Zarr attrs for TiTiler CRS detection

TiTiler defaults to EPSG:4326 when ds.rio.crs is None, causing validation
errors on UTM coordinates. Extract helper to embed proj:epsg in attrs so
rioxarray reconstructs CRS after Zarr I/O.

Author: wietzesuijker

src/eopf_geozarr/conversion/geozarr.py

@@ -352,6 +352,13 @@ def iterative_copy(
     return dt_result if isinstance(dt_result, xr.DataTree) else xr.DataTree(dt_result)
 
 
+def _embed_crs_in_attrs(ds: xr.Dataset) -> None:
+    """Embed EPSG code in dataset attrs for Zarr persistence (modifies in-place)."""
+    if epsg_code := ds.rio.crs.to_epsg():
+        ds.attrs["proj:epsg"] = epsg_code
+        ds.attrs["proj:code"] = f"EPSG:{epsg_code}"
+
+
 def prepare_dataset_with_crs_info(
     ds: xr.Dataset, reference_crs: str | None = None
 ) -> xr.Dataset:
@@ -380,6 +387,7 @@ def prepare_dataset_with_crs_info(
         print(f"  Adding CRS information: {reference_crs}")
         ds = ds.rio.write_crs(reference_crs)
         ds.attrs["grid_mapping"] = "spatial_ref"
+        _embed_crs_in_attrs(ds)
 
         # Ensure spatial_ref variable has proper attributes
         if "spatial_ref" in ds:
@@ -1003,9 +1011,10 @@ def create_overview_dataset_all_vars(
     # Create overview dataset
     overview_ds = xr.Dataset(overview_data_vars, coords=overview_coords)
 
-    # Set CRS using rioxarray first
+    # Set CRS using rioxarray
     overview_ds.rio.write_crs(native_crs, inplace=True)
     overview_ds.attrs["grid_mapping"] = grid_mapping_var_name
+    _embed_crs_in_attrs(overview_ds)
 
     # Add grid_mapping variable after setting CRS
     # TODO: refactor? grid mapping attributes and variables are handled

test_crs_fix.py

@@ -0,0 +1,117 @@
+#!/usr/bin/env python3
+"""
+Quick validation script to verify CRS metadata is embedded in Zarr attributes.
+Tests the fix for TiTiler 500 error caused by missing CRS.
+"""
+
+import tempfile
+from pathlib import Path
+
+import numpy as np
+import xarray as xr
+
+from eopf_geozarr.conversion.geozarr import prepare_dataset_with_crs_info
+
+
+def test_crs_metadata_embedding() -> None:
+    """Test that prepare_dataset_with_crs_info embeds EPSG in attrs."""
+
+    # Create a minimal test dataset with spatial coordinates
+    x = np.arange(0, 110000, 10)  # 10m resolution, 10km extent
+    y = np.arange(5300000, 5200000, -10)  # UTM coordinates
+    data = np.random.rand(len(y), len(x))
+
+    ds = xr.Dataset(
+        {
+            "reflectance": (
+                ["y", "x"],
+                data,
+                {"standard_name": "toa_bidirectional_reflectance"},
+            ),
+        },
+        coords={
+            "x": (["x"], x, {"standard_name": "projection_x_coordinate", "units": "m"}),
+            "y": (["y"], y, {"standard_name": "projection_y_coordinate", "units": "m"}),
+        },
+    )
+
+    # Apply CRS preparation (this is what happens before to_zarr calls)
+    reference_crs = "EPSG:32632"  # UTM Zone 32N
+    ds_prepared = prepare_dataset_with_crs_info(ds, reference_crs=reference_crs)
+
+    # Drop encoding to avoid conflicts (this is done in real conversion code)
+    ds_prepared = ds_prepared.drop_encoding()
+
+    # Verify CRS metadata is in attrs (this persists to Zarr)
+    assert "proj:epsg" in ds_prepared.attrs, "proj:epsg not found in dataset attrs"
+    assert "proj:code" in ds_prepared.attrs, "proj:code not found in dataset attrs"
+    assert ds_prepared.attrs["proj:epsg"] == 32632, (
+        f"Expected EPSG:32632, got {ds_prepared.attrs['proj:epsg']}"
+    )
+    assert ds_prepared.attrs["proj:code"] == "EPSG:32632", (
+        f"Expected 'EPSG:32632', got {ds_prepared.attrs['proj:code']}"
+    )
+
+    print("✓ CRS metadata correctly embedded in dataset attrs")
+
+    # Test that it persists through Zarr write/read cycle
+    with tempfile.TemporaryDirectory() as tmpdir:
+        zarr_path = Path(tmpdir) / "test.zarr"
+
+        # Write to Zarr v3
+        ds_prepared.to_zarr(
+            zarr_path,
+            mode="w",
+            zarr_format=3,
+            consolidated=True,
+        )
+
+        # Read back and verify
+        ds_read = xr.open_zarr(zarr_path, zarr_format=3, consolidated=True)
+
+        assert "proj:epsg" in ds_read.attrs, "proj:epsg lost after Zarr write/read"
+        assert ds_read.attrs["proj:epsg"] == 32632, (
+            f"EPSG changed after write/read: {ds_read.attrs['proj:epsg']}"
+        )
+        assert ds_read.attrs["proj:code"] == "EPSG:32632", (
+            f"proj:code changed: {ds_read.attrs['proj:code']}"
+        )
+
+        print("✓ CRS metadata persists through Zarr write/read cycle")
+
+        # Verify rioxarray can read the CRS
+        assert ds_read.rio.crs is not None, "rioxarray cannot read CRS from Zarr"
+        assert ds_read.rio.crs.to_epsg() == 32632, (
+            f"Expected EPSG:32632 from rio.crs, got {ds_read.rio.crs.to_epsg()}"
+        )
+
+        print("✓ rioxarray successfully reads CRS from persisted metadata")
+
+        # Check bounds are in UTM meters (not lat/lon)
+        bounds = ds_read.rio.bounds()
+        # UTM coordinates are typically in range [0, 1000000] for easting, [0, 10000000] for northing
+        assert bounds[0] < 200000 and bounds[2] < 200000, (
+            f"X bounds don't look like UTM coordinates: {bounds}"
+        )
+        assert bounds[1] > 5000000 and bounds[3] > 5000000, (
+            f"Y bounds don't look like UTM coordinates: {bounds}"
+        )
+        print(f"✓ Bounds in UTM meters: {bounds}")
+
+
+if __name__ == "__main__":
+    try:
+        test_crs_metadata_embedding()
+        print("\n✅ All CRS metadata tests passed!")
+        print(
+            "The fix ensures TiTiler will read native CRS instead of defaulting to EPSG:4326"
+        )
+    except AssertionError as e:
+        print(f"\n❌ Test failed: {e}")
+        exit(1)
+    except Exception as e:
+        print(f"\n❌ Unexpected error: {e}")
+        import traceback
+
+        traceback.print_exc()
+        exit(1)