Skip to main content
`Ctrl`+`K`

# Creating a STAC of Landsat data

In this tutorial we create a STAC of Landsat data provided by Microsoft’s Planetary Computer. There’s a lot of Landsat scenes, so we’ll only take a subset of scenes that are from a specific year and over a specific location. We’ll translate existing metadata about each scene to STAC information, utilizing the `eo`, `view`, `proj`, `raster` and `classification` extensions. Finally we’ll write out the STAC catalog to our local machine, allowing us to use stac-browser to preview the images.

## Requirements

To run this tutorial you’ll need to have installed PySTAC with the validation extra and the Planetary Computer package. To do this, use:

    pip install 'pystac[validation]' planetary-computer

Also to run this notebook you’ll need jupyter installed locally as well. If you’re running in a docker container, make sure that port `5555` is exposed if you want to run the server at the end of the notebook.

    [1]:

    import json
    from datetime import datetime
    from functools import partial
    from os.path import dirname, join
    from typing import Any, Callable, Dict, List, Optional, Tuple, Union, cast
    from urllib.parse import urlparse, urlunparse

    import planetary_computer as pc
    from dateutil.parser import parse
    from typing_extensions import TypedDict

    import pystac
    from pystac.extensions.classification import (
        Bitfield,
        Classification,
        ClassificationExtension,
    )
    from pystac.extensions.eo import Band as EOBand
    from pystac.extensions.eo import EOExtension
    from pystac.extensions.projection import ProjectionExtension
    from pystac.extensions.raster import RasterBand, RasterExtension
    from pystac.extensions.view import ViewExtension

### Identify target scenes

The Planetary Computer provides a STAC API that we could use to search for data within an area and time of interest, but since this notebook is intended to be a tutorial on creating STAC in the first place, doing so would put the cart ahead of the horse. Instead, we supply a list of metadata files for Landsat-8 and Landsat-9 scenes covering the center of Philadelphia, Pennsylvania in autumn of 2022 that we will work with:

    [2]:

    location_name = "Philly"
    scene_mtls = [
        "https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC08_L2SP_014032_20221219_20230113_02_T1/LC08_L2SP_014032_20221219_20230113_02_T1_MTL.xml",
        "https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC09_L2SP_014032_20221211_20221213_02_T2/LC09_L2SP_014032_20221211_20221213_02_T2_MTL.xml",
        "https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC08_L2SP_014032_20221203_20221212_02_T2/LC08_L2SP_014032_20221203_20221212_02_T2_MTL.xml",
        "https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC09_L2SP_014032_20221125_20230320_02_T2/LC09_L2SP_014032_20221125_20230320_02_T2_MTL.xml",
        "https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC08_L2SP_014032_20221117_20221128_02_T1/LC08_L2SP_014032_20221117_20221128_02_T1_MTL.xml",
        "https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC09_L2SP_014032_20221109_20221111_02_T1/LC09_L2SP_014032_20221109_20221111_02_T1_MTL.xml",
        "https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC08_L2SP_014032_20221101_20221114_02_T1/LC08_L2SP_014032_20221101_20221114_02_T1_MTL.xml",
        "https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC09_L2SP_014032_20221024_20221026_02_T2/LC09_L2SP_014032_20221024_20221026_02_T2_MTL.xml",
        "https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC09_L2SP_014032_20221008_20221010_02_T1/LC09_L2SP_014032_20221008_20221010_02_T1_MTL.xml",
    ]

## Read metadata from the MTL file

Landsat metadata is contained in an `MTL` file that comes in either `.txt` or `.xml` formats. We’ll rely on the XML version since it is more consistently available. This will require that we provide some facility for parsing the XML into a more usable format:

    [3]:

    # Taken from https://stackoverflow.com/questions/2148119/how-to-convert-an-xml-string-to-a-dictionary
    from xml.etree import cElementTree as ElementTree

    class XmlListConfig(list):
        def __init__(self, aList):
            for element in aList:
                if element:
                    if len(element) == 1 or element[0].tag != element[1].tag:
                        self.append(XmlDictConfig(element))
                    elif element[0].tag == element[1].tag:
                        self.append(XmlListConfig(element))
                elif element.text:
                    text = element.text.strip()
                    if text:
                        self.append(text)

    class XmlDictConfig(dict):
        def __init__(self, parent_element):
            if parent_element.items():
                self.update(dict(parent_element.items()))
            for element in parent_element:
                if element:
                    if len(element) == 1 or element[0].tag != element[1].tag:
                        aDict = XmlDictConfig(element)
                    else:
                        aDict = {element[0].tag: XmlListConfig(element)}
                    if element.items():
                        aDict.update(dict(element.items()))
                    self.update({element.tag: aDict})
                elif element.items():
                    self.update({element.tag: dict(element.items())})
                else:
                    self.update({element.tag: element.text})

We can then use these classes to get the MTL file for our scene. Notice we use `pystac.STAC_IO.read_text`; this is the method that PySTAC uses to read text as it crawls a STAC. It can read from the local filesystem or HTTP/HTTPS by default. Also, it can be extended to read from other sources such as cloud providers—see the documentation here. For now we’ll use it directly as an easy way to read a text file from an HTTPS source.

    [4]:

    stac_io = pystac.StacIO.default()

Since we’re reading our files from the Planetary Computer’s blob storage, we’re also going to have to take the additional step of signing our requests using the `planetary-computer` package’s `sign()` function. The raw URL is passed in, and the result has a shared access token applied. See the planetary-computer Python package for more details. We’ll see the use of `pc.sign()` throughout the code below, and it will be necessary for asset HREFs to be passed through this function by the user of the catalog as well.

    [5]:

    def get_metadata(xml_url: str) -> Dict[str, Any]:
        result = XmlDictConfig(ElementTree.XML(stac_io.read_text(pc.sign(xml_url))))
        result["ORIGINAL_URL"] = (
            xml_url  # Include the original URL in the metadata for use later
        )
        return result

