第07章 - 静态 STAC 目录构建

import pystac

cat = pystac.Catalog(id="demo", description="demo", catalog_type=pystac.CatalogType.SELF_CONTAINED)
print(cat.catalog_type)
# CatalogType.SELF_CONTAINED

import rasterio
from rasterio.warp import transform_bounds
from rasterio.crs import CRS

def cog_to_props(path: str) -> dict:
    """从 COG 读出 STAC Item 需要的几何与投影字段。"""
    with rasterio.open(path) as src:
        # WGS84 的 bbox 与 geometry（STAC 要求经纬度）
        if src.crs is None:
            raise ValueError(f"{path}: 缺少 CRS，无法生成 STAC Item")
        left, bottom, right, top = transform_bounds(
            src.crs, CRS.from_epsg(4326),
            *src.bounds, densify_pts=21,
        )
        bbox = [left, bottom, right, top]

        # 多边形足迹：COG 是矩形网格，四个角即可
        geometry = {
            "type": "Polygon",
            "coordinates": [[
                [left, bottom], [right, bottom],
                [right, top], [left, top], [left, bottom],
            ]],
        }

        # proj:transform 用 rasterio 的顺序，9 个元素
        t = src.transform
        proj_transform = [
            t.a, t.b, t.c,
            t.d, t.e, t.f,
            0.0, 0.0, 1.0,
        ]

        props = {
            "proj:code": f"EPSG:{src.crs.to_epsg()}" if src.crs.to_epsg() else None,
            "proj:shape": [src.height, src.width],
            "proj:transform": proj_transform,
            "proj:bbox": list(src.bounds),
            "gsd": abs(src.transform.a),
        }
        # 清理 None，避免写出非法字段
        props = {k: v for k, v in props.items() if v is not None}
        return {"bbox": bbox, "geometry": geometry, "properties": props}

import datetime as dt
import pystac

def make_collection(collection_id: str, title: str, items: list[pystac.Item]) -> pystac.Collection:
    """用一批 Item 反算 extent，构造 Collection。"""
    bboxes = [it.bbox for it in items if it.bbox]
    west = min(b[0] for b in bboxes)
    south = min(b[1] for b in bboxes)
    east = max(b[2] for b in bboxes)
    north = max(b[3] for b in bboxes)

    times = [it.datetime for it in items if it.datetime]
    start = min(times)
    end = max(times)

    extent = pystac.Extent(
        spatial=pystac.SpatialExtent([[west, south, east, north]]),
        temporal=pystac.TemporalExtent([[start, end]]),
    )

    return pystac.Collection(
        id=collection_id,
        description=f"{title}：{len(items)} 景，{start:%Y-%m-%d} 至 {end:%Y-%m-%d}",
        extent=extent,
        title=title,
        license="CC-BY-4.0",
        providers=[
            pystac.Provider(
                name="示例数据提供方",
                roles=[pystac.ProviderRole.PRODUCER, pystac.ProviderRole.HOST],
                url="https://example.com",
            )
        ],
        keywords=["landsat", "surface-reflectance"],
        summaries={
            "platform": ["landsat-8", "landsat-9"],
            "gsd": pystac.Summaries({"minimum": 30, "maximum": 30}),
        },
    )

import pystac
from pystac.utils import datetime_to_str

def make_item(item_id: str, cog_path: str, href: str, acquisition: dt.datetime) -> pystac.Item:
    meta = cog_to_props(cog_path)
    bbox = meta["bbox"]
    geometry = meta["geometry"]

    item = pystac.Item(
        id=item_id,
        geometry=geometry,
        bbox=bbox,
        datetime=acquisition,
        properties={
            "platform": "landsat-8",
            "instruments": ["oli", "tirs"],
            "gsd": 30,
            "created": datetime_to_str(dt.datetime.now(dt.timezone.utc)),
            **meta["properties"],
        },
    )

    item.add_asset(
        "visual",
        pystac.Asset(
            href=href,
            media_type="image/tiff; application=geotiff; profile=cloud-optimized",
            roles=["data", "visual"],
            title="真彩色合成（B04/B03/B02）",
            extra_fields={
                "bands": [
                    {"name": "B04", "eo:common_name": "red",
                     "data_type": "uint16", "nodata": 0,
                     "raster:scale": 0.0000275, "raster:offset": -0.2},
                    {"name": "B03", "eo:common_name": "green",
                     "data_type": "uint16", "nodata": 0,
                     "raster:scale": 0.0000275, "raster:offset": -0.2},
                    {"name": "B02", "eo:common_name": "blue",
                     "data_type": "uint16", "nodata": 0,
                     "raster:scale": 0.0000275, "raster:offset": -0.2},
                ]
            },
        ),
    )

    item.add_asset(
        "thumbnail",
        pystac.Asset(
            href=href.replace(".tif", "_thumb.png"),
            media_type="image/png",
            roles=["thumbnail"],
            title="缩略图",
        ),
    )

    return item

# 分目录存放，链接仍在 Collection 上
item.set_self_href(f"{base}/items/{date:%Y/%m}/{item_id}.json")

cat = pystac.Catalog(id="demo", description="示例影像目录")
cat.add_child(collection)
collection.add_items(items)

# 1. 算 href：按 id 生成路径
cat.normalize_hrefs("https://cdn.example.com/stac")

# 2. 落盘：写本地目录
cat.save(catalog_type=pystac.CatalogType.ABSOLUTE_PUBLISHED, dest_href="./dist")

from pystac.layout import AsIsLayoutStrategy

cat.normalize_hrefs("https://cdn.example.com/stac", strategy=AsIsLayoutStrategy())

#!/usr/bin/env python3
"""从 COG 目录生成静态 STAC 目录。"""
import argparse
import datetime as dt
from pathlib import Path

import pystac

from stac_meta import cog_to_props, make_collection, make_item   # 前面定义的函数


def build(src_dir: Path, out_dir: Path, base_url: str, collection_id: str) -> None:
    cogs = sorted(src_dir.rglob("*.tif"))
    if not cogs:
        raise SystemExit(f"no .tif under {src_dir}")

    items = []
    for cog in cogs:
        # 用文件名当 id，生产上应该用数据提供方的原始 id
        item_id = cog.stem
        # 从文件名解析时间，这里假设形如 LC08_123032_20230501.tif
        date_part = item_id.split("_")[-1]
        acq = dt.datetime.strptime(date_part, "%Y%m%d").replace(tzinfo=dt.timezone.utc)

        href = f"{base_url}/imagery/{cog.parent.name}/{cog.name}"
        items.append(make_item(item_id, str(cog), href, acq))

    collection = make_collection(collection_id, "示例 Landsat 集合", items)
    collection.add_items(items)

    root = pystac.Catalog(
        id="demo-catalog",
        description="示例影像目录",
        title="示例影像目录",
    )
    root.add_child(collection)

    root.normalize_hrefs(f"{base_url}/stac")
    root.save(catalog_type=pystac.CatalogType.ABSOLUTE_PUBLISHED, dest_href=str(out_dir))

    n = root.validate_all()
    print(f"wrote {len(items)} items, validated {n}")


if __name__ == "__main__":
    ap = argparse.ArgumentParser()
    ap.add_argument("src", type=Path)
    ap.add_argument("out", type=Path)
    ap.add_argument("--base-url", required=True)
    ap.add_argument("--collection", default="demo-landsat")
    args = ap.parse_args()
    build(args.src, args.out, args.base_url.rstrip("/"), args.collection)