Description

Step purpose

Purpose: to crop a particular band defined as a common band name (such as the "green" or "nir" band) from a satellite image acquired by either Sentinel-2 or Landsat-9.

This step is highlighted below:

graph TB style B stroke:#f66,stroke-width:3px style C stroke:#f66,stroke-width:3px subgraph Process STAC item A[STAC Item] == STAC Item URL ==> B A[STAC Item] == STAC Item URL ==> C A[STAC Item] -.-> F subgraph scatter on bands B["crop(green)"]; C["crop(nir)"]; end B["crop(green)"] == crop_green.tif ==> D[Normalized difference]; C["crop(nir)"] == crop_green.tif ==> D[Normalized difference]; D -.-> E[Otsu threshold] end E -.-> F[Create STAC Catalog] F -.-> G[(storage)]

Code

The crop.py script is a command-line tool that takes as input

a SpatioTemporal Asset Catalog (STAC) Item
a bounding box area of interest (AOI), an EPSG code
a common band name as input

and then crops the specified band from the asset associated with the common band name to the specified AOI.

It uses various Python libraries like pystac, rasterio, pyproj, shapely, and loguru.

Here is an overview of the script's functionality:

It defines a function aoi2box to convert an AOI expressed as a bounding box string into a list of floats.
It defines a function get_asset to retrieve the asset of a STAC Item that is defined with a common band name. It iterates through the assets and checks if a band has the specified common name.
It defines a command-line interface using click, with options for providing the input STAC Item URL, AOI, EPSG code, and common band name.
The crop function is the main entry point. It reads the STAC Item specified by the input URL and retrieves the asset associated with the common band name. It then crops the asset to the specified AOI using the rasterio library.
It transforms the bounding box coordinates to match the EPSG code provided.
It performs the cropping using the rasterio.mask.mask function.
It writes the cropped image to a GeoTIFF file with a filename like "crop_bandname.tif."

The script is executable as a command-line tool as its usage is:

Usage: app.py [OPTIONS]

  Crops a STAC Item asset defined with its common band name

Options:
  --input-item TEXT  STAC Item URL or staged STAC catalog  [required]
  --aoi TEXT         Area of interest expressed as a bounding box  [required]
  --epsg TEXT        EPSG code  [required]
  --band TEXT        Common band name  [required]
  --help             Show this message and exit.

To use this script, you would typically run it from the command line, providing the necessary input options such as the STAC Item URL, AOI, EPSG code, and common band name:

python app.py \
  --input-item "https://earth-search.aws.element84.com/v0/collections/sentinel-s2-l2a-cogs/items/S2B_10TFK_20210713_0_L2A" \
  --aoi "-121.399,39.834,-120.74,40.472" \
  --epsg "EPSG:4326" \
  --band "green"

It will then crop the specified band from the STAC asset and save it as a GeoTIFF file.

The Python code is provided here:

water-bodies/command-line-tools/crop/app.py
import os
import click
import pystac
import rasterio
from rasterio.mask import mask
from pyproj import Transformer
from shapely import box
from loguru import logger


def aoi2box(aoi):
    """Converts an area of interest expressed as a bounding box to a list of floats"""
    return [float(c) for c in aoi.split(",")]


def get_asset(item, common_name):
    """Returns the asset of a STAC Item defined with its common band name"""
    for _, asset in item.get_assets().items():
        if not "data" in asset.to_dict()["roles"]:
            continue

        eo_asset = pystac.extensions.eo.AssetEOExtension(asset)
        if not eo_asset.bands:
            continue
        for b in eo_asset.bands:
            if (
                "common_name" in b.properties.keys()
                and b.properties["common_name"] == common_name
            ):
                return asset


@click.command(
    short_help="Crop",
    help="Crops a STAC Item asset defined with its common band name",
)
@click.option(
    "--input-item",
    "item_url",
    help="STAC Item URL or staged STAC catalog",
    required=True,
)
@click.option(
    "--aoi",
    "aoi",
    help="Area of interest expressed as a bounding box",
    required=True,
)
@click.option(
    "--epsg",
    "epsg",
    help="EPSG code",
    required=True,
)
@click.option(
    "--band",
    "band",
    help="Common band name",
    required=True,
)
def crop(item_url, aoi, band, epsg):

    if os.path.isdir(item_url):
        catalog = pystac.read_file(os.path.join(item_url, "catalog.json"))
        item = next(catalog.get_items())
    else:
        item = pystac.read_file(item_url)

    logger.info(f"Read {item.id} from {item.get_self_href()}")

    asset = get_asset(item, band)
    logger.info(f"Read asset {band} from {asset.get_absolute_href()}")

    if not asset:
        msg = f"Common band name {band} not found in the assets"
        logger.error(msg)
        raise ValueError(msg)

    bbox = aoi2box(aoi)

    with rasterio.open(asset.get_absolute_href()) as src:

        transformer = Transformer.from_crs(epsg, src.crs, always_xy=True)

        minx, miny = transformer.transform(bbox[0], bbox[1])
        maxx, maxy = transformer.transform(bbox[2], bbox[3])

        transformed_bbox = box(minx, miny, maxx, maxy)

        logger.info(f"Crop {asset.get_absolute_href()}")

        out_image, out_transform = rasterio.mask.mask(
            src, [transformed_bbox], crop=True
        )
        out_meta = src.meta.copy()

        out_meta.update(
            {
                "height": out_image.shape[1],
                "width": out_image.shape[2],
                "transform": out_transform,
                "dtype": "uint16",
                "driver": "COG",
                "tiled": True,
                "compress": "lzw",
                "blockxsize": 256,
                "blockysize": 256,
            }
        )

        with rasterio.open(f"crop_{band}.tif", "w", **out_meta) as dst_dataset:
            logger.info(f"Write crop_{band}.tif")
            dst_dataset.write(out_image)

    logger.info("Done!")


if __name__ == "__main__":
    crop()