Build STAC Metadata For A COG

from pathlib import Path


def find_repo_root(start: Path = Path.cwd()) -> Path:
    """Find the repository root when the kernel starts in a subfolder."""
    start = start.resolve()
    for path in (start, *start.parents):
        if (path / ".git").exists() or (path / "downloaded_data").exists():
            return path
    return start


REPO_ROOT = find_repo_root()
DATA_DIR = REPO_ROOT / "downloaded_data"
DATA_DIR.mkdir(exist_ok=True)

print(f"Repository root: {REPO_ROOT}")
print(f"Data directory: {DATA_DIR}")

import datetime as dt
import json

import pystac
import rasterio
from pystac.extensions.projection import ProjectionExtension
from pystac.extensions.raster import RasterBand, RasterExtension
from rasterio.warp import transform_bounds
from shapely.geometry import box, mapping

COG_PATH = DATA_DIR / "basal_melt_map_racmo_firn_air_corrected_cog.tif"
PUBLIC_HREF = "https://s3.waw4-1.cloudferro.com/EarthCODE/OSCAssets/polar_cube_datasets/basal_melt/basal_melt_map_racmo_firn_air_corrected.tif"

if not COG_PATH.exists():
    raise FileNotFoundError("Create the COG with 1_cog_from_tiff.ipynb before running this notebook.")

with rasterio.open(COG_PATH) as src:
    bounds_native = src.bounds
    bbox_wgs84 = transform_bounds(src.crs, "EPSG:4326", *bounds_native, densify_pts=21)
    geometry_wgs84 = mapping(box(*bbox_wgs84))
    epsg = src.crs.to_epsg()
    transform = [src.transform.a, src.transform.b, src.transform.c, src.transform.d, src.transform.e, src.transform.f]
    shape = [src.height, src.width]
    data_type = src.dtypes[0]
    nodata = src.nodata
    resolution = abs(src.res[0])

print(f"EPSG: {epsg}")
print(f"WGS84 bbox: {bbox_wgs84}")
print(f"Shape: {shape}")
print(f"Data type: {data_type}, nodata: {nodata}")

item = pystac.Item(
    id="antarctic-ice-shelf-basal-melt-cog",
    geometry=geometry_wgs84,
    bbox=list(bbox_wgs84),
    datetime=dt.datetime(2021, 1, 1, tzinfo=dt.timezone.utc),
    properties={
        "title": "Antarctic ice shelf basal melt rate COG",
        "description": "Cloud Optimized GeoTIFF example for Antarctic ice shelf basal melt rates.",
        "license": "cc-by-4.0",
        "start_datetime": "1997-01-01T00:00:00Z",
        "end_datetime": "2021-12-31T23:59:59Z",
    },
)

ProjectionExtension.add_to(item)
proj = ProjectionExtension.ext(item)
proj.epsg = epsg
proj.shape = shape
proj.bbox = list(bounds_native)
proj.transform = transform

asset = pystac.Asset(
    href=PUBLIC_HREF,
    media_type=pystac.MediaType.COG,
    roles=["data"],
    title="Basal melt rate COG",
)
RasterExtension.add_to(item)
item.add_asset("basal_melt_rate", asset)

raster = RasterExtension.ext(item.assets["basal_melt_rate"])
raster.bands = [
    RasterBand.create(
        nodata=nodata,
        data_type=data_type,
        spatial_resolution=resolution,
        unit="m yr-1",
    )
]

item.validate()
item.to_dict()