Release version 1.13.0 with new feature: added DatacubeService.enrich_raster_metadata() for enriching STAC items with projection and raster band metadata. Updated project version in pyproject.toml and uv.lock. (#65) (#230)

GitOrigin-RevId: 41fe3bdf621592b16ffac4d48127fa3cbba3fd79

Co-authored-by: Iman Shafiei <iman.shafiei@earthdaily.com>

Author: internal-python-client-app[bot]

CHANGELOG.md

@@ -4,6 +4,12 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [1.13.0]
+
+### Added
+
+- `DatacubeService.enrich_raster_metadata()` for client-side enrichment of STAC items with projection (`proj:code`, `proj:transform`, `proj:shape`, `gsd`) and raster band (`raster:bands`) metadata by reading asset files via rasterio
+
 ## [1.12.0]
 
 ### Added

earthdaily/datacube/_datacube_service.py

@@ -4,10 +4,40 @@
 
 from earthdaily.datacube._builder import _deduplicate_items, build_datacube
 from earthdaily.datacube._datacube import Datacube
+from earthdaily.datacube._enrichment import RasterMetadataEnricher
 from earthdaily.datacube.constants import DEFAULT_DTYPE, DEFAULT_ENGINE, DEFAULT_HREF_PATH, DEFAULT_NODATA
 
 
 class DatacubeService:
+    def enrich_raster_metadata(
+        self,
+        item: Item,
+        *,
+        force: bool = False,
+    ) -> Item:
+        """
+        Enrich a STAC item with projection and raster band metadata by reading asset files.
+
+        Opens each GeoTIFF asset with rasterio to extract CRS, transform, shape,
+        and per-band information (data type, nodata, scale, offset), then writes
+        these fields back into the STAC item structure as proj: and raster:bands
+        extension properties.
+
+        Parameters
+        ----------
+        item : Item
+            STAC item to enrich. The input item is not mutated.
+        force : bool, default False
+            If False, skip assets that already carry projection/raster metadata.
+            If True, overwrite existing values.
+
+        Returns
+        -------
+        Item
+            A new item with enriched metadata.
+        """
+        return RasterMetadataEnricher.enrich_item(item, force=force)
+
     def create(
         self,
         items: list[Item],

earthdaily/datacube/_enrichment.py

@@ -0,0 +1,221 @@
+from __future__ import annotations
+
+from typing import Any
+
+import rasterio
+from pystac import Item
+from rasterio.transform import from_gcps
+
+from earthdaily._eds_logging import LoggerConfig
+from earthdaily.datacube._builder import _replace_item_hrefs
+from earthdaily.datacube.constants import DEFAULT_HREF_PATH
+
+logger = LoggerConfig(logger_name=__name__).get_logger()
+
+
+class RasterMetadataEnricher:
+    """Enriches STAC items with projection and raster band metadata by reading asset files via rasterio."""
+
+    @staticmethod
+    def enrich_item(
+        item: Item,
+        *,
+        force: bool = False,
+    ) -> Item:
+        """
+        Enrich a STAC item with projection and raster band metadata.
+
+        Opens each GeoTIFF asset with rasterio to extract CRS, transform, shape,
+        and per-band information, then writes it back into the STAC item structure.
+
+        Parameters
+        ----------
+        item : Item
+            STAC item to enrich. The input item is not mutated.
+        force : bool, default False
+            If False, skip assets that already carry projection/raster metadata.
+            If True, overwrite existing values.
+
+        Returns
+        -------
+        Item
+            A **new** item with enriched metadata.
+        """
+        item_dict = item.to_dict(transform_hrefs=False)
+        resolved = _replace_item_hrefs([Item.from_dict(item_dict)], DEFAULT_HREF_PATH)
+        item_dict = resolved[0].to_dict(transform_hrefs=False)
+
+        RasterMetadataEnricher._enrich_item(item_dict, force)
+
+        return Item.from_dict(item_dict)
+
+    @staticmethod
+    def _enrich_item(item: dict[str, Any], force: bool) -> None:
+        all_gsd_values: list[float] = []
+        all_epsg_values: list[int] = []
+        all_epsg_code_values: list[str] = []
+
+        current_assets: dict[str, dict] = item.get("assets", {})
+        for asset_key, asset in current_assets.items():
+            if not asset.get("type", "").startswith("image/tiff") or "href" not in asset:
+                continue
+
+            has_projection = bool(asset.get("proj:code") or asset.get("proj:epsg") or asset.get("proj:wkt2"))
+            has_raster = bool(asset.get("raster:bands"))
+
+            if force or not has_projection or not has_raster:
+                result = RasterMetadataEnricher._process_asset(asset)
+                if result is None:
+                    continue
+
+                proj_metadata, raster_metadata = result
+                if proj_metadata is None:
+                    logger.warning("CRS retrieval failed for item %s, asset %s.", item.get("id"), asset_key)
+                    continue
+
+                raster_data = None if not force and has_raster else raster_metadata
+                RasterMetadataEnricher._update_asset_metadata(asset, proj_metadata, raster_data, force)
+            else:
+                RasterMetadataEnricher._copy_properties_to_asset(asset, item.get("properties", {}))
+
+            gsd = asset.get("gsd")
+            epsg = asset.get("proj:epsg")
+            epsg_code = asset.get("proj:code")
+            if gsd:
+                all_gsd_values.append(gsd)
+            if epsg:
+                all_epsg_values.append(epsg)
+            if epsg_code:
+                all_epsg_code_values.append(epsg_code)
+
+        RasterMetadataEnricher._update_item_properties(
+            item, all_gsd_values, all_epsg_values, all_epsg_code_values, force
+        )
+
+    @staticmethod
+    def _extract_metadata(dataset: Any) -> tuple[dict | None, dict | None]:
+        crs = dataset.crs
+        transform = dataset.transform
+
+        if crs is None:
+            gcps, gcp_crs = dataset.gcps
+            if gcps and gcp_crs:
+                transform = from_gcps(gcps)
+                crs = gcp_crs
+            else:
+                return None, None
+
+        epsg = crs.to_epsg()
+        if epsg:
+            crs_info = {"proj:code": f"EPSG:{epsg}", "proj:epsg": epsg}
+        else:
+            crs_info = {"proj:wkt2": crs.to_wkt()}
+
+        proj_metadata = {
+            "proj:transform": [
+                float(transform.a),
+                float(transform.b),
+                float(transform.c),
+                float(transform.d),
+                float(transform.e),
+                float(transform.f),
+            ],
+            "proj:shape": [dataset.height, dataset.width],
+            "gsd": abs(transform.a),
+            **crs_info,
+        }
+
+        bands = []
+        for i in range(dataset.count):
+            band_data: dict[str, Any] = {
+                "data_type": str(dataset.dtypes[i]),
+                "spatial_resolution": abs(transform.a),
+            }
+
+            if dataset.nodata is not None:
+                if isinstance(dataset.nodata, (int, float)) or dataset.nodata in ("nan", "inf", "-inf"):
+                    band_data["nodata"] = dataset.nodata
+
+            if hasattr(dataset, "units") and dataset.units and i < len(dataset.units) and dataset.units[i]:
+                band_data["unit"] = dataset.units[i]
+
+            if (
+                hasattr(dataset, "scales")
+                and dataset.scales
+                and i < len(dataset.scales)
+                and dataset.scales[i] is not None
+            ):
+                band_data["scale"] = dataset.scales[i]
+
+            if (
+                hasattr(dataset, "offsets")
+                and dataset.offsets
+                and i < len(dataset.offsets)
+                and dataset.offsets[i] is not None
+            ):
+                band_data["offset"] = dataset.offsets[i]
+
+            bands.append(band_data)
+
+        raster_metadata = {"raster:bands": bands}
+        return proj_metadata, raster_metadata
+
+    @staticmethod
+    def _process_asset(asset: dict[str, Any]) -> tuple[dict | None, dict | None] | None:
+        try:
+            with rasterio.open(asset["href"]) as dataset:
+                return RasterMetadataEnricher._extract_metadata(dataset)
+        except Exception as e:
+            logger.warning("Failed to open asset %s: %s", asset.get("href", "<no href>"), e)
+            return None
+
+    @staticmethod
+    def _update_asset_metadata(
+        asset: dict[str, Any],
+        proj_metadata: dict | None,
+        raster_metadata: dict | None,
+        force: bool,
+    ) -> None:
+        if proj_metadata:
+            for key, value in proj_metadata.items():
+                if key not in asset or force:
+                    asset[key] = value
+
+        if raster_metadata and ("raster:bands" not in asset or force):
+            asset.update(raster_metadata)
+
+    @staticmethod
+    def _copy_properties_to_asset(asset: dict[str, Any], item_properties: dict[str, Any]) -> None:
+        proj_attrs = ["proj:code", "proj:epsg", "proj:wkt2", "proj:transform", "proj:shape", "gsd"]
+        for attr in proj_attrs:
+            if attr in item_properties and attr not in asset:
+                asset[attr] = item_properties[attr]
+
+        if "raster:bands" in item_properties and "raster:bands" not in asset:
+            asset["raster:bands"] = item_properties["raster:bands"]
+
+    @staticmethod
+    def _update_item_properties(
+        item: dict[str, Any],
+        all_gsd_values: list[float],
+        all_epsg_values: list[int],
+        all_epsg_code_values: list[str],
+        force: bool,
+    ) -> None:
+        properties = item.get("properties", {})
+
+        unique_gsd = list(set(filter(None, all_gsd_values)))
+        unique_epsg = list(set(filter(None, all_epsg_values)))
+        unique_code = list(set(filter(None, all_epsg_code_values)))
+
+        if len(unique_gsd) == 1:
+            if "gsd" not in properties or properties["gsd"] is None or force:
+                properties["gsd"] = unique_gsd[0]
+
+        if len(unique_code) == 1:
+            if "proj:code" not in properties or properties["proj:code"] is None or force:
+                properties["proj:code"] = unique_code[0]
+
+        if len(unique_epsg) == 1:
+            if "proj:epsg" not in properties or properties["proj:epsg"] is None or force:
+                properties["proj:epsg"] = unique_epsg[0]

examples/enrich_raster_metadata_example.py

@@ -0,0 +1,81 @@
+"""
+Enrich Raster Metadata Example
+================================
+
+Demonstrates how to enrich STAC items with projection and raster band
+metadata using ``client.datacube.enrich_raster_metadata()``.
+
+This is useful when a collection's STAC items lack ``proj:`` or
+``raster:bands`` fields and you need them for downstream analysis
+(e.g., datacube creation, reprojection, spatial alignment).
+
+The example uses the ``sentinel-2-l2a`` collection.
+
+Requirements:
+- pip install 'earthdaily[platform,datacube]'
+- Set EDS credentials as environment variables
+"""
+
+try:
+    from dotenv import load_dotenv
+
+    load_dotenv()
+except ImportError:
+    pass
+
+from earthdaily import EDSClient, EDSConfig
+
+
+def _print_asset_metadata(item, label):
+    """Print projection and raster metadata for each tiff asset."""
+    print(f"\n--- {label}: {item.id} ---")
+    print(f"  properties.gsd: {item.properties.get('gsd')}")
+    print(f"  properties.proj:code: {item.properties.get('proj:code')}")
+
+    for name, asset in item.assets.items():
+        if not asset.media_type or not asset.media_type.startswith("image/tiff"):
+            continue
+        ef = asset.extra_fields
+        print(f"  asset '{name}':")
+        print(f"    proj:code      = {ef.get('proj:code')}")
+        print(f"    proj:shape     = {ef.get('proj:shape')}")
+        print(f"    proj:transform = {ef.get('proj:transform')}")
+        print(f"    gsd            = {ef.get('gsd')}")
+        bands = ef.get("raster:bands")
+        if bands:
+            print(f"    raster:bands   = {len(bands)} band(s)")
+            for i, b in enumerate(bands):
+                print(f"      [{i}] data_type={b.get('data_type')}, nodata={b.get('nodata')}")
+        else:
+            print("    raster:bands   = None")
+
+
+def main():
+    config = EDSConfig()
+    client = EDSClient(config)
+
+    print("Searching sentinel-2-l2a ...")
+    items = client.platform.search(
+        collections=["sentinel-2-l2a"],
+        max_items=2,
+    )
+    print(f"Found {len(items)} item(s)")
+
+    if not items:
+        print("No items found.")
+        return
+
+    for item in items:
+        _print_asset_metadata(item, "BEFORE enrichment")
+
+    print("\n\nEnriching raster metadata ...")
+    enriched = [client.datacube.enrich_raster_metadata(item) for item in items]
+
+    for item in enriched:
+        _print_asset_metadata(item, "AFTER enrichment")
+
+    print("\nDone.")
+
+
+if __name__ == "__main__":
+    main()

pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name = "earthdaily"
-version = "1.12.0"
+version = "1.13.0"
 description = "EarthDaily Python Client"
 authors = [
     { name = "EarthDaily", email = "info@earthdaily.com" }