Let’s read the MTL file for the first scene and see what it looks like.

    [6]:

    metadata = get_metadata(scene_mtls[0])
    print(json.dumps(metadata, indent=4))

    {
        "PRODUCT_CONTENTS": {
            "ORIGIN": "Image courtesy of the U.S. Geological Survey",
            "DIGITAL_OBJECT_IDENTIFIER": "https://doi.org/10.5066/P9OGBGM6",
            "LANDSAT_PRODUCT_ID": "LC08_L2SP_014032_20221219_20230113_02_T1",
            "PROCESSING_LEVEL": "L2SP",
            "COLLECTION_NUMBER": "02",
            "COLLECTION_CATEGORY": "T1",
            "OUTPUT_FORMAT": "GEOTIFF",
            "FILE_NAME_BAND_1": "LC08_L2SP_014032_20221219_20230113_02_T1_SR_B1.TIF",
            "FILE_NAME_BAND_2": "LC08_L2SP_014032_20221219_20230113_02_T1_SR_B2.TIF",
            "FILE_NAME_BAND_3": "LC08_L2SP_014032_20221219_20230113_02_T1_SR_B3.TIF",
            "FILE_NAME_BAND_4": "LC08_L2SP_014032_20221219_20230113_02_T1_SR_B4.TIF",
            "FILE_NAME_BAND_5": "LC08_L2SP_014032_20221219_20230113_02_T1_SR_B5.TIF",
            "FILE_NAME_BAND_6": "LC08_L2SP_014032_20221219_20230113_02_T1_SR_B6.TIF",
            "FILE_NAME_BAND_7": "LC08_L2SP_014032_20221219_20230113_02_T1_SR_B7.TIF",
            "FILE_NAME_BAND_ST_B10": "LC08_L2SP_014032_20221219_20230113_02_T1_ST_B10.TIF",
            "FILE_NAME_THERMAL_RADIANCE": "LC08_L2SP_014032_20221219_20230113_02_T1_ST_TRAD.TIF",
            "FILE_NAME_UPWELL_RADIANCE": "LC08_L2SP_014032_20221219_20230113_02_T1_ST_URAD.TIF",
            "FILE_NAME_DOWNWELL_RADIANCE": "LC08_L2SP_014032_20221219_20230113_02_T1_ST_DRAD.TIF",
            "FILE_NAME_ATMOSPHERIC_TRANSMITTANCE": "LC08_L2SP_014032_20221219_20230113_02_T1_ST_ATRAN.TIF",
            "FILE_NAME_EMISSIVITY": "LC08_L2SP_014032_20221219_20230113_02_T1_ST_EMIS.TIF",
            "FILE_NAME_EMISSIVITY_STDEV": "LC08_L2SP_014032_20221219_20230113_02_T1_ST_EMSD.TIF",
            "FILE_NAME_CLOUD_DISTANCE": "LC08_L2SP_014032_20221219_20230113_02_T1_ST_CDIST.TIF",
            "FILE_NAME_QUALITY_L2_AEROSOL": "LC08_L2SP_014032_20221219_20230113_02_T1_SR_QA_AEROSOL.TIF",
            "FILE_NAME_QUALITY_L2_SURFACE_TEMPERATURE": "LC08_L2SP_014032_20221219_20230113_02_T1_ST_QA.TIF",
            "FILE_NAME_QUALITY_L1_PIXEL": "LC08_L2SP_014032_20221219_20230113_02_T1_QA_PIXEL.TIF",
            "FILE_NAME_QUALITY_L1_RADIOMETRIC_SATURATION": "LC08_L2SP_014032_20221219_20230113_02_T1_QA_RADSAT.TIF",
            "FILE_NAME_ANGLE_COEFFICIENT": "LC08_L2SP_014032_20221219_20230113_02_T1_ANG.txt",
            "FILE_NAME_METADATA_ODL": "LC08_L2SP_014032_20221219_20230113_02_T1_MTL.txt",
            "FILE_NAME_METADATA_XML": "LC08_L2SP_014032_20221219_20230113_02_T1_MTL.xml",
            "DATA_TYPE_BAND_1": "UINT16",
            "DATA_TYPE_BAND_2": "UINT16",
            "DATA_TYPE_BAND_3": "UINT16",
            "DATA_TYPE_BAND_4": "UINT16",
            "DATA_TYPE_BAND_5": "UINT16",
            "DATA_TYPE_BAND_6": "UINT16",
            "DATA_TYPE_BAND_7": "UINT16",
            "DATA_TYPE_BAND_ST_B10": "UINT16",
            "DATA_TYPE_THERMAL_RADIANCE": "INT16",
            "DATA_TYPE_UPWELL_RADIANCE": "INT16",
            "DATA_TYPE_DOWNWELL_RADIANCE": "INT16",
            "DATA_TYPE_ATMOSPHERIC_TRANSMITTANCE": "INT16",
            "DATA_TYPE_EMISSIVITY": "INT16",
            "DATA_TYPE_EMISSIVITY_STDEV": "INT16",
            "DATA_TYPE_CLOUD_DISTANCE": "INT16",
            "DATA_TYPE_QUALITY_L2_AEROSOL": "UINT8",
            "DATA_TYPE_QUALITY_L2_SURFACE_TEMPERATURE": "INT16",
            "DATA_TYPE_QUALITY_L1_PIXEL": "UINT16",
            "DATA_TYPE_QUALITY_L1_RADIOMETRIC_SATURATION": "UINT16"
        },
        "IMAGE_ATTRIBUTES": {
            "SPACECRAFT_ID": "LANDSAT_8",
            "SENSOR_ID": "OLI_TIRS",
            "WRS_TYPE": "2",
            "WRS_PATH": "14",
            "WRS_ROW": "32",
            "NADIR_OFFNADIR": "NADIR",
            "TARGET_WRS_PATH": "14",
            "TARGET_WRS_ROW": "32",
            "DATE_ACQUIRED": "2022-12-19",
            "SCENE_CENTER_TIME": "15:40:17.7299160Z",
            "STATION_ID": "LGN",
            "CLOUD_COVER": "43.42",
            "CLOUD_COVER_LAND": "48.41",
            "IMAGE_QUALITY_OLI": "9",
            "IMAGE_QUALITY_TIRS": "9",
            "SATURATION_BAND_1": "N",
            "SATURATION_BAND_2": "Y",
            "SATURATION_BAND_3": "N",
            "SATURATION_BAND_4": "Y",
            "SATURATION_BAND_5": "Y",
            "SATURATION_BAND_6": "Y",
            "SATURATION_BAND_7": "Y",
            "SATURATION_BAND_8": "N",
            "SATURATION_BAND_9": "N",
            "ROLL_ANGLE": "-0.001",
            "SUN_AZIMUTH": "160.86021018",
            "SUN_ELEVATION": "23.81656674",
            "EARTH_SUN_DISTANCE": "0.9839500",
            "TRUNCATION_OLI": "UPPER",
            "TIRS_SSM_MODEL": "FINAL",
            "TIRS_SSM_POSITION_STATUS": "ESTIMATED"
        },
        "PROJECTION_ATTRIBUTES": {
            "MAP_PROJECTION": "UTM",
            "DATUM": "WGS84",
            "ELLIPSOID": "WGS84",
            "UTM_ZONE": "18",
            "GRID_CELL_SIZE_REFLECTIVE": "30.00",
            "GRID_CELL_SIZE_THERMAL": "30.00",
            "REFLECTIVE_LINES": "7861",
            "REFLECTIVE_SAMPLES": "7731",
            "THERMAL_LINES": "7861",
            "THERMAL_SAMPLES": "7731",
            "ORIENTATION": "NORTH_UP",
            "CORNER_UL_LAT_PRODUCT": "41.38441",
            "CORNER_UL_LON_PRODUCT": "-76.26178",
            "CORNER_UR_LAT_PRODUCT": "41.38140",
            "CORNER_UR_LON_PRODUCT": "-73.48833",
            "CORNER_LL_LAT_PRODUCT": "39.26052",
            "CORNER_LL_LON_PRODUCT": "-76.22284",
            "CORNER_LR_LAT_PRODUCT": "39.25773",
            "CORNER_LR_LON_PRODUCT": "-73.53498",
            "CORNER_UL_PROJECTION_X_PRODUCT": "394500.000",
            "CORNER_UL_PROJECTION_Y_PRODUCT": "4582200.000",
            "CORNER_UR_PROJECTION_X_PRODUCT": "626400.000",
            "CORNER_UR_PROJECTION_Y_PRODUCT": "4582200.000",
            "CORNER_LL_PROJECTION_X_PRODUCT": "394500.000",
            "CORNER_LL_PROJECTION_Y_PRODUCT": "4346400.000",
            "CORNER_LR_PROJECTION_X_PRODUCT": "626400.000",
            "CORNER_LR_PROJECTION_Y_PRODUCT": "4346400.000"
        },
        "LEVEL2_PROCESSING_RECORD": {
            "ORIGIN": "Image courtesy of the U.S. Geological Survey",
            "DIGITAL_OBJECT_IDENTIFIER": "https://doi.org/10.5066/P9OGBGM6",
            "REQUEST_ID": "1626123_00008",
            "LANDSAT_PRODUCT_ID": "LC08_L2SP_014032_20221219_20230113_02_T1",
            "PROCESSING_LEVEL": "L2SP",
            "OUTPUT_FORMAT": "GEOTIFF",
            "DATE_PRODUCT_GENERATED": "2023-01-13T02:53:40Z",
            "PROCESSING_SOFTWARE_VERSION": "LPGS_16.1.0",
            "ALGORITHM_SOURCE_SURFACE_REFLECTANCE": "LaSRC_1.5.0",
            "DATA_SOURCE_OZONE": "MODIS",
            "DATA_SOURCE_PRESSURE": "Calculated",
            "DATA_SOURCE_WATER_VAPOR": "MODIS",
            "DATA_SOURCE_AIR_TEMPERATURE": "MODIS",
            "ALGORITHM_SOURCE_SURFACE_TEMPERATURE": "st_1.3.0",
            "DATA_SOURCE_REANALYSIS": "GEOS-5 FP-IT"
        },
        "LEVEL2_SURFACE_REFLECTANCE_PARAMETERS": {
            "REFLECTANCE_MAXIMUM_BAND_1": "1.602213",
            "REFLECTANCE_MINIMUM_BAND_1": "-0.199972",
            "REFLECTANCE_MAXIMUM_BAND_2": "1.602213",
            "REFLECTANCE_MINIMUM_BAND_2": "-0.199972",
            "REFLECTANCE_MAXIMUM_BAND_3": "1.602213",
            "REFLECTANCE_MINIMUM_BAND_3": "-0.199972",
            "REFLECTANCE_MAXIMUM_BAND_4": "1.602213",
            "REFLECTANCE_MINIMUM_BAND_4": "-0.199972",
            "REFLECTANCE_MAXIMUM_BAND_5": "1.602213",
            "REFLECTANCE_MINIMUM_BAND_5": "-0.199972",
            "REFLECTANCE_MAXIMUM_BAND_6": "1.602213",
            "REFLECTANCE_MINIMUM_BAND_6": "-0.199972",
            "REFLECTANCE_MAXIMUM_BAND_7": "1.602213",
            "REFLECTANCE_MINIMUM_BAND_7": "-0.199972",
            "QUANTIZE_CAL_MAX_BAND_1": "65535",
            "QUANTIZE_CAL_MIN_BAND_1": "1",
            "QUANTIZE_CAL_MAX_BAND_2": "65535",
            "QUANTIZE_CAL_MIN_BAND_2": "1",
            "QUANTIZE_CAL_MAX_BAND_3": "65535",
            "QUANTIZE_CAL_MIN_BAND_3": "1",
            "QUANTIZE_CAL_MAX_BAND_4": "65535",
            "QUANTIZE_CAL_MIN_BAND_4": "1",
            "QUANTIZE_CAL_MAX_BAND_5": "65535",
            "QUANTIZE_CAL_MIN_BAND_5": "1",
            "QUANTIZE_CAL_MAX_BAND_6": "65535",
            "QUANTIZE_CAL_MIN_BAND_6": "1",
            "QUANTIZE_CAL_MAX_BAND_7": "65535",
            "QUANTIZE_CAL_MIN_BAND_7": "1",
            "REFLECTANCE_MULT_BAND_1": "2.75e-05",
            "REFLECTANCE_MULT_BAND_2": "2.75e-05",
            "REFLECTANCE_MULT_BAND_3": "2.75e-05",
            "REFLECTANCE_MULT_BAND_4": "2.75e-05",
            "REFLECTANCE_MULT_BAND_5": "2.75e-05",
            "REFLECTANCE_MULT_BAND_6": "2.75e-05",
            "REFLECTANCE_MULT_BAND_7": "2.75e-05",
            "REFLECTANCE_ADD_BAND_1": "-0.2",
            "REFLECTANCE_ADD_BAND_2": "-0.2",
            "REFLECTANCE_ADD_BAND_3": "-0.2",
            "REFLECTANCE_ADD_BAND_4": "-0.2",
            "REFLECTANCE_ADD_BAND_5": "-0.2",
            "REFLECTANCE_ADD_BAND_6": "-0.2",
            "REFLECTANCE_ADD_BAND_7": "-0.2"
        },
        "LEVEL2_SURFACE_TEMPERATURE_PARAMETERS": {
            "TEMPERATURE_MAXIMUM_BAND_ST_B10": "372.999941",
            "TEMPERATURE_MINIMUM_BAND_ST_B10": "149.003418",
            "QUANTIZE_CAL_MAXIMUM_BAND_ST_B10": "65535",
            "QUANTIZE_CAL_MINIMUM_BAND_ST_B10": "1",
            "TEMPERATURE_MULT_BAND_ST_B10": "0.00341802",
            "TEMPERATURE_ADD_BAND_ST_B10": "149.0"
        },
        "LEVEL1_PROCESSING_RECORD": {
            "ORIGIN": "Image courtesy of the U.S. Geological Survey",
            "DIGITAL_OBJECT_IDENTIFIER": "https://doi.org/10.5066/P975CC9B",
            "REQUEST_ID": "1626123_00008",
            "LANDSAT_SCENE_ID": "LC80140322022353LGN00",
            "LANDSAT_PRODUCT_ID": "LC08_L1TP_014032_20221219_20230113_02_T1",
            "PROCESSING_LEVEL": "L1TP",
            "COLLECTION_CATEGORY": "T1",
            "OUTPUT_FORMAT": "GEOTIFF",
            "DATE_PRODUCT_GENERATED": "2023-01-13T02:38:55Z",
            "PROCESSING_SOFTWARE_VERSION": "LPGS_16.1.0",
            "FILE_NAME_BAND_1": "LC08_L1TP_014032_20221219_20230113_02_T1_B1.TIF",
            "FILE_NAME_BAND_2": "LC08_L1TP_014032_20221219_20230113_02_T1_B2.TIF",
            "FILE_NAME_BAND_3": "LC08_L1TP_014032_20221219_20230113_02_T1_B3.TIF",
            "FILE_NAME_BAND_4": "LC08_L1TP_014032_20221219_20230113_02_T1_B4.TIF",
            "FILE_NAME_BAND_5": "LC08_L1TP_014032_20221219_20230113_02_T1_B5.TIF",
            "FILE_NAME_BAND_6": "LC08_L1TP_014032_20221219_20230113_02_T1_B6.TIF",
            "FILE_NAME_BAND_7": "LC08_L1TP_014032_20221219_20230113_02_T1_B7.TIF",
            "FILE_NAME_BAND_8": "LC08_L1TP_014032_20221219_20230113_02_T1_B8.TIF",
            "FILE_NAME_BAND_9": "LC08_L1TP_014032_20221219_20230113_02_T1_B9.TIF",
            "FILE_NAME_BAND_10": "LC08_L1TP_014032_20221219_20230113_02_T1_B10.TIF",
            "FILE_NAME_BAND_11": "LC08_L1TP_014032_20221219_20230113_02_T1_B11.TIF",
            "FILE_NAME_QUALITY_L1_PIXEL": "LC08_L1TP_014032_20221219_20230113_02_T1_QA_PIXEL.TIF",
            "FILE_NAME_QUALITY_L1_RADIOMETRIC_SATURATION": "LC08_L1TP_014032_20221219_20230113_02_T1_QA_RADSAT.TIF",
            "FILE_NAME_ANGLE_COEFFICIENT": "LC08_L1TP_014032_20221219_20230113_02_T1_ANG.txt",
            "FILE_NAME_ANGLE_SENSOR_AZIMUTH_BAND_4": "LC08_L1TP_014032_20221219_20230113_02_T1_VAA.TIF",
            "FILE_NAME_ANGLE_SENSOR_ZENITH_BAND_4": "LC08_L1TP_014032_20221219_20230113_02_T1_VZA.TIF",
            "FILE_NAME_ANGLE_SOLAR_AZIMUTH_BAND_4": "LC08_L1TP_014032_20221219_20230113_02_T1_SAA.TIF",
            "FILE_NAME_ANGLE_SOLAR_ZENITH_BAND_4": "LC08_L1TP_014032_20221219_20230113_02_T1_SZA.TIF",
            "FILE_NAME_METADATA_ODL": "LC08_L1TP_014032_20221219_20230113_02_T1_MTL.txt",
            "FILE_NAME_METADATA_XML": "LC08_L1TP_014032_20221219_20230113_02_T1_MTL.xml",
            "FILE_NAME_CPF": "LC08CPF_20221001_20221231_02.03",
            "FILE_NAME_BPF_OLI": "LO8BPF20221219152831_20221219170353.01",
            "FILE_NAME_BPF_TIRS": "LT8BPF20221215135451_20221222101440.01",
            "FILE_NAME_RLUT": "LC08RLUT_20150303_20431231_02_01.h5",
            "DATA_SOURCE_TIRS_STRAY_LIGHT_CORRECTION": "TIRS",
            "DATA_SOURCE_ELEVATION": "GLS2000",
            "GROUND_CONTROL_POINTS_VERSION": "5",
            "GROUND_CONTROL_POINTS_MODEL": "462",
            "GEOMETRIC_RMSE_MODEL": "8.179",
            "GEOMETRIC_RMSE_MODEL_Y": "7.213",
            "GEOMETRIC_RMSE_MODEL_X": "3.856",
            "GROUND_CONTROL_POINTS_VERIFY": "120",
            "GEOMETRIC_RMSE_VERIFY": "7.426"
        },
        "LEVEL1_MIN_MAX_RADIANCE": {
            "RADIANCE_MAXIMUM_BAND_1": "785.06079",
            "RADIANCE_MINIMUM_BAND_1": "-64.83057",
            "RADIANCE_MAXIMUM_BAND_2": "803.91187",
            "RADIANCE_MINIMUM_BAND_2": "-66.38730",
            "RADIANCE_MAXIMUM_BAND_3": "740.79791",
            "RADIANCE_MINIMUM_BAND_3": "-61.17533",
            "RADIANCE_MAXIMUM_BAND_4": "624.68250",
            "RADIANCE_MINIMUM_BAND_4": "-51.58648",
            "RADIANCE_MAXIMUM_BAND_5": "382.27454",
            "RADIANCE_MINIMUM_BAND_5": "-31.56836",
            "RADIANCE_MAXIMUM_BAND_6": "95.06820",
            "RADIANCE_MINIMUM_BAND_6": "-7.85076",
            "RADIANCE_MAXIMUM_BAND_7": "32.04307",
            "RADIANCE_MINIMUM_BAND_7": "-2.64613",
            "RADIANCE_MAXIMUM_BAND_8": "706.96869",
            "RADIANCE_MINIMUM_BAND_8": "-58.38170",
            "RADIANCE_MAXIMUM_BAND_9": "149.40157",
            "RADIANCE_MINIMUM_BAND_9": "-12.33763",
            "RADIANCE_MAXIMUM_BAND_10": "22.00180",
            "RADIANCE_MINIMUM_BAND_10": "0.10033",
            "RADIANCE_MAXIMUM_BAND_11": "22.00180",
            "RADIANCE_MINIMUM_BAND_11": "0.10033"
        },
        "LEVEL1_MIN_MAX_REFLECTANCE": {
            "REFLECTANCE_MAXIMUM_BAND_1": "1.210700",
            "REFLECTANCE_MINIMUM_BAND_1": "-0.099980",
            "REFLECTANCE_MAXIMUM_BAND_2": "1.210700",
            "REFLECTANCE_MINIMUM_BAND_2": "-0.099980",
            "REFLECTANCE_MAXIMUM_BAND_3": "1.210700",
            "REFLECTANCE_MINIMUM_BAND_3": "-0.099980",
            "REFLECTANCE_MAXIMUM_BAND_4": "1.210700",
            "REFLECTANCE_MINIMUM_BAND_4": "-0.099980",
            "REFLECTANCE_MAXIMUM_BAND_5": "1.210700",
            "REFLECTANCE_MINIMUM_BAND_5": "-0.099980",
            "REFLECTANCE_MAXIMUM_BAND_6": "1.210700",
            "REFLECTANCE_MINIMUM_BAND_6": "-0.099980",
            "REFLECTANCE_MAXIMUM_BAND_7": "1.210700",
            "REFLECTANCE_MINIMUM_BAND_7": "-0.099980",
            "REFLECTANCE_MAXIMUM_BAND_8": "1.210700",
            "REFLECTANCE_MINIMUM_BAND_8": "-0.099980",
            "REFLECTANCE_MAXIMUM_BAND_9": "1.210700",
            "REFLECTANCE_MINIMUM_BAND_9": "-0.099980"
        },
        "LEVEL1_MIN_MAX_PIXEL_VALUE": {
            "QUANTIZE_CAL_MAX_BAND_1": "65535",
            "QUANTIZE_CAL_MIN_BAND_1": "1",
            "QUANTIZE_CAL_MAX_BAND_2": "65535",
            "QUANTIZE_CAL_MIN_BAND_2": "1",
            "QUANTIZE_CAL_MAX_BAND_3": "65535",
            "QUANTIZE_CAL_MIN_BAND_3": "1",
            "QUANTIZE_CAL_MAX_BAND_4": "65535",
            "QUANTIZE_CAL_MIN_BAND_4": "1",
            "QUANTIZE_CAL_MAX_BAND_5": "65535",
            "QUANTIZE_CAL_MIN_BAND_5": "1",
            "QUANTIZE_CAL_MAX_BAND_6": "65535",
            "QUANTIZE_CAL_MIN_BAND_6": "1",
            "QUANTIZE_CAL_MAX_BAND_7": "65535",
            "QUANTIZE_CAL_MIN_BAND_7": "1",
            "QUANTIZE_CAL_MAX_BAND_8": "65535",
            "QUANTIZE_CAL_MIN_BAND_8": "1",
            "QUANTIZE_CAL_MAX_BAND_9": "65535",
            "QUANTIZE_CAL_MIN_BAND_9": "1",
            "QUANTIZE_CAL_MAX_BAND_10": "65535",
            "QUANTIZE_CAL_MIN_BAND_10": "1",
            "QUANTIZE_CAL_MAX_BAND_11": "65535",
            "QUANTIZE_CAL_MIN_BAND_11": "1"
        },
        "LEVEL1_RADIOMETRIC_RESCALING": {
            "RADIANCE_MULT_BAND_1": "1.2969E-02",
            "RADIANCE_MULT_BAND_2": "1.3280E-02",
            "RADIANCE_MULT_BAND_3": "1.2238E-02",
            "RADIANCE_MULT_BAND_4": "1.0319E-02",
            "RADIANCE_MULT_BAND_5": "6.3149E-03",
            "RADIANCE_MULT_BAND_6": "1.5705E-03",
            "RADIANCE_MULT_BAND_7": "5.2933E-04",
            "RADIANCE_MULT_BAND_8": "1.1679E-02",
            "RADIANCE_MULT_BAND_9": "2.4680E-03",
            "RADIANCE_MULT_BAND_10": "3.3420E-04",
            "RADIANCE_MULT_BAND_11": "3.3420E-04",
            "RADIANCE_ADD_BAND_1": "-64.84355",
            "RADIANCE_ADD_BAND_2": "-66.40058",
            "RADIANCE_ADD_BAND_3": "-61.18757",
            "RADIANCE_ADD_BAND_4": "-51.59680",
            "RADIANCE_ADD_BAND_5": "-31.57467",
            "RADIANCE_ADD_BAND_6": "-7.85233",
            "RADIANCE_ADD_BAND_7": "-2.64666",
            "RADIANCE_ADD_BAND_8": "-58.39338",
            "RADIANCE_ADD_BAND_9": "-12.34010",
            "RADIANCE_ADD_BAND_10": "0.10000",
            "RADIANCE_ADD_BAND_11": "0.10000",
            "REFLECTANCE_MULT_BAND_1": "2.0000E-05",
            "REFLECTANCE_MULT_BAND_2": "2.0000E-05",
            "REFLECTANCE_MULT_BAND_3": "2.0000E-05",
            "REFLECTANCE_MULT_BAND_4": "2.0000E-05",
            "REFLECTANCE_MULT_BAND_5": "2.0000E-05",
            "REFLECTANCE_MULT_BAND_6": "2.0000E-05",
            "REFLECTANCE_MULT_BAND_7": "2.0000E-05",
            "REFLECTANCE_MULT_BAND_8": "2.0000E-05",
            "REFLECTANCE_MULT_BAND_9": "2.0000E-05",
            "REFLECTANCE_ADD_BAND_1": "-0.100000",
            "REFLECTANCE_ADD_BAND_2": "-0.100000",
            "REFLECTANCE_ADD_BAND_3": "-0.100000",
            "REFLECTANCE_ADD_BAND_4": "-0.100000",
            "REFLECTANCE_ADD_BAND_5": "-0.100000",
            "REFLECTANCE_ADD_BAND_6": "-0.100000",
            "REFLECTANCE_ADD_BAND_7": "-0.100000",
            "REFLECTANCE_ADD_BAND_8": "-0.100000",
            "REFLECTANCE_ADD_BAND_9": "-0.100000"
        },
        "LEVEL1_THERMAL_CONSTANTS": {
            "K1_CONSTANT_BAND_10": "774.8853",
            "K2_CONSTANT_BAND_10": "1321.0789",
            "K1_CONSTANT_BAND_11": "480.8883",
            "K2_CONSTANT_BAND_11": "1201.1442"
        },
        "LEVEL1_PROJECTION_PARAMETERS": {
            "MAP_PROJECTION": "UTM",
            "DATUM": "WGS84",
            "ELLIPSOID": "WGS84",
            "UTM_ZONE": "18",
            "GRID_CELL_SIZE_PANCHROMATIC": "15.00",
            "GRID_CELL_SIZE_REFLECTIVE": "30.00",
            "GRID_CELL_SIZE_THERMAL": "30.00",
            "ORIENTATION": "NORTH_UP",
            "RESAMPLING_OPTION": "CUBIC_CONVOLUTION"
        },
        "ORIGINAL_URL": "https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC08_L2SP_014032_20221219_20230113_02_T1/LC08_L2SP_014032_20221219_20230113_02_T1_MTL.xml"
    }

There are a number of files referred to by this metadata file which are in the same tree in the cloud. We must provide an easy means for creating a URL for these sidecar files. We can then use `partial` to create a helper function to turn a file name into a URL.

    [7]:

    def download_sidecar(metadata: Dict[str, Any], filename: str) -> str:
        parsed = urlparse(metadata["ORIGINAL_URL"])
        return urlunparse(parsed._replace(path=join(dirname(parsed.path), filename)))

    [8]:

    download_url = partial(download_sidecar, metadata)

### Create a STAC Item from a scene

Now that we have metadata for the scene let’s use it to create a STAC Item.

We can use the `help` method to see the signature of the `__init__` method on `pystac.Item`. You can also call `help` directly on `pystac.Item` for broader documentation, or check the API docs for Item here.

    [9]:

    help(pystac.Item.__init__)

    Help on function __init__ in module pystac.item:

    __init__(self, id: 'str', geometry: 'Optional[Dict[str, Any]]', bbox: 'Optional[List[float]]', datetime: 'Optional[Datetime]', properties: 'Dict[str, Any]', start_datetime: 'Optional[Datetime]' = None, end_datetime: 'Optional[Datetime]' = None, stac_extensions: 'Optional[List[str]]' = None, href: 'Optional[str]' = None, collection: 'Optional[Union[str, Collection]]' = None, extra_fields: 'Optional[Dict[str, Any]]' = None, assets: 'Optional[Dict[str, Asset]]' = None)
        Initialize self.  See help(type(self)) for accurate signature.

We can see we’ll need at least an `id`, `geometry`, `bbox`, and `datetime`. Properties is required, but can be an empty dictionary that we fill out on the Item once it’s created.

> Caution! The `Optional` type hint is used when None can be provided in place of a meaningful argument; it does not indicate that the argument does not need to be supplied—that is only true if a default value is indicated in the type signature.

#### Item `id`

For the Item’s `id`, we’ll use the scene ID. We’ll chop off the last 5 characters as they are repeated for each ID and so aren’t necessary:

    [10]:

    def get_item_id(metadata: Dict[str, Any]) -> str:
        return cast(str, metadata["LEVEL1_PROCESSING_RECORD"]["LANDSAT_SCENE_ID"][:-5])

    [11]:

    item_id = get_item_id(metadata)
    item_id

    [11]:

    'LC80140322022353'

#### Item `datetime`

Here we parse the datetime of the Item from two metadata fields that describe the date and time the scene was captured:

    [12]:

    def get_datetime(metadata: Dict[str, Any]) -> datetime:
        return parse(
            "%sT%s"
            % (
                metadata["IMAGE_ATTRIBUTES"]["DATE_ACQUIRED"],
                metadata["IMAGE_ATTRIBUTES"]["SCENE_CENTER_TIME"],
            )
        )

    [13]:

    item_datetime = get_datetime(metadata)
    item_datetime

    [13]:

    datetime.datetime(2022, 12, 19, 15, 40, 17, 729916, tzinfo=tzutc())

#### Item `bbox`

Here we read in the bounding box information from the scene and transform it into the format of the Item’s `bbox` property:

    [14]:

    def get_bbox(metadata: Dict[str, Any]) -> List[float]:
        metadata = metadata["PROJECTION_ATTRIBUTES"]
        coords = [
            [
                [
                    float(metadata["CORNER_UL_LON_PRODUCT"]),
                    float(metadata["CORNER_UL_LAT_PRODUCT"]),
                ],
                [
                    float(metadata["CORNER_UR_LON_PRODUCT"]),
                    float(metadata["CORNER_UR_LAT_PRODUCT"]),
                ],
                [
                    float(metadata["CORNER_LR_LON_PRODUCT"]),
                    float(metadata["CORNER_LR_LAT_PRODUCT"]),
                ],
                [
                    float(metadata["CORNER_LL_LON_PRODUCT"]),
                    float(metadata["CORNER_LL_LAT_PRODUCT"]),
                ],
                [
                    float(metadata["CORNER_UL_LON_PRODUCT"]),
                    float(metadata["CORNER_UL_LAT_PRODUCT"]),
                ],
            ]
        ]
        lats = [c[1] for c in coords[0]]
        lons = [c[0] for c in coords[0]]
        return [min(lons), min(lats), max(lons), max(lats)]

    [15]:

    item_bbox = get_bbox(metadata)
    item_bbox

    [15]:

    [-76.26178, 39.25773, -73.48833, 41.38441]

#### Item `geometry`

Getting the geometry of the scene is a little more tricky. The bounding box will be a axis-aligned rectangle of the area the scene occupies, but will not represent the true footprint of the image - Landsat scenes are “tilted” according the the coordinate reference system, so there will be areas in the corner where no image data exists. When constructing a STAC Item it’s best if you have the Item geometry represent the true footprint of the assets.

To get the footprint of the scene we’ll read in another metadata file that lives alongside the MTL - the `ANG.txt` file. This function uses the ANG file and the bbox to construct the GeoJSON polygon that represents the footprint of the scene:

    [16]:

    def get_geometry(metadata: Dict[str, Any], bbox: List[float]) -> Dict[str, Any]:
        url = download_sidecar(
            metadata, metadata["PRODUCT_CONTENTS"]["FILE_NAME_ANGLE_COEFFICIENT"]
        )
        sz = []
        coords = []
        ang_text = stac_io.read_text(pc.sign(url))
        if not ang_text.startswith("GROUP"):
            raise ValueError(f"ANG file for url {url} is incorrectly formatted")
        for line in ang_text.split("\n"):
            if "BAND01_NUM_L1T_LINES" in line or "BAND01_NUM_L1T_SAMPS" in line:
                sz.append(float(line.split("=")[1]))
            if (
                "BAND01_L1T_IMAGE_CORNER_LINES" in line
                or "BAND01_L1T_IMAGE_CORNER_SAMPS" in line
            ):
                coords.append(
                    [float(v) for v in line.split("=")[1].strip().strip("()").split(",")]
                )
            if len(coords) == 2:
                break
        dlon = bbox[2] - bbox[0]
        dlat = bbox[3] - bbox[1]
        lons = [c / sz[1] * dlon + bbox[0] for c in coords[1]]
        lats = [((sz[0] - c) / sz[0]) * dlat + bbox[1] for c in coords[0]]
        coordinates = [
            [
                [lons[0], lats[0]],
                [lons[1], lats[1]],
                [lons[2], lats[2]],
                [lons[3], lats[3]],
                [lons[0], lats[0]],
            ]
        ]

        return {"type": "Polygon", "coordinates": coordinates}

    [17]:

    item_geometry = get_geometry(metadata, item_bbox)
    print(json.dumps(item_geometry, indent=2))

    {
      "type": "Polygon",
      "coordinates": [
        [
          [
            -75.71075270108336,
            41.3823086369878
          ],
          [
            -73.48924866988654,
            40.980654308234485
          ],
          [
            -74.0425618957281,
            39.25823722657151
          ],
          [
            -76.26093009667797,
            39.66800780107756
          ],
          [
            -75.71075270108336,
            41.3823086369878
          ]
        ]
      ]
    }

This would be a good time to check our work - we can print out the GeoJSON and use geojson.io to check and make sure we’re using scenes that overlap our location. If this footprint is somewhere unexpected in the world, make sure the Lat/Long coordinates are correct and in the right order!

#### Create the item

Now that we have the required attributes for an Item we can create it:

    [18]:

    item = pystac.Item(
        id=item_id,
        datetime=item_datetime,
        geometry=item_geometry,
        bbox=item_bbox,
        properties={},
    )

PySTAC has a `validate` method on STAC objects, which you can use to make sure you’re constructing things correctly. If there’s an issue the following line will throw an exception:

    [19]:

    item.validate()

    [19]:

    ['https://schemas.stacspec.org/v1.0.0/item-spec/json-schema/item.json']

#### Add Ground Sample Distance to common metadata

We’ll add the Ground Sample Distance that is defined as part of the Item Common Metadata. We define this on the Item level as 30 meters, which is the GSD for most of the Landsat bands. However, if some bands have a different resolution; we can account for this by setting the GSD explicitly for each of those bands below.

    [20]:

    item.common_metadata.gsd = 30.0

#### Adding the EO extension

STAC has a rich set of extensions that allow STAC objects to encode information that is not part of the core spec but is used widely and standardized. These extensions allow us to augment STAC objects with additional structured metadata that describe referenced data with semantically-meaningful fields. An example of this is the eo extension, which captures fields needed for electro-optical data, like center wavelength and full-width half maximum values.

This notebook will also rely on other extensions; but as they will apply to different objects, not just the item itself, they will be invoked later.

For now, we will enable the EO extension for this item by using the `ext` property provided by the extension object:

    [21]:

    eo_ext = EOExtension.ext(item, add_if_missing=True)

#### Add cloud cover

Here we add cloud cover from the metadata as part of the `eo` extension.

    [22]:

    def get_cloud_cover(metadata: Dict[str, Any]) -> float:
        return float(metadata["IMAGE_ATTRIBUTES"]["CLOUD_COVER"])

    [23]:

    eo_ext.cloud_cover = get_cloud_cover(metadata)
    eo_ext.cloud_cover

    [23]:

    43.42

## Adding assets

STAC Items contain a list of Assets, which are a list of files that relate to the Item. In our case we’ll be cataloging each file related to the scene, including the Landsat band files as well as the metadata files associated with the scene.

Each asset will have a name, some basic properties, and then possibly some properties defined by the various extensions in use (`eo`, `raster`, and `classification`). So, we begin by defining a type alias for this package of information and some helper functions for creating them:

    [24]:

    class BandInfo(TypedDict):
        name: str
        asset_fields: Dict[str, str]
        extensions: List[Union[EOBand, RasterBand, List[Bitfield]]]

    def eo_band_info(
        common_name: str,
        href: str,
        name: str,
        description: str,
        center: float,
        fwhm: float,
        default_raster_band: Optional[RasterBand] = None,
    ):
        raster_band = (
            RasterBand.create(
                spatial_resolution=30.0,
                scale=0.0000275,
                nodata=0,
                offset=-0.2,
                data_type="uint16",
            )
            if default_raster_band is None
            else default_raster_band
        )
        return {
            "name": common_name,
            "asset_fields": {
                "href": href,
                "media_type": str(pystac.media_type.MediaType.COG),
            },
            "extensions": [
                EOBand.create(
                    name=name,
                    common_name=common_name,
                    description=description,
                    center_wavelength=center,
                    full_width_half_max=fwhm,
                ),
                raster_band,
            ],
        }

    def plain_band_info(name: str, href: str, title: str, ext: RasterBand):
        return {
            "name": name,
            "asset_fields": {
                "href": href,
                "media_type": str(pystac.media_type.MediaType.COG),
                "title": title,
            },
            "extensions": [ext],
        }

Some common raster band information definitions will also be useful.

    [25]:

    thermal_raster_band = RasterBand.create(
        spatial_resolution=30.0,
        scale=0.00341802,
        nodata=0,
        offset=149.0,
        data_type="uint6",
        unit="kelvin",
    )
    radiance_raster_band = RasterBand.create(
        unit="watt/steradian/square_meter/micrometer",
        scale=1e-3,
        nodata=-9999,
        data_type="uint16",
        spatial_resolution=30.0,
    )
    emissivity_transmission_raster_band = RasterBand.create(
        scale=1e-4,
        nodata=-9999,
        data_type="int16",
        spatial_resolution=30.0,
    )

Several QA bands are provided that utilize bit-wise masks which we can define using the classification extension. Because these definitions can be verbose, we provide some additional helper functions to minimize the length of their definition.

    [26]:

    def create_bitfield(
        offset: int,
        length: int,
        name: str,
        field_names_descriptions: List[Tuple[str, str]],
        description: Optional[str] = None,
    ) -> Bitfield:
        return Bitfield.create(
            offset=offset,
            length=length,
            name=name,
            description=description,
            classes=[
                Classification.create(value=i, name=n, description=d)
                for (i, (n, d)) in enumerate(field_names_descriptions)
            ],
        )

    def create_qa_bitfield(
        offset: int,
        class_name: str,
        description: Optional[Union[str, Tuple[str, str]]] = None,
    ) -> Bitfield:
        if description is None:
            descr0 = f"{class_name.replace('_', ' ').capitalize()} confidence is not high"
            descr1 = f"High confidence {class_name.replace('_', ' ')}"
        elif isinstance(description, str):
            descr0 = f"{description.capitalize()} confidence is not high"
            descr1 = f"High confidence {description.lower()}"
        else:
            descr0 = description[0]
            descr1 = description[1]

        return create_bitfield(
            offset, 1, class_name, [(f"not_{class_name}", descr0), (class_name, descr1)]
        )

    def create_confidence_bitfield(
        offset: int, class_name: str, use_medium: bool = False
    ) -> Bitfield:
        label = class_name.replace("_", " ")
        return create_bitfield(
            offset,
            2,
            f"{class_name}_confidence",
            [
                ("not_set", "No confidence level set"),
                ("low", f"Low confidence {label}"),
                (
                    ("medium", f"Medium confidence {label}")
                    if use_medium
                    else ("reserved", "Reserved - value not used")
                ),
                ("high", f"High confidence {label}"),
            ],
        )

We can now create the `BandInfo` definitions for the Landsat scenes. This begins with the definition of a function to convert metadata into a list of `BandInfo` records, which is lengthy but ultimately straightforward.

    [27]:

    def landsat_band_info(
        metadata: Dict[str, Any], downloader: Callable[[str], str]
    ) -> List[BandInfo]:
        product_contents = metadata["PRODUCT_CONTENTS"]
        return [
            {
                "name": "ang",
                "asset_fields": {
                    "href": downloader(product_contents["FILE_NAME_ANGLE_COEFFICIENT"]),
                    "media_type": "text/plain",
                    "title": "Angle coefficients",
                },
                "extensions": [],
            },
            {
                "name": "mtl.txt",
                "asset_fields": {
                    "href": downloader(product_contents["FILE_NAME_METADATA_ODL"]),
                    "media_type": "text/plain",
                    "title": "Product metadata",
                },
                "extensions": [],
            },
            eo_band_info(
                "coastal",
                downloader(product_contents["FILE_NAME_BAND_1"]),
                "OLI_B1",
                "Coastal/Aerosol (Operational Land Imager)",
                0.44,
                0.02,
            ),
            eo_band_info(
                "blue",
                downloader(product_contents["FILE_NAME_BAND_2"]),
                "OLI_B2",
                "Visible blue (Operational Land Imager)",
                0.48,
                0.06,
            ),
            eo_band_info(
                "green",
                downloader(product_contents["FILE_NAME_BAND_3"]),
                "OLI_B3",
                "Visible green (Operational Land Imager)",
                0.56,
                0.06,
            ),
            eo_band_info(
                "red",
                downloader(product_contents["FILE_NAME_BAND_4"]),
                "OLI_B4",
                "Visible red (Operational Land Imager)",
                0.65,
                0.04,
            ),
            eo_band_info(
                "nir08",
                downloader(product_contents["FILE_NAME_BAND_5"]),
                "OLI_B5",
                "Near infrared (Operational Land Imager)",
                0.87,
                0.03,
            ),
            eo_band_info(
                "swir16",
                downloader(product_contents["FILE_NAME_BAND_6"]),
                "OLI_B6",
                "Short-wave infrared (Operational Land Imager)",
                1.61,
                0.09,
            ),
            eo_band_info(
                "swir22",
                downloader(product_contents["FILE_NAME_BAND_7"]),
                "OLI_B7",
                "Short-wave infrared (Operational Land Imager)",
                2.2,
                0.19,
            ),
            eo_band_info(
                "lwir11",
                downloader(product_contents["FILE_NAME_BAND_ST_B10"]),
                "TIRS_B10",
                "Long-wave infrared (Thermal InfraRed Sensor)",
                10.9,
                0.59,
                thermal_raster_band,
            ),
            plain_band_info(
                "trad",
                downloader(product_contents["FILE_NAME_THERMAL_RADIANCE"]),
                "Thermal radiance",
                radiance_raster_band,
            ),
            plain_band_info(
                "urad",
                downloader(product_contents["FILE_NAME_UPWELL_RADIANCE"]),
                "Upwelled radiance",
                radiance_raster_band,
            ),
            plain_band_info(
                "drad",
                downloader(product_contents["FILE_NAME_DOWNWELL_RADIANCE"]),
                "Downwelled radiance",
                radiance_raster_band,
            ),
            plain_band_info(
                "emis",
                downloader(product_contents["FILE_NAME_EMISSIVITY"]),
                "Emissivity",
                emissivity_transmission_raster_band,
            ),
            plain_band_info(
                "emsd",
                downloader(product_contents["FILE_NAME_EMISSIVITY_STDEV"]),
                "Emissivity standard deviation",
                emissivity_transmission_raster_band,
            ),
            plain_band_info(
                "atran",
                downloader(product_contents["FILE_NAME_ATMOSPHERIC_TRANSMITTANCE"]),
                "Atmospheric transmission",
                emissivity_transmission_raster_band,
            ),
            plain_band_info(
                "cdist",
                downloader(product_contents["FILE_NAME_CLOUD_DISTANCE"]),
                "Cloud distance",
                RasterBand.create(
                    unit="kilometer",
                    scale=1e-2,
                    nodata=-9999,
                    data_type="uint16",
                    spatial_resolution=30.0,
                ),
            ),
            {
                "name": "qa",
                "asset_fields": {
                    "href": downloader(
                        product_contents["FILE_NAME_QUALITY_L2_SURFACE_TEMPERATURE"]
                    ),
                    "title": "Surface Temperature Quality Assessment Band",
                },
                "extensions": [
                    RasterBand.create(
                        unit="kelvin",
                        scale=1e-2,
                        nodata=-9999,
                        data_type="int16",
                        spatial_resolution=30,
                    )
                ],
            },
            {
                "name": "qa_pixel",
                "asset_fields": {
                    "href": downloader(product_contents["FILE_NAME_QUALITY_L1_PIXEL"]),
                    "media_type": str(pystac.media_type.MediaType.COG),
                    "title": "Pixel quality assessment",
                },
                "extensions": [
                    [
                        create_qa_bitfield(0, "fill", ("Image data", "Fill data")),
                        create_qa_bitfield(
                            1,
                            "dilated_cloud",
                            ("Cloud is not dilated or no cloud", "Dilated cloud"),
                        ),
                        create_qa_bitfield(2, "cirrus"),
                        create_qa_bitfield(3, "cloud"),
                        create_qa_bitfield(4, "cloud_shadow"),
                        create_qa_bitfield(5, "snow"),
                        create_qa_bitfield(6, "clear"),
                        create_qa_bitfield(7, "water"),
                        create_confidence_bitfield(8, "cloud", True),
                        create_confidence_bitfield(10, "cloud_shadow"),
                        create_confidence_bitfield(12, "snow"),
                        create_confidence_bitfield(14, "cirrus"),
                    ]
                ],
            },
            {
                "name": "qa_radsat",
                "asset_fields": {
                    "href": downloader(
                        product_contents["FILE_NAME_QUALITY_L1_RADIOMETRIC_SATURATION"]
                    ),
                    "media_type": str(pystac.media_type.MediaType.COG),
                    "description": (
                        "Collection 2 Level-1 Radiometric Saturation and "
                        "Terrain Occlusion Quality Assessment "
                        "Band (QA_RADSAT)"
                    ),
                },
                "extensions": [
                    [
                        Bitfield.create(
                            offset=i - 1,
                            length=1,
                            description=f"Band {i} radiometric saturation",
                            classes=[
                                Classification.create(
                                    0, f"Band {i} not saturated", "not_saturated"
                                ),
                                Classification.create(
                                    1, f"Band {i} saturated", "saturated"
                                ),
                            ],
                        )
                        for i in [1, 2, 3, 4, 5, 6, 7, 9]
                    ]
                    + [
                        Bitfield.create(
                            offset=11,
                            length=1,
                            description=(
                                "Terrain not visible from sensor due to intervening terrain"
                            ),
                            classes=[
                                Classification.create(
                                    0, "Terrain is not occluded", "not_occluded"
                                ),
                                Classification.create(1, "Terrain is occluded", "occluded"),
                            ],
                        )
                    ]
                ],
            },
            {
                "name": "qa_aerosol",
                "asset_fields": {
                    "href": downloader(product_contents["FILE_NAME_QUALITY_L2_AEROSOL"]),
                    "media_type": str(pystac.media_type.MediaType.COG),
                    "title": "Aerosol Quality Assessment Band",
                },
                "extensions": [
                    [
                        create_bitfield(
                            0,
                            1,
                            "fill",
                            [("not_fill", "Pixel is not fill"), ("fill", "Pixel is fill")],
                            "Image or fill data",
                        ),
                        create_bitfield(
                            1,
                            1,
                            "retrieval",
                            [
                                ("not_valid", "Pixel retrieval is not valid"),
                                ("valid", "Pixel retrieval is valid"),
                            ],
                            "Valid aerosol retrieval",
                        ),
                        create_bitfield(
                            2,
                            1,
                            "water",
                            [
                                ("not_water", "Pixel is not water"),
                                ("water", "Pixel is water"),
                            ],
                            "Water mask",
                        ),
                        create_bitfield(
                            5,
                            1,
                            "interpolated",
                            [
                                ("not_interpolated", "Pixel is not interpolated"),
                                ("interpolated", "Pixel is interpolated"),
                            ],
                            "Aerosol interpolation",
                        ),
                        create_bitfield(
                            6,
                            2,
                            "level",
                            [
                                ("climatology", "No aerosol correction applied"),
                                ("low", "Low aerosol level"),
                                ("medium", "Medium aerosol level"),
                                ("high", "High aerosol level"),
                            ],
                            "Aerosol level",
                        ),
                    ]
                ],
            },
        ]

For illustration purposes, we can look at the band info records for an example scene:

    [28]:

    landsat_band_info(metadata, download_url)

    [28]:

    [{'name': 'ang',
      'asset_fields': {'href': 'https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC08_L2SP_014032_20221219_20230113_02_T1/LC08_L2SP_014032_20221219_20230113_02_T1_ANG.txt',
       'media_type': 'text/plain',
       'title': 'Angle coefficients'},
      'extensions': []},
     {'name': 'mtl.txt',
      'asset_fields': {'href': 'https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC08_L2SP_014032_20221219_20230113_02_T1/LC08_L2SP_014032_20221219_20230113_02_T1_MTL.txt',
       'media_type': 'text/plain',
       'title': 'Product metadata'},
      'extensions': []},
     {'name': 'coastal',
      'asset_fields': {'href': 'https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC08_L2SP_014032_20221219_20230113_02_T1/LC08_L2SP_014032_20221219_20230113_02_T1_SR_B1.TIF',
       'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized'},
      'extensions': [<Band name=OLI_B1>]},
     {'name': 'blue',
      'asset_fields': {'href': 'https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC08_L2SP_014032_20221219_20230113_02_T1/LC08_L2SP_014032_20221219_20230113_02_T1_SR_B2.TIF',
       'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized'},
      'extensions': [<Band name=OLI_B2>]},
     {'name': 'green',
      'asset_fields': {'href': 'https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC08_L2SP_014032_20221219_20230113_02_T1/LC08_L2SP_014032_20221219_20230113_02_T1_SR_B3.TIF',
       'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized'},
      'extensions': [<Band name=OLI_B3>]},
     {'name': 'red',
      'asset_fields': {'href': 'https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC08_L2SP_014032_20221219_20230113_02_T1/LC08_L2SP_014032_20221219_20230113_02_T1_SR_B4.TIF',
       'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized'},
      'extensions': [<Band name=OLI_B4>]},
     {'name': 'nir08',
      'asset_fields': {'href': 'https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC08_L2SP_014032_20221219_20230113_02_T1/LC08_L2SP_014032_20221219_20230113_02_T1_SR_B5.TIF',
       'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized'},
      'extensions': [<Band name=OLI_B5>]},
     {'name': 'swir16',
      'asset_fields': {'href': 'https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC08_L2SP_014032_20221219_20230113_02_T1/LC08_L2SP_014032_20221219_20230113_02_T1_SR_B6.TIF',
       'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized'},
      'extensions': [<Band name=OLI_B6>]},
     {'name': 'swir22',
      'asset_fields': {'href': 'https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC08_L2SP_014032_20221219_20230113_02_T1/LC08_L2SP_014032_20221219_20230113_02_T1_SR_B7.TIF',
       'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized'},
      'extensions': [<Band name=OLI_B7>]},
     {'name': 'lwir11',
      'asset_fields': {'href': 'https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC08_L2SP_014032_20221219_20230113_02_T1/LC08_L2SP_014032_20221219_20230113_02_T1_ST_B10.TIF',
       'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized'},
      'extensions': [<Band name=TIRS_B10>]},
     {'name': 'trad',
      'asset_fields': {'href': 'https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC08_L2SP_014032_20221219_20230113_02_T1/LC08_L2SP_014032_20221219_20230113_02_T1_ST_TRAD.TIF',
       'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized',
       'title': 'Thermal radiance'},
      'extensions': [<Raster Band>]},
     {'name': 'urad',
      'asset_fields': {'href': 'https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC08_L2SP_014032_20221219_20230113_02_T1/LC08_L2SP_014032_20221219_20230113_02_T1_ST_URAD.TIF',
       'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized',
       'title': 'Upwelled radiance'},
      'extensions': [<Raster Band>]},
     {'name': 'drad',
      'asset_fields': {'href': 'https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC08_L2SP_014032_20221219_20230113_02_T1/LC08_L2SP_014032_20221219_20230113_02_T1_ST_DRAD.TIF',
       'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized',
       'title': 'Downwelled radiance'},
      'extensions': [<Raster Band>]},
     {'name': 'emis',
      'asset_fields': {'href': 'https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC08_L2SP_014032_20221219_20230113_02_T1/LC08_L2SP_014032_20221219_20230113_02_T1_ST_EMIS.TIF',
       'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized',
       'title': 'Emissivity'},
      'extensions': [<Raster Band>]},
     {'name': 'emsd',
      'asset_fields': {'href': 'https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC08_L2SP_014032_20221219_20230113_02_T1/LC08_L2SP_014032_20221219_20230113_02_T1_ST_EMSD.TIF',
       'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized',
       'title': 'Emissivity standard deviation'},
      'extensions': [<Raster Band>]},
     {'name': 'atran',
      'asset_fields': {'href': 'https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC08_L2SP_014032_20221219_20230113_02_T1/LC08_L2SP_014032_20221219_20230113_02_T1_ST_ATRAN.TIF',
       'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized',
       'title': 'Atmospheric transmission'},
      'extensions': [<Raster Band>]},
     {'name': 'cdist',
      'asset_fields': {'href': 'https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC08_L2SP_014032_20221219_20230113_02_T1/LC08_L2SP_014032_20221219_20230113_02_T1_ST_CDIST.TIF',
       'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized',
       'title': 'Cloud distance'},
      'extensions': [<Raster Band>]},
     {'name': 'qa',
      'asset_fields': {'href': 'https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC08_L2SP_014032_20221219_20230113_02_T1/LC08_L2SP_014032_20221219_20230113_02_T1_ST_QA.TIF',
       'title': 'Surface Temperature Quality Assessment Band'},
      'extensions': [<Raster Band>]},
     {'name': 'qa_pixel',
      'asset_fields': {'href': 'https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC08_L2SP_014032_20221219_20230113_02_T1/LC08_L2SP_014032_20221219_20230113_02_T1_QA_PIXEL.TIF',
       'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized',
       'title': 'Pixel quality assessment'},
      'extensions': [[<Bitfield offset=0 length=1 classes=[<Classification value=0 description=Image data>, <Classification value=1 description=Fill data>]>,
        <Bitfield offset=1 length=1 classes=[<Classification value=0 description=Cloud is not dilated or no cloud>, <Classification value=1 description=Dilated cloud>]>,
        <Bitfield offset=2 length=1 classes=[<Classification value=0 description=Cirrus confidence is not high>, <Classification value=1 description=High confidence cirrus>]>,
        <Bitfield offset=3 length=1 classes=[<Classification value=0 description=Cloud confidence is not high>, <Classification value=1 description=High confidence cloud>]>,
        <Bitfield offset=4 length=1 classes=[<Classification value=0 description=Cloud shadow confidence is not high>, <Classification value=1 description=High confidence cloud shadow>]>,
        <Bitfield offset=5 length=1 classes=[<Classification value=0 description=Snow confidence is not high>, <Classification value=1 description=High confidence snow>]>,
        <Bitfield offset=6 length=1 classes=[<Classification value=0 description=Clear confidence is not high>, <Classification value=1 description=High confidence clear>]>,
        <Bitfield offset=7 length=1 classes=[<Classification value=0 description=Water confidence is not high>, <Classification value=1 description=High confidence water>]>,
        <Bitfield offset=8 length=2 classes=[<Classification value=0 description=No confidence level set>, <Classification value=1 description=Low confidence cloud>, <Classification value=2 description=Medium confidence cloud>, <Classification value=3 description=High confidence cloud>]>,
        <Bitfield offset=10 length=2 classes=[<Classification value=0 description=No confidence level set>, <Classification value=1 description=Low confidence cloud shadow>, <Classification value=2 description=Reserved - value not used>, <Classification value=3 description=High confidence cloud shadow>]>,
        <Bitfield offset=12 length=2 classes=[<Classification value=0 description=No confidence level set>, <Classification value=1 description=Low confidence snow>, <Classification value=2 description=Reserved - value not used>, <Classification value=3 description=High confidence snow>]>,
        <Bitfield offset=14 length=2 classes=[<Classification value=0 description=No confidence level set>, <Classification value=1 description=Low confidence cirrus>, <Classification value=2 description=Reserved - value not used>, <Classification value=3 description=High confidence cirrus>]>]]},
     {'name': 'qa_radsat',
      'asset_fields': {'href': 'https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC08_L2SP_014032_20221219_20230113_02_T1/LC08_L2SP_014032_20221219_20230113_02_T1_QA_RADSAT.TIF',
       'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized',
       'description': 'Collection 2 Level-1 Radiometric Saturation and Terrain Occlusion Quality Assessment Band (QA_RADSAT)'},
      'extensions': [[<Bitfield offset=0 length=1 classes=[<Classification value=0 description=Band 1 not saturated>, <Classification value=1 description=Band 1 saturated>]>,
        <Bitfield offset=1 length=1 classes=[<Classification value=0 description=Band 2 not saturated>, <Classification value=1 description=Band 2 saturated>]>,
        <Bitfield offset=2 length=1 classes=[<Classification value=0 description=Band 3 not saturated>, <Classification value=1 description=Band 3 saturated>]>,
        <Bitfield offset=3 length=1 classes=[<Classification value=0 description=Band 4 not saturated>, <Classification value=1 description=Band 4 saturated>]>,
        <Bitfield offset=4 length=1 classes=[<Classification value=0 description=Band 5 not saturated>, <Classification value=1 description=Band 5 saturated>]>,
        <Bitfield offset=5 length=1 classes=[<Classification value=0 description=Band 6 not saturated>, <Classification value=1 description=Band 6 saturated>]>,
        <Bitfield offset=6 length=1 classes=[<Classification value=0 description=Band 7 not saturated>, <Classification value=1 description=Band 7 saturated>]>,
        <Bitfield offset=8 length=1 classes=[<Classification value=0 description=Band 9 not saturated>, <Classification value=1 description=Band 9 saturated>]>,
        <Bitfield offset=11 length=1 classes=[<Classification value=0 description=Terrain is not occluded>, <Classification value=1 description=Terrain is occluded>]>]]},
     {'name': 'qa_aerosol',
      'asset_fields': {'href': 'https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC08_L2SP_014032_20221219_20230113_02_T1/LC08_L2SP_014032_20221219_20230113_02_T1_SR_QA_AEROSOL.TIF',
       'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized',
       'title': 'Aerosol Quality Assessment Band'},
      'extensions': [[<Bitfield offset=0 length=1 classes=[<Classification value=0 description=Pixel is not fill>, <Classification value=1 description=Pixel is fill>]>,
        <Bitfield offset=1 length=1 classes=[<Classification value=0 description=Pixel retrieval is not valid>, <Classification value=1 description=Pixel retrieval is valid>]>,
        <Bitfield offset=2 length=1 classes=[<Classification value=0 description=Pixel is not water>, <Classification value=1 description=Pixel is water>]>,
        <Bitfield offset=5 length=1 classes=[<Classification value=0 description=Pixel is not interpolated>, <Classification value=1 description=Pixel is interpolated>]>,
        <Bitfield offset=6 length=2 classes=[<Classification value=0 description=No aerosol correction applied>, <Classification value=1 description=Low aerosol level>, <Classification value=2 description=Medium aerosol level>, <Classification value=3 description=High aerosol level>]>]]}]

With this information we can now define a method that adds all the relevant assets for a scene to an item:

    [29]:

    def add_assets(item: pystac.Item, band_info: List[BandInfo]) -> None:
        for band in band_info:
            asset = pystac.Asset(**band["asset_fields"])
            asset.set_owner(item)
            for ext_data in band["extensions"]:
                if isinstance(ext_data, EOBand):
                    EOExtension.ext(asset, add_if_missing=True).bands = [ext_data]
                elif isinstance(ext_data, RasterBand):
                    RasterExtension.ext(asset, add_if_missing=True).bands = [ext_data]
                elif isinstance(ext_data, list):
                    ClassificationExtension.ext(
                        asset, add_if_missing=True
                    ).bitfields = ext_data
            item.add_asset(band["name"], asset)

    [30]:

    add_assets(item, landsat_band_info(metadata, download_url))

We can examine the item to ensure that the assets appear as expected.

    [31]:

    item.assets["nir08"].to_dict()

    [31]:

    {'href': 'https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC08_L2SP_014032_20221219_20230113_02_T1/LC08_L2SP_014032_20221219_20230113_02_T1_SR_B5.TIF',
     'type': 'image/tiff; application=geotiff; profile=cloud-optimized',
     'eo:bands': [{'name': 'OLI_B5',
       'common_name': 'nir08',
       'description': 'Near infrared (Operational Land Imager)',
       'center_wavelength': 0.87,
       'full_width_half_max': 0.03}]}

    [32]:

    item.assets["qa_aerosol"].to_dict()

    [32]:

    {'href': 'https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC08_L2SP_014032_20221219_20230113_02_T1/LC08_L2SP_014032_20221219_20230113_02_T1_SR_QA_AEROSOL.TIF',
     'type': 'image/tiff; application=geotiff; profile=cloud-optimized',
     'title': 'Aerosol Quality Assessment Band',
     'classification:bitfields': [{'offset': 0,
       'length': 1,
       'classes': [{'value': 0,
         'name': 'not_fill',
         'description': 'Pixel is not fill'},
        {'value': 1, 'name': 'fill', 'description': 'Pixel is fill'}],
       'description': 'Image or fill data',
       'name': 'fill'},
      {'offset': 1,
       'length': 1,
       'classes': [{'value': 0,
         'name': 'not_valid',
         'description': 'Pixel retrieval is not valid'},
        {'value': 1, 'name': 'valid', 'description': 'Pixel retrieval is valid'}],
       'description': 'Valid aerosol retrieval',
       'name': 'retrieval'},
      {'offset': 2,
       'length': 1,
       'classes': [{'value': 0,
         'name': 'not_water',
         'description': 'Pixel is not water'},
        {'value': 1, 'name': 'water', 'description': 'Pixel is water'}],
       'description': 'Water mask',
       'name': 'water'},
      {'offset': 5,
       'length': 1,
       'classes': [{'value': 0,
         'name': 'not_interpolated',
         'description': 'Pixel is not interpolated'},
        {'value': 1,
         'name': 'interpolated',
         'description': 'Pixel is interpolated'}],
       'description': 'Aerosol interpolation',
       'name': 'interpolated'},
      {'offset': 6,
       'length': 2,
       'classes': [{'value': 0,
         'name': 'climatology',
         'description': 'No aerosol correction applied'},
        {'value': 1, 'name': 'low', 'description': 'Low aerosol level'},
        {'value': 2, 'name': 'medium', 'description': 'Medium aerosol level'},
        {'value': 3, 'name': 'high', 'description': 'High aerosol level'}],
       'description': 'Aerosol level',
       'name': 'level'}]}

## Add projection information

We can specify the EPSG code for the scene as part of the projection extension. The below method, adapted from stactools, figures out the correct UTM Zone EPSG:

    [33]:

    def get_epsg(metadata: Dict[str, Any], min_lat: float, max_lat: float) -> Optional[int]:
        if "UTM_ZONE" in metadata["PROJECTION_ATTRIBUTES"]:
            utm_zone_integer = metadata["PROJECTION_ATTRIBUTES"]["UTM_ZONE"].zfill(2)
            return int(f"326{utm_zone_integer}")
        else:
            lat_ts = metadata["PROJECTION_ATTRIBUTES"]["TRUE_SCALE_LAT"]
            if lat_ts == "-71.00000":
                # Antarctic
                return 3031
            elif lat_ts == "71.00000":
                # Arctic
                return 3995
            else:
                raise ValueError(
                    f"Unexpected value for PROJECTION_ATTRIBUTES/TRUE_SCALE_LAT: {lat_ts} "
                )

    [34]:

    proj_ext = ProjectionExtension.ext(item, add_if_missing=True)
    assert item.bbox is not None
    proj_ext.epsg = get_epsg(metadata, item.bbox[1], item.bbox[3])
    proj_ext.epsg

    [34]:

    32618

## Add view geometry information

The View Geometry extension specifies information related to angles of sensors and other radiance angles that affect the view of resulting data. The Landsat metadata specifies two of these parameters, so we add them to our Item:

    [35]:

    def get_view_info(metadata: Dict[str, Any]) -> Dict[str, float]:
        return {
            "sun_azimuth": float(metadata["IMAGE_ATTRIBUTES"]["SUN_AZIMUTH"]),
            "sun_elevation": float(metadata["IMAGE_ATTRIBUTES"]["SUN_ELEVATION"]),
        }

    [36]:

    view_ext = ViewExtension.ext(item, add_if_missing=True)
    view_info = get_view_info(metadata)
    view_ext.sun_azimuth = view_info["sun_azimuth"]
    view_ext.sun_elevation = view_info["sun_elevation"]
    item.properties

    [36]:

    {'datetime': '2022-12-19T15:40:17.729916Z',
     'gsd': 30.0,
     'eo:cloud_cover': 43.42,
     'proj:epsg': 32618,
     'view:sun_azimuth': 160.86021018,
     'view:sun_elevation': 23.81656674}

Now that we’ve added all the metadata to the Item, let’s check the validator to make sure we’ve specified everything correctly. The validation logic will take into account the new extensions that have been enabled and validate against the proper schemas for those extensions.

    [37]:

    item.validate()

    [37]:

    ['https://schemas.stacspec.org/v1.0.0/item-spec/json-schema/item.json',
     'https://stac-extensions.github.io/eo/v1.1.0/schema.json',
     'https://stac-extensions.github.io/raster/v1.1.0/schema.json',
     'https://stac-extensions.github.io/classification/v1.1.0/schema.json',
     'https://stac-extensions.github.io/projection/v1.1.0/schema.json',
     'https://stac-extensions.github.io/view/v1.0.0/schema.json']

### Building the Collection

Now that we know how to build an Item for a scene, let’s build the Collection that will contain all the Items.

If we look at the `__init__` method for `pystac.Collection`, we can see what properties are required:

    [38]:

    help(pystac.Collection.__init__)

    Help on function __init__ in module pystac.collection:

    __init__(self, id: 'str', description: 'str', extent: 'Extent', title: 'Optional[str]' = None, stac_extensions: 'Optional[List[str]]' = None, href: 'Optional[str]' = None, extra_fields: 'Optional[Dict[str, Any]]' = None, catalog_type: 'Optional[CatalogType]' = None, license: 'str' = 'proprietary', keywords: 'Optional[List[str]]' = None, providers: 'Optional[List[Provider]]' = None, summaries: 'Optional[Summaries]' = None, assets: 'Optional[Dict[str, Asset]]' = None)
        Initialize self.  See help(type(self)) for accurate signature.

#### Collection `id`

We’ll use the location name we defined above in the ID for our Collection:

    [39]:

    collection_id = "{}-landsat-collection-2".format(location_name.lower())
    collection_id

    [39]:

    'philly-landsat-collection-2'

#### Collection `title`

Here we set a simple title for our collection.

    [40]:

    collection_title = "2022 Landsat images over {}".format(location_name.lower())
    collection_title

    [40]:

    '2022 Landsat images over philly'

#### Collection `description`

Here we give a brief description of the Collection. If this were a real Collection that were being published, I’d recommend putting a much more detailed description to ensure anyone using your STAC knows what they are working with!

Notice we are using Markdown to write the description. The `description` field can be Markdown to help tools that render information about STAC to display the information in a more readable way.

    [41]:

    collection_description = """### {} Landsat Collection 2

    A collection of Landsat scenes around {} in 2022.
    """.format(location_name, location_name)
    print(collection_description)

    ### Philly Landsat Collection 2

    A collection of Landsat scenes around Philly in 2022.

#### Collection `extent`

A Collection specifies the spatial and temporal extent of all the item it contains. Since Landsat spans the globe, we’ll simply put a global extent here. We’ll also specify an open-ended time interval.

Towards the end of the notebook, we’ll use a method to easily scope this down to cover the times and space the Items occupy once we’ve added all the items.

    [42]:

    spatial_extent = pystac.SpatialExtent([[-180.0, -90.0, 180.0, 90.0]])
    temporal_extent = pystac.TemporalExtent([[datetime(2013, 6, 1), None]])
    collection_extent = pystac.Extent(spatial_extent, temporal_extent)

    [43]:

    collection = pystac.Collection(
        id=collection_id,
        title=collection_title,
        description=collection_description,
        extent=collection_extent,
    )

We can now look at our Collection as a `dict` to check our values.

    [44]:

    collection.to_dict()

    [44]:

    {'type': 'Collection',
     'id': 'philly-landsat-collection-2',
     'stac_version': '1.0.0',
     'description': '### Philly Landsat Collection 2\n\nA collection of Landsat scenes around Philly in 2022.\n',
     'links': [],
     'title': '2022 Landsat images over philly',
     'extent': {'spatial': {'bbox': [[-180.0, -90.0, 180.0, 90.0]]},
      'temporal': {'interval': [['2013-06-01T00:00:00Z', None]]}},
     'license': 'proprietary'}

#### Set the license

Notice the `license` above is `proprietary`. This is the default in PySTAC if no license is specified; however Landsat is certainly not proprietary (thankfully!), so let’s change the license to the correct SPDX string for public domain data:

    [45]:

    collection_license = "PDDL-1.0"

#### Set the providers

A collection will specify the providers of the data, including what role they have played. We can set our provider information by instantiating `pystac.Provider` objects:

    [46]:

    collection.providers = [
        pystac.Provider(
            name="NASA",
            roles=[pystac.ProviderRole.PRODUCER, pystac.ProviderRole.LICENSOR],
            url="https://landsat.gsfc.nasa.gov/",
        ),
        pystac.Provider(
            name="USGS",
            roles=[
                pystac.ProviderRole.PROCESSOR,
                pystac.ProviderRole.PRODUCER,
                pystac.ProviderRole.LICENSOR,
            ],
            url="https://www.usgs.gov/landsat-missions/landsat-collection-2-level-2-science-products",
        ),
        pystac.Provider(
            name="Microsoft",
            roles=[pystac.ProviderRole.HOST],
            url="https://planetarycomputer.microsoft.com",
        ),
    ]

## Create items for each scene

We created an Item for a single scene above. This method consolidates that logic into a single method that can construct an Item from a scene, so we can create an Item for every scene in our subset:

    [47]:

    def item_from_metadata(mtl_xml_url: str) -> pystac.Item:
        metadata = get_metadata(mtl_xml_url)
        download_url = partial(download_sidecar, metadata)

        bbox = get_bbox(metadata)
        item = pystac.Item(
            id=get_item_id(metadata),
            datetime=get_datetime(metadata),
            geometry=get_geometry(metadata, bbox),
            bbox=bbox,
            properties={},
        )

        item.common_metadata.gsd = 30.0

        item_eo_ext = EOExtension.ext(item, add_if_missing=True)
        item_eo_ext.cloud_cover = get_cloud_cover(metadata)

        add_assets(item, landsat_band_info(metadata, download_url))

        item_proj_ext = ProjectionExtension.ext(item, add_if_missing=True)
        assert item.bbox is not None
        item_proj_ext.epsg = get_epsg(metadata, item.bbox[1], item.bbox[3])

        item_view_ext = ViewExtension.ext(item, add_if_missing=True)
        view_info = get_view_info(metadata)
        item_view_ext.sun_azimuth = view_info["sun_azimuth"]
        item_view_ext.sun_elevation = view_info["sun_elevation"]

        item.validate()

        return item

Here we create an item per scene and add it to our collection. Since this is reading multiple metadata files per scene from the internet, it may take a little bit to run:

    [48]:

    for url in scene_mtls:
        try:
            item = item_from_metadata(url)
            collection.add_item(item)
        except Exception as e:
            print(e)

    ANG file for url https://landsateuwest.blob.core.windows.net/landsat-c2/level-2/standard/oli-tirs/2022/014/032/LC09_L2SP_014032_20221024_20221026_02_T2/LC09_L2SP_014032_20221024_20221026_02_T2_ANG.txt is incorrectly formatted

## Reset collection extent based on items

Now that we’ve added all the item we can use the `update_extent_from_items` method on the Collection to set the extent based on the contained items:

    [49]:

    collection.update_extent_from_items()
    collection.extent.to_dict()

    [49]:

    {'spatial': {'bbox': [[-76.29048, 39.25502, -73.48833, 41.38441]]},
     'temporal': {'interval': [['2022-10-08T15:40:14.577173Z',
        '2022-12-19T15:40:17.729916Z']]}}

## Set the HREFs for everything in the catalog

We’ve been building up our Collection and Items in memory. This has been convenient as it allows us not to think about file paths as we construct our Catalog. However, a STAC is not valid without any HREFs!

We can use the `normalize_hrefs` method to set all the HREFs in the entire STAC based on a root directory. This will use the STAC Best Practices recommendations for STAC file layout for each Catalog, Collection and Item in the STAC.

Here we use that method and set the root directory to a subdirectory of our user’s `home` directory:

    [50]:

    from pathlib import Path

    root_path = str(Path.home() / "{}-landsat-stac".format(location_name))

    collection.normalize_hrefs(root_path)

Now that we have all the Collection’s data set and HREFs in place we can validate the entire STAC using `validate_all`, which recursively crawls through a catalog and validates every STAC object in the catalog:

    [51]:

    collection.validate_all()

    [51]:

    8

## Write the catalog locally

Now that we have our complete, validated STAC in memory, let’s write it out. This is as simple as calling `save` on the Collection. We need to specify the type of catalog in order to property write out links - these types are described again in the STAC Best Practices documentation.

We’ll use the “self contained” type, which uses relative paths and does not specify absolute “self” links to any object. This makes the catalog more portable, as it remains valid even if you copy the STAC to new locations.

    [52]:

    collection.save(pystac.CatalogType.SELF_CONTAINED)

Now that we’ve written our STAC out we probably want to view it. We can use the `describe` method to print out a simple representation of the catalog:

    [53]:

    collection.describe()

    * <Collection id=philly-landsat-collection-2>
      * <Item id=LC80140322022353>
      * <Item id=LC90140322022345>
      * <Item id=LC80140322022337>
      * <Item id=LC90140322022329>
      * <Item id=LC80140322022321>
      * <Item id=LC90140322022313>
      * <Item id=LC80140322022305>
      * <Item id=LC90140322022281>

We can also use the `to_dict` method on individual STAC objects in order to see the data, as we’ve been doing in the tutorial:

    [54]:

    collection.to_dict()

    [54]:

    {'type': 'Collection',
     'id': 'philly-landsat-collection-2',
     'stac_version': '1.0.0',
     'description': '### Philly Landsat Collection 2\n\nA collection of Landsat scenes around Philly in 2022.\n',
     'links': [{'rel': 'root',
       'href': './collection.json',
       'type': 'application/json',
       'title': '2022 Landsat images over philly'},
      {'rel': 'item',
       'href': './LC80140322022353/LC80140322022353.json',
       'type': 'application/json'},
      {'rel': 'item',
       'href': './LC90140322022345/LC90140322022345.json',
       'type': 'application/json'},
      {'rel': 'item',
       'href': './LC80140322022337/LC80140322022337.json',
       'type': 'application/json'},
      {'rel': 'item',
       'href': './LC90140322022329/LC90140322022329.json',
       'type': 'application/json'},
      {'rel': 'item',
       'href': './LC80140322022321/LC80140322022321.json',
       'type': 'application/json'},
      {'rel': 'item',
       'href': './LC90140322022313/LC90140322022313.json',
       'type': 'application/json'},
      {'rel': 'item',
       'href': './LC80140322022305/LC80140322022305.json',
       'type': 'application/json'},
      {'rel': 'item',
       'href': './LC90140322022281/LC90140322022281.json',
       'type': 'application/json'},
      {'rel': 'self',
       'href': '/Users/pjh/Philly-landsat-stac/collection.json',
       'type': 'application/json'}],
     'title': '2022 Landsat images over philly',
     'extent': {'spatial': {'bbox': [[-76.29048, 39.25502, -73.48833, 41.38441]]},
      'temporal': {'interval': [['2022-10-08T15:40:14.577173Z',
         '2022-12-19T15:40:17.729916Z']]}},
     'license': 'proprietary',
     'providers': [{'name': 'NASA',
       'roles': [<ProviderRole.PRODUCER: 'producer'>,
        <ProviderRole.LICENSOR: 'licensor'>],
       'url': 'https://landsat.gsfc.nasa.gov/'},
      {'name': 'USGS',
       'roles': [<ProviderRole.PROCESSOR: 'processor'>,
        <ProviderRole.PRODUCER: 'producer'>,
        <ProviderRole.LICENSOR: 'licensor'>],
       'url': 'https://www.usgs.gov/landsat-missions/landsat-collection-2-level-2-science-products'},
      {'name': 'Microsoft',
       'roles': [<ProviderRole.HOST: 'host'>],
       'url': 'https://planetarycomputer.microsoft.com'}]}

However, if we want to browse our STAC more interactively, we can serve the Collection from a local webserver and then browse the Collection with stac-browser.

We can use this simple Python server (copied from this gist) to serve our our directory at port 5555:

    [ ]:

    import os
    from http.server import HTTPServer, SimpleHTTPRequestHandler

    os.chdir(root_path)

    class CORSRequestHandler(SimpleHTTPRequestHandler):
        def end_headers(self) -> None:
            self.send_header("Access-Control-Allow-Origin", "*")
            self.send_header("Access-Control-Allow-Methods", "GET")
            self.send_header("Cache-Control", "no-store, no-cache, must-revalidate")
            return super(CORSRequestHandler, self).end_headers()

    with HTTPServer(("localhost", 5555), CORSRequestHandler) as httpd:
        httpd.serve_forever()

Now we can browse our STAC Collection with stac-browser in a few different ways: 1. Follow the instructions for starting a stac-browser instance and point it at `http://localhost:5555/collection.json` to serve out the STAC. 2. If you want to avoid setting up your own stac-browser instance, you can use the STAC Browser Demo hosted by Radiant Earth: https://radiantearth.github.io/stac-browser/#/http://localhost:5555/collection.json

To quit the server, use the `Kernel` -> `Interrupt` menu option.

## Acknowledgements

Credit to sat-stac-landsat from which a lot of this code was based.