{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Creating a STAC of Landsat data\n", "\n", "In this tutorial we create a STAC of Landsat data provided by Microsoft's [Planetary Computer](https://planetarycomputer.microsoft.com/dataset/landsat-c2-l2). 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](https://github.com/radiantearth/stac-browser) to preview the images.\n", "\n", "### Requirements\n", "\n", "To run this tutorial you'll need to have installed PySTAC with the validation extra and the Planetary Computer package. To do this, use:\n", "\n", "```\n", "pip install 'pystac[validation]' planetary-computer\n", "```\n", "\n", "Also to run this notebook you'll need [jupyter](https://jupyter.org/) 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." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import json\n", "from datetime import datetime\n", "from functools import partial\n", "from os.path import dirname, join\n", "from typing import Any, Callable, Dict, List, Optional, Tuple, Union, cast\n", "from urllib.parse import urlparse, urlunparse\n", "\n", "import planetary_computer as pc\n", "from dateutil.parser import parse\n", "from typing_extensions import TypedDict\n", "\n", "import pystac\n", "from pystac.extensions.classification import (\n", " Bitfield,\n", " Classification,\n", " ClassificationExtension,\n", ")\n", "from pystac.extensions.eo import Band as EOBand\n", "from pystac.extensions.eo import EOExtension\n", "from pystac.extensions.projection import ProjectionExtension\n", "from pystac.extensions.raster import RasterBand, RasterExtension\n", "from pystac.extensions.view import ViewExtension" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Identify target scenes\n", "\n", "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:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "location_name = \"Philly\"\n", "scene_mtls = [\n", " \"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\",\n", " \"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\",\n", " \"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\",\n", " \"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\",\n", " \"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\",\n", " \"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\",\n", " \"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\",\n", " \"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\",\n", " \"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\",\n", "]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Read metadata from the MTL file\n", "\n", "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:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Taken from https://stackoverflow.com/questions/2148119/how-to-convert-an-xml-string-to-a-dictionary\n", "from xml.etree import cElementTree as ElementTree\n", "\n", "\n", "class XmlListConfig(list):\n", " def __init__(self, aList):\n", " for element in aList:\n", " if element:\n", " if len(element) == 1 or element[0].tag != element[1].tag:\n", " self.append(XmlDictConfig(element))\n", " elif element[0].tag == element[1].tag:\n", " self.append(XmlListConfig(element))\n", " elif element.text:\n", " text = element.text.strip()\n", " if text:\n", " self.append(text)\n", "\n", "\n", "class XmlDictConfig(dict):\n", " def __init__(self, parent_element):\n", " if parent_element.items():\n", " self.update(dict(parent_element.items()))\n", " for element in parent_element:\n", " if element:\n", " if len(element) == 1 or element[0].tag != element[1].tag:\n", " aDict = XmlDictConfig(element)\n", " else:\n", " aDict = {element[0].tag: XmlListConfig(element)}\n", " if element.items():\n", " aDict.update(dict(element.items()))\n", " self.update({element.tag: aDict})\n", " elif element.items():\n", " self.update({element.tag: dict(element.items())})\n", " else:\n", " self.update({element.tag: element.text})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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](https://pystac.readthedocs.io/en/latest/concepts.html#using-stac-io). For now we'll use it directly as an easy way to read a text file from an HTTPS source." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "stac_io = pystac.StacIO.default()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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](https://github.com/microsoft/planetary-computer-sdk-for-python) 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." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def get_metadata(xml_url: str) -> Dict[str, Any]:\n", " result = XmlDictConfig(ElementTree.XML(stac_io.read_text(pc.sign(xml_url))))\n", " result[\"ORIGINAL_URL\"] = (\n", " xml_url # Include the original URL in the metadata for use later\n", " )\n", " return result" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's read the MTL file for the first scene and see what it looks like." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{\n", " \"PRODUCT_CONTENTS\": {\n", " \"ORIGIN\": \"Image courtesy of the U.S. Geological Survey\",\n", " \"DIGITAL_OBJECT_IDENTIFIER\": \"https://doi.org/10.5066/P9OGBGM6\",\n", " \"LANDSAT_PRODUCT_ID\": \"LC08_L2SP_014032_20221219_20230113_02_T1\",\n", " \"PROCESSING_LEVEL\": \"L2SP\",\n", " \"COLLECTION_NUMBER\": \"02\",\n", " \"COLLECTION_CATEGORY\": \"T1\",\n", " \"OUTPUT_FORMAT\": \"GEOTIFF\",\n", " \"FILE_NAME_BAND_1\": \"LC08_L2SP_014032_20221219_20230113_02_T1_SR_B1.TIF\",\n", " \"FILE_NAME_BAND_2\": \"LC08_L2SP_014032_20221219_20230113_02_T1_SR_B2.TIF\",\n", " \"FILE_NAME_BAND_3\": \"LC08_L2SP_014032_20221219_20230113_02_T1_SR_B3.TIF\",\n", " \"FILE_NAME_BAND_4\": \"LC08_L2SP_014032_20221219_20230113_02_T1_SR_B4.TIF\",\n", " \"FILE_NAME_BAND_5\": \"LC08_L2SP_014032_20221219_20230113_02_T1_SR_B5.TIF\",\n", " \"FILE_NAME_BAND_6\": \"LC08_L2SP_014032_20221219_20230113_02_T1_SR_B6.TIF\",\n", " \"FILE_NAME_BAND_7\": \"LC08_L2SP_014032_20221219_20230113_02_T1_SR_B7.TIF\",\n", " \"FILE_NAME_BAND_ST_B10\": \"LC08_L2SP_014032_20221219_20230113_02_T1_ST_B10.TIF\",\n", " \"FILE_NAME_THERMAL_RADIANCE\": \"LC08_L2SP_014032_20221219_20230113_02_T1_ST_TRAD.TIF\",\n", " \"FILE_NAME_UPWELL_RADIANCE\": \"LC08_L2SP_014032_20221219_20230113_02_T1_ST_URAD.TIF\",\n", " \"FILE_NAME_DOWNWELL_RADIANCE\": \"LC08_L2SP_014032_20221219_20230113_02_T1_ST_DRAD.TIF\",\n", " \"FILE_NAME_ATMOSPHERIC_TRANSMITTANCE\": \"LC08_L2SP_014032_20221219_20230113_02_T1_ST_ATRAN.TIF\",\n", " \"FILE_NAME_EMISSIVITY\": \"LC08_L2SP_014032_20221219_20230113_02_T1_ST_EMIS.TIF\",\n", " \"FILE_NAME_EMISSIVITY_STDEV\": \"LC08_L2SP_014032_20221219_20230113_02_T1_ST_EMSD.TIF\",\n", " \"FILE_NAME_CLOUD_DISTANCE\": \"LC08_L2SP_014032_20221219_20230113_02_T1_ST_CDIST.TIF\",\n", " \"FILE_NAME_QUALITY_L2_AEROSOL\": \"LC08_L2SP_014032_20221219_20230113_02_T1_SR_QA_AEROSOL.TIF\",\n", " \"FILE_NAME_QUALITY_L2_SURFACE_TEMPERATURE\": \"LC08_L2SP_014032_20221219_20230113_02_T1_ST_QA.TIF\",\n", " \"FILE_NAME_QUALITY_L1_PIXEL\": \"LC08_L2SP_014032_20221219_20230113_02_T1_QA_PIXEL.TIF\",\n", " \"FILE_NAME_QUALITY_L1_RADIOMETRIC_SATURATION\": \"LC08_L2SP_014032_20221219_20230113_02_T1_QA_RADSAT.TIF\",\n", " \"FILE_NAME_ANGLE_COEFFICIENT\": \"LC08_L2SP_014032_20221219_20230113_02_T1_ANG.txt\",\n", " \"FILE_NAME_METADATA_ODL\": \"LC08_L2SP_014032_20221219_20230113_02_T1_MTL.txt\",\n", " \"FILE_NAME_METADATA_XML\": \"LC08_L2SP_014032_20221219_20230113_02_T1_MTL.xml\",\n", " \"DATA_TYPE_BAND_1\": \"UINT16\",\n", " \"DATA_TYPE_BAND_2\": \"UINT16\",\n", " \"DATA_TYPE_BAND_3\": \"UINT16\",\n", " \"DATA_TYPE_BAND_4\": \"UINT16\",\n", " \"DATA_TYPE_BAND_5\": \"UINT16\",\n", " \"DATA_TYPE_BAND_6\": \"UINT16\",\n", " \"DATA_TYPE_BAND_7\": \"UINT16\",\n", " \"DATA_TYPE_BAND_ST_B10\": \"UINT16\",\n", " \"DATA_TYPE_THERMAL_RADIANCE\": \"INT16\",\n", " \"DATA_TYPE_UPWELL_RADIANCE\": \"INT16\",\n", " \"DATA_TYPE_DOWNWELL_RADIANCE\": \"INT16\",\n", " \"DATA_TYPE_ATMOSPHERIC_TRANSMITTANCE\": \"INT16\",\n", " \"DATA_TYPE_EMISSIVITY\": \"INT16\",\n", " \"DATA_TYPE_EMISSIVITY_STDEV\": \"INT16\",\n", " \"DATA_TYPE_CLOUD_DISTANCE\": \"INT16\",\n", " \"DATA_TYPE_QUALITY_L2_AEROSOL\": \"UINT8\",\n", " \"DATA_TYPE_QUALITY_L2_SURFACE_TEMPERATURE\": \"INT16\",\n", " \"DATA_TYPE_QUALITY_L1_PIXEL\": \"UINT16\",\n", " \"DATA_TYPE_QUALITY_L1_RADIOMETRIC_SATURATION\": \"UINT16\"\n", " },\n", " \"IMAGE_ATTRIBUTES\": {\n", " \"SPACECRAFT_ID\": \"LANDSAT_8\",\n", " \"SENSOR_ID\": \"OLI_TIRS\",\n", " \"WRS_TYPE\": \"2\",\n", " \"WRS_PATH\": \"14\",\n", " \"WRS_ROW\": \"32\",\n", " \"NADIR_OFFNADIR\": \"NADIR\",\n", " \"TARGET_WRS_PATH\": \"14\",\n", " \"TARGET_WRS_ROW\": \"32\",\n", " \"DATE_ACQUIRED\": \"2022-12-19\",\n", " \"SCENE_CENTER_TIME\": \"15:40:17.7299160Z\",\n", " \"STATION_ID\": \"LGN\",\n", " \"CLOUD_COVER\": \"43.42\",\n", " \"CLOUD_COVER_LAND\": \"48.41\",\n", " \"IMAGE_QUALITY_OLI\": \"9\",\n", " \"IMAGE_QUALITY_TIRS\": \"9\",\n", " \"SATURATION_BAND_1\": \"N\",\n", " \"SATURATION_BAND_2\": \"Y\",\n", " \"SATURATION_BAND_3\": \"N\",\n", " \"SATURATION_BAND_4\": \"Y\",\n", " \"SATURATION_BAND_5\": \"Y\",\n", " \"SATURATION_BAND_6\": \"Y\",\n", " \"SATURATION_BAND_7\": \"Y\",\n", " \"SATURATION_BAND_8\": \"N\",\n", " \"SATURATION_BAND_9\": \"N\",\n", " \"ROLL_ANGLE\": \"-0.001\",\n", " \"SUN_AZIMUTH\": \"160.86021018\",\n", " \"SUN_ELEVATION\": \"23.81656674\",\n", " \"EARTH_SUN_DISTANCE\": \"0.9839500\",\n", " \"TRUNCATION_OLI\": \"UPPER\",\n", " \"TIRS_SSM_MODEL\": \"FINAL\",\n", " \"TIRS_SSM_POSITION_STATUS\": \"ESTIMATED\"\n", " },\n", " \"PROJECTION_ATTRIBUTES\": {\n", " \"MAP_PROJECTION\": \"UTM\",\n", " \"DATUM\": \"WGS84\",\n", " \"ELLIPSOID\": \"WGS84\",\n", " \"UTM_ZONE\": \"18\",\n", " \"GRID_CELL_SIZE_REFLECTIVE\": \"30.00\",\n", " \"GRID_CELL_SIZE_THERMAL\": \"30.00\",\n", " \"REFLECTIVE_LINES\": \"7861\",\n", " \"REFLECTIVE_SAMPLES\": \"7731\",\n", " \"THERMAL_LINES\": \"7861\",\n", " \"THERMAL_SAMPLES\": \"7731\",\n", " \"ORIENTATION\": \"NORTH_UP\",\n", " \"CORNER_UL_LAT_PRODUCT\": \"41.38441\",\n", " \"CORNER_UL_LON_PRODUCT\": \"-76.26178\",\n", " \"CORNER_UR_LAT_PRODUCT\": \"41.38140\",\n", " \"CORNER_UR_LON_PRODUCT\": \"-73.48833\",\n", " \"CORNER_LL_LAT_PRODUCT\": \"39.26052\",\n", " \"CORNER_LL_LON_PRODUCT\": \"-76.22284\",\n", " \"CORNER_LR_LAT_PRODUCT\": \"39.25773\",\n", " \"CORNER_LR_LON_PRODUCT\": \"-73.53498\",\n", " \"CORNER_UL_PROJECTION_X_PRODUCT\": \"394500.000\",\n", " \"CORNER_UL_PROJECTION_Y_PRODUCT\": \"4582200.000\",\n", " \"CORNER_UR_PROJECTION_X_PRODUCT\": \"626400.000\",\n", " \"CORNER_UR_PROJECTION_Y_PRODUCT\": \"4582200.000\",\n", " \"CORNER_LL_PROJECTION_X_PRODUCT\": \"394500.000\",\n", " \"CORNER_LL_PROJECTION_Y_PRODUCT\": \"4346400.000\",\n", " \"CORNER_LR_PROJECTION_X_PRODUCT\": \"626400.000\",\n", " \"CORNER_LR_PROJECTION_Y_PRODUCT\": \"4346400.000\"\n", " },\n", " \"LEVEL2_PROCESSING_RECORD\": {\n", " \"ORIGIN\": \"Image courtesy of the U.S. Geological Survey\",\n", " \"DIGITAL_OBJECT_IDENTIFIER\": \"https://doi.org/10.5066/P9OGBGM6\",\n", " \"REQUEST_ID\": \"1626123_00008\",\n", " \"LANDSAT_PRODUCT_ID\": \"LC08_L2SP_014032_20221219_20230113_02_T1\",\n", " \"PROCESSING_LEVEL\": \"L2SP\",\n", " \"OUTPUT_FORMAT\": \"GEOTIFF\",\n", " \"DATE_PRODUCT_GENERATED\": \"2023-01-13T02:53:40Z\",\n", " \"PROCESSING_SOFTWARE_VERSION\": \"LPGS_16.1.0\",\n", " \"ALGORITHM_SOURCE_SURFACE_REFLECTANCE\": \"LaSRC_1.5.0\",\n", " \"DATA_SOURCE_OZONE\": \"MODIS\",\n", " \"DATA_SOURCE_PRESSURE\": \"Calculated\",\n", " \"DATA_SOURCE_WATER_VAPOR\": \"MODIS\",\n", " \"DATA_SOURCE_AIR_TEMPERATURE\": \"MODIS\",\n", " \"ALGORITHM_SOURCE_SURFACE_TEMPERATURE\": \"st_1.3.0\",\n", " \"DATA_SOURCE_REANALYSIS\": \"GEOS-5 FP-IT\"\n", " },\n", " \"LEVEL2_SURFACE_REFLECTANCE_PARAMETERS\": {\n", " \"REFLECTANCE_MAXIMUM_BAND_1\": \"1.602213\",\n", " \"REFLECTANCE_MINIMUM_BAND_1\": \"-0.199972\",\n", " \"REFLECTANCE_MAXIMUM_BAND_2\": \"1.602213\",\n", " \"REFLECTANCE_MINIMUM_BAND_2\": \"-0.199972\",\n", " \"REFLECTANCE_MAXIMUM_BAND_3\": \"1.602213\",\n", " \"REFLECTANCE_MINIMUM_BAND_3\": \"-0.199972\",\n", " \"REFLECTANCE_MAXIMUM_BAND_4\": \"1.602213\",\n", " \"REFLECTANCE_MINIMUM_BAND_4\": \"-0.199972\",\n", " \"REFLECTANCE_MAXIMUM_BAND_5\": \"1.602213\",\n", " \"REFLECTANCE_MINIMUM_BAND_5\": \"-0.199972\",\n", " \"REFLECTANCE_MAXIMUM_BAND_6\": \"1.602213\",\n", " \"REFLECTANCE_MINIMUM_BAND_6\": \"-0.199972\",\n", " \"REFLECTANCE_MAXIMUM_BAND_7\": \"1.602213\",\n", " \"REFLECTANCE_MINIMUM_BAND_7\": \"-0.199972\",\n", " \"QUANTIZE_CAL_MAX_BAND_1\": \"65535\",\n", " \"QUANTIZE_CAL_MIN_BAND_1\": \"1\",\n", " \"QUANTIZE_CAL_MAX_BAND_2\": \"65535\",\n", " \"QUANTIZE_CAL_MIN_BAND_2\": \"1\",\n", " \"QUANTIZE_CAL_MAX_BAND_3\": \"65535\",\n", " \"QUANTIZE_CAL_MIN_BAND_3\": \"1\",\n", " \"QUANTIZE_CAL_MAX_BAND_4\": \"65535\",\n", " \"QUANTIZE_CAL_MIN_BAND_4\": \"1\",\n", " \"QUANTIZE_CAL_MAX_BAND_5\": \"65535\",\n", " \"QUANTIZE_CAL_MIN_BAND_5\": \"1\",\n", " \"QUANTIZE_CAL_MAX_BAND_6\": \"65535\",\n", " \"QUANTIZE_CAL_MIN_BAND_6\": \"1\",\n", " \"QUANTIZE_CAL_MAX_BAND_7\": \"65535\",\n", " \"QUANTIZE_CAL_MIN_BAND_7\": \"1\",\n", " \"REFLECTANCE_MULT_BAND_1\": \"2.75e-05\",\n", " \"REFLECTANCE_MULT_BAND_2\": \"2.75e-05\",\n", " \"REFLECTANCE_MULT_BAND_3\": \"2.75e-05\",\n", " \"REFLECTANCE_MULT_BAND_4\": \"2.75e-05\",\n", " \"REFLECTANCE_MULT_BAND_5\": \"2.75e-05\",\n", " \"REFLECTANCE_MULT_BAND_6\": \"2.75e-05\",\n", " \"REFLECTANCE_MULT_BAND_7\": \"2.75e-05\",\n", " \"REFLECTANCE_ADD_BAND_1\": \"-0.2\",\n", " \"REFLECTANCE_ADD_BAND_2\": \"-0.2\",\n", " \"REFLECTANCE_ADD_BAND_3\": \"-0.2\",\n", " \"REFLECTANCE_ADD_BAND_4\": \"-0.2\",\n", " \"REFLECTANCE_ADD_BAND_5\": \"-0.2\",\n", " \"REFLECTANCE_ADD_BAND_6\": \"-0.2\",\n", " \"REFLECTANCE_ADD_BAND_7\": \"-0.2\"\n", " },\n", " \"LEVEL2_SURFACE_TEMPERATURE_PARAMETERS\": {\n", " \"TEMPERATURE_MAXIMUM_BAND_ST_B10\": \"372.999941\",\n", " \"TEMPERATURE_MINIMUM_BAND_ST_B10\": \"149.003418\",\n", " \"QUANTIZE_CAL_MAXIMUM_BAND_ST_B10\": \"65535\",\n", " \"QUANTIZE_CAL_MINIMUM_BAND_ST_B10\": \"1\",\n", " \"TEMPERATURE_MULT_BAND_ST_B10\": \"0.00341802\",\n", " \"TEMPERATURE_ADD_BAND_ST_B10\": \"149.0\"\n", " },\n", " \"LEVEL1_PROCESSING_RECORD\": {\n", " \"ORIGIN\": \"Image courtesy of the U.S. Geological Survey\",\n", " \"DIGITAL_OBJECT_IDENTIFIER\": \"https://doi.org/10.5066/P975CC9B\",\n", " \"REQUEST_ID\": \"1626123_00008\",\n", " \"LANDSAT_SCENE_ID\": \"LC80140322022353LGN00\",\n", " \"LANDSAT_PRODUCT_ID\": \"LC08_L1TP_014032_20221219_20230113_02_T1\",\n", " \"PROCESSING_LEVEL\": \"L1TP\",\n", " \"COLLECTION_CATEGORY\": \"T1\",\n", " \"OUTPUT_FORMAT\": \"GEOTIFF\",\n", " \"DATE_PRODUCT_GENERATED\": \"2023-01-13T02:38:55Z\",\n", " \"PROCESSING_SOFTWARE_VERSION\": \"LPGS_16.1.0\",\n", " \"FILE_NAME_BAND_1\": \"LC08_L1TP_014032_20221219_20230113_02_T1_B1.TIF\",\n", " \"FILE_NAME_BAND_2\": \"LC08_L1TP_014032_20221219_20230113_02_T1_B2.TIF\",\n", " \"FILE_NAME_BAND_3\": \"LC08_L1TP_014032_20221219_20230113_02_T1_B3.TIF\",\n", " \"FILE_NAME_BAND_4\": \"LC08_L1TP_014032_20221219_20230113_02_T1_B4.TIF\",\n", " \"FILE_NAME_BAND_5\": \"LC08_L1TP_014032_20221219_20230113_02_T1_B5.TIF\",\n", " \"FILE_NAME_BAND_6\": \"LC08_L1TP_014032_20221219_20230113_02_T1_B6.TIF\",\n", " \"FILE_NAME_BAND_7\": \"LC08_L1TP_014032_20221219_20230113_02_T1_B7.TIF\",\n", " \"FILE_NAME_BAND_8\": \"LC08_L1TP_014032_20221219_20230113_02_T1_B8.TIF\",\n", " \"FILE_NAME_BAND_9\": \"LC08_L1TP_014032_20221219_20230113_02_T1_B9.TIF\",\n", " \"FILE_NAME_BAND_10\": \"LC08_L1TP_014032_20221219_20230113_02_T1_B10.TIF\",\n", " \"FILE_NAME_BAND_11\": \"LC08_L1TP_014032_20221219_20230113_02_T1_B11.TIF\",\n", " \"FILE_NAME_QUALITY_L1_PIXEL\": \"LC08_L1TP_014032_20221219_20230113_02_T1_QA_PIXEL.TIF\",\n", " \"FILE_NAME_QUALITY_L1_RADIOMETRIC_SATURATION\": \"LC08_L1TP_014032_20221219_20230113_02_T1_QA_RADSAT.TIF\",\n", " \"FILE_NAME_ANGLE_COEFFICIENT\": \"LC08_L1TP_014032_20221219_20230113_02_T1_ANG.txt\",\n", " \"FILE_NAME_ANGLE_SENSOR_AZIMUTH_BAND_4\": \"LC08_L1TP_014032_20221219_20230113_02_T1_VAA.TIF\",\n", " \"FILE_NAME_ANGLE_SENSOR_ZENITH_BAND_4\": \"LC08_L1TP_014032_20221219_20230113_02_T1_VZA.TIF\",\n", " \"FILE_NAME_ANGLE_SOLAR_AZIMUTH_BAND_4\": \"LC08_L1TP_014032_20221219_20230113_02_T1_SAA.TIF\",\n", " \"FILE_NAME_ANGLE_SOLAR_ZENITH_BAND_4\": \"LC08_L1TP_014032_20221219_20230113_02_T1_SZA.TIF\",\n", " \"FILE_NAME_METADATA_ODL\": \"LC08_L1TP_014032_20221219_20230113_02_T1_MTL.txt\",\n", " \"FILE_NAME_METADATA_XML\": \"LC08_L1TP_014032_20221219_20230113_02_T1_MTL.xml\",\n", " \"FILE_NAME_CPF\": \"LC08CPF_20221001_20221231_02.03\",\n", " \"FILE_NAME_BPF_OLI\": \"LO8BPF20221219152831_20221219170353.01\",\n", " \"FILE_NAME_BPF_TIRS\": \"LT8BPF20221215135451_20221222101440.01\",\n", " \"FILE_NAME_RLUT\": \"LC08RLUT_20150303_20431231_02_01.h5\",\n", " \"DATA_SOURCE_TIRS_STRAY_LIGHT_CORRECTION\": \"TIRS\",\n", " \"DATA_SOURCE_ELEVATION\": \"GLS2000\",\n", " \"GROUND_CONTROL_POINTS_VERSION\": \"5\",\n", " \"GROUND_CONTROL_POINTS_MODEL\": \"462\",\n", " \"GEOMETRIC_RMSE_MODEL\": \"8.179\",\n", " \"GEOMETRIC_RMSE_MODEL_Y\": \"7.213\",\n", " \"GEOMETRIC_RMSE_MODEL_X\": \"3.856\",\n", " \"GROUND_CONTROL_POINTS_VERIFY\": \"120\",\n", " \"GEOMETRIC_RMSE_VERIFY\": \"7.426\"\n", " },\n", " \"LEVEL1_MIN_MAX_RADIANCE\": {\n", " \"RADIANCE_MAXIMUM_BAND_1\": \"785.06079\",\n", " \"RADIANCE_MINIMUM_BAND_1\": \"-64.83057\",\n", " \"RADIANCE_MAXIMUM_BAND_2\": \"803.91187\",\n", " \"RADIANCE_MINIMUM_BAND_2\": \"-66.38730\",\n", " \"RADIANCE_MAXIMUM_BAND_3\": \"740.79791\",\n", " \"RADIANCE_MINIMUM_BAND_3\": \"-61.17533\",\n", " \"RADIANCE_MAXIMUM_BAND_4\": \"624.68250\",\n", " \"RADIANCE_MINIMUM_BAND_4\": \"-51.58648\",\n", " \"RADIANCE_MAXIMUM_BAND_5\": \"382.27454\",\n", " \"RADIANCE_MINIMUM_BAND_5\": \"-31.56836\",\n", " \"RADIANCE_MAXIMUM_BAND_6\": \"95.06820\",\n", " \"RADIANCE_MINIMUM_BAND_6\": \"-7.85076\",\n", " \"RADIANCE_MAXIMUM_BAND_7\": \"32.04307\",\n", " \"RADIANCE_MINIMUM_BAND_7\": \"-2.64613\",\n", " \"RADIANCE_MAXIMUM_BAND_8\": \"706.96869\",\n", " \"RADIANCE_MINIMUM_BAND_8\": \"-58.38170\",\n", " \"RADIANCE_MAXIMUM_BAND_9\": \"149.40157\",\n", " \"RADIANCE_MINIMUM_BAND_9\": \"-12.33763\",\n", " \"RADIANCE_MAXIMUM_BAND_10\": \"22.00180\",\n", " \"RADIANCE_MINIMUM_BAND_10\": \"0.10033\",\n", " \"RADIANCE_MAXIMUM_BAND_11\": \"22.00180\",\n", " \"RADIANCE_MINIMUM_BAND_11\": \"0.10033\"\n", " },\n", " \"LEVEL1_MIN_MAX_REFLECTANCE\": {\n", " \"REFLECTANCE_MAXIMUM_BAND_1\": \"1.210700\",\n", " \"REFLECTANCE_MINIMUM_BAND_1\": \"-0.099980\",\n", " \"REFLECTANCE_MAXIMUM_BAND_2\": \"1.210700\",\n", " \"REFLECTANCE_MINIMUM_BAND_2\": \"-0.099980\",\n", " \"REFLECTANCE_MAXIMUM_BAND_3\": \"1.210700\",\n", " \"REFLECTANCE_MINIMUM_BAND_3\": \"-0.099980\",\n", " \"REFLECTANCE_MAXIMUM_BAND_4\": \"1.210700\",\n", " \"REFLECTANCE_MINIMUM_BAND_4\": \"-0.099980\",\n", " \"REFLECTANCE_MAXIMUM_BAND_5\": \"1.210700\",\n", " \"REFLECTANCE_MINIMUM_BAND_5\": \"-0.099980\",\n", " \"REFLECTANCE_MAXIMUM_BAND_6\": \"1.210700\",\n", " \"REFLECTANCE_MINIMUM_BAND_6\": \"-0.099980\",\n", " \"REFLECTANCE_MAXIMUM_BAND_7\": \"1.210700\",\n", " \"REFLECTANCE_MINIMUM_BAND_7\": \"-0.099980\",\n", " \"REFLECTANCE_MAXIMUM_BAND_8\": \"1.210700\",\n", " \"REFLECTANCE_MINIMUM_BAND_8\": \"-0.099980\",\n", " \"REFLECTANCE_MAXIMUM_BAND_9\": \"1.210700\",\n", " \"REFLECTANCE_MINIMUM_BAND_9\": \"-0.099980\"\n", " },\n", " \"LEVEL1_MIN_MAX_PIXEL_VALUE\": {\n", " \"QUANTIZE_CAL_MAX_BAND_1\": \"65535\",\n", " \"QUANTIZE_CAL_MIN_BAND_1\": \"1\",\n", " \"QUANTIZE_CAL_MAX_BAND_2\": \"65535\",\n", " \"QUANTIZE_CAL_MIN_BAND_2\": \"1\",\n", " \"QUANTIZE_CAL_MAX_BAND_3\": \"65535\",\n", " \"QUANTIZE_CAL_MIN_BAND_3\": \"1\",\n", " \"QUANTIZE_CAL_MAX_BAND_4\": \"65535\",\n", " \"QUANTIZE_CAL_MIN_BAND_4\": \"1\",\n", " \"QUANTIZE_CAL_MAX_BAND_5\": \"65535\",\n", " \"QUANTIZE_CAL_MIN_BAND_5\": \"1\",\n", " \"QUANTIZE_CAL_MAX_BAND_6\": \"65535\",\n", " \"QUANTIZE_CAL_MIN_BAND_6\": \"1\",\n", " \"QUANTIZE_CAL_MAX_BAND_7\": \"65535\",\n", " \"QUANTIZE_CAL_MIN_BAND_7\": \"1\",\n", " \"QUANTIZE_CAL_MAX_BAND_8\": \"65535\",\n", " \"QUANTIZE_CAL_MIN_BAND_8\": \"1\",\n", " \"QUANTIZE_CAL_MAX_BAND_9\": \"65535\",\n", " \"QUANTIZE_CAL_MIN_BAND_9\": \"1\",\n", " \"QUANTIZE_CAL_MAX_BAND_10\": \"65535\",\n", " \"QUANTIZE_CAL_MIN_BAND_10\": \"1\",\n", " \"QUANTIZE_CAL_MAX_BAND_11\": \"65535\",\n", " \"QUANTIZE_CAL_MIN_BAND_11\": \"1\"\n", " },\n", " \"LEVEL1_RADIOMETRIC_RESCALING\": {\n", " \"RADIANCE_MULT_BAND_1\": \"1.2969E-02\",\n", " \"RADIANCE_MULT_BAND_2\": \"1.3280E-02\",\n", " \"RADIANCE_MULT_BAND_3\": \"1.2238E-02\",\n", " \"RADIANCE_MULT_BAND_4\": \"1.0319E-02\",\n", " \"RADIANCE_MULT_BAND_5\": \"6.3149E-03\",\n", " \"RADIANCE_MULT_BAND_6\": \"1.5705E-03\",\n", " \"RADIANCE_MULT_BAND_7\": \"5.2933E-04\",\n", " \"RADIANCE_MULT_BAND_8\": \"1.1679E-02\",\n", " \"RADIANCE_MULT_BAND_9\": \"2.4680E-03\",\n", " \"RADIANCE_MULT_BAND_10\": \"3.3420E-04\",\n", " \"RADIANCE_MULT_BAND_11\": \"3.3420E-04\",\n", " \"RADIANCE_ADD_BAND_1\": \"-64.84355\",\n", " \"RADIANCE_ADD_BAND_2\": \"-66.40058\",\n", " \"RADIANCE_ADD_BAND_3\": \"-61.18757\",\n", " \"RADIANCE_ADD_BAND_4\": \"-51.59680\",\n", " \"RADIANCE_ADD_BAND_5\": \"-31.57467\",\n", " \"RADIANCE_ADD_BAND_6\": \"-7.85233\",\n", " \"RADIANCE_ADD_BAND_7\": \"-2.64666\",\n", " \"RADIANCE_ADD_BAND_8\": \"-58.39338\",\n", " \"RADIANCE_ADD_BAND_9\": \"-12.34010\",\n", " \"RADIANCE_ADD_BAND_10\": \"0.10000\",\n", " \"RADIANCE_ADD_BAND_11\": \"0.10000\",\n", " \"REFLECTANCE_MULT_BAND_1\": \"2.0000E-05\",\n", " \"REFLECTANCE_MULT_BAND_2\": \"2.0000E-05\",\n", " \"REFLECTANCE_MULT_BAND_3\": \"2.0000E-05\",\n", " \"REFLECTANCE_MULT_BAND_4\": \"2.0000E-05\",\n", " \"REFLECTANCE_MULT_BAND_5\": \"2.0000E-05\",\n", " \"REFLECTANCE_MULT_BAND_6\": \"2.0000E-05\",\n", " \"REFLECTANCE_MULT_BAND_7\": \"2.0000E-05\",\n", " \"REFLECTANCE_MULT_BAND_8\": \"2.0000E-05\",\n", " \"REFLECTANCE_MULT_BAND_9\": \"2.0000E-05\",\n", " \"REFLECTANCE_ADD_BAND_1\": \"-0.100000\",\n", " \"REFLECTANCE_ADD_BAND_2\": \"-0.100000\",\n", " \"REFLECTANCE_ADD_BAND_3\": \"-0.100000\",\n", " \"REFLECTANCE_ADD_BAND_4\": \"-0.100000\",\n", " \"REFLECTANCE_ADD_BAND_5\": \"-0.100000\",\n", " \"REFLECTANCE_ADD_BAND_6\": \"-0.100000\",\n", " \"REFLECTANCE_ADD_BAND_7\": \"-0.100000\",\n", " \"REFLECTANCE_ADD_BAND_8\": \"-0.100000\",\n", " \"REFLECTANCE_ADD_BAND_9\": \"-0.100000\"\n", " },\n", " \"LEVEL1_THERMAL_CONSTANTS\": {\n", " \"K1_CONSTANT_BAND_10\": \"774.8853\",\n", " \"K2_CONSTANT_BAND_10\": \"1321.0789\",\n", " \"K1_CONSTANT_BAND_11\": \"480.8883\",\n", " \"K2_CONSTANT_BAND_11\": \"1201.1442\"\n", " },\n", " \"LEVEL1_PROJECTION_PARAMETERS\": {\n", " \"MAP_PROJECTION\": \"UTM\",\n", " \"DATUM\": \"WGS84\",\n", " \"ELLIPSOID\": \"WGS84\",\n", " \"UTM_ZONE\": \"18\",\n", " \"GRID_CELL_SIZE_PANCHROMATIC\": \"15.00\",\n", " \"GRID_CELL_SIZE_REFLECTIVE\": \"30.00\",\n", " \"GRID_CELL_SIZE_THERMAL\": \"30.00\",\n", " \"ORIENTATION\": \"NORTH_UP\",\n", " \"RESAMPLING_OPTION\": \"CUBIC_CONVOLUTION\"\n", " },\n", " \"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\"\n", "}\n" ] } ], "source": [ "metadata = get_metadata(scene_mtls[0])\n", "print(json.dumps(metadata, indent=4))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "def download_sidecar(metadata: Dict[str, Any], filename: str) -> str:\n", " parsed = urlparse(metadata[\"ORIGINAL_URL\"])\n", " return urlunparse(parsed._replace(path=join(dirname(parsed.path), filename)))" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "download_url = partial(download_sidecar, metadata)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Create a STAC Item from a scene\n", "\n", "Now that we have metadata for the scene let's use it to create a [STAC Item](https://github.com/radiantearth/stac-spec/blob/master/item-spec/item-spec.md).\n", "\n", "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](https://pystac.readthedocs.io/en/latest/api.html#item)." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Help on function __init__ in module pystac.item:\n", "\n", "__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)\n", " Initialize self. See help(type(self)) for accurate signature.\n", "\n" ] } ], "source": [ "help(pystac.Item.__init__)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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.\n", "\n", "> 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." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Item `id`\n", "\n", "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: " ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "def get_item_id(metadata: Dict[str, Any]) -> str:\n", " return cast(str, metadata[\"LEVEL1_PROCESSING_RECORD\"][\"LANDSAT_SCENE_ID\"][:-5])" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'LC80140322022353'" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "item_id = get_item_id(metadata)\n", "item_id" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Item `datetime`\n", "\n", "Here we parse the datetime of the Item from two metadata fields that describe the date and time the scene was captured:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "def get_datetime(metadata: Dict[str, Any]) -> datetime:\n", " return parse(\n", " \"%sT%s\"\n", " % (\n", " metadata[\"IMAGE_ATTRIBUTES\"][\"DATE_ACQUIRED\"],\n", " metadata[\"IMAGE_ATTRIBUTES\"][\"SCENE_CENTER_TIME\"],\n", " )\n", " )" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "datetime.datetime(2022, 12, 19, 15, 40, 17, 729916, tzinfo=tzutc())" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "item_datetime = get_datetime(metadata)\n", "item_datetime" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Item `bbox`\n", "\n", "Here we read in the bounding box information from the scene and transform it into the format of the Item's `bbox` property:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "def get_bbox(metadata: Dict[str, Any]) -> List[float]:\n", " metadata = metadata[\"PROJECTION_ATTRIBUTES\"]\n", " coords = [\n", " [\n", " [\n", " float(metadata[\"CORNER_UL_LON_PRODUCT\"]),\n", " float(metadata[\"CORNER_UL_LAT_PRODUCT\"]),\n", " ],\n", " [\n", " float(metadata[\"CORNER_UR_LON_PRODUCT\"]),\n", " float(metadata[\"CORNER_UR_LAT_PRODUCT\"]),\n", " ],\n", " [\n", " float(metadata[\"CORNER_LR_LON_PRODUCT\"]),\n", " float(metadata[\"CORNER_LR_LAT_PRODUCT\"]),\n", " ],\n", " [\n", " float(metadata[\"CORNER_LL_LON_PRODUCT\"]),\n", " float(metadata[\"CORNER_LL_LAT_PRODUCT\"]),\n", " ],\n", " [\n", " float(metadata[\"CORNER_UL_LON_PRODUCT\"]),\n", " float(metadata[\"CORNER_UL_LAT_PRODUCT\"]),\n", " ],\n", " ]\n", " ]\n", " lats = [c[1] for c in coords[0]]\n", " lons = [c[0] for c in coords[0]]\n", " return [min(lons), min(lats), max(lons), max(lats)]" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[-76.26178, 39.25773, -73.48833, 41.38441]" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "item_bbox = get_bbox(metadata)\n", "item_bbox" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Item `geometry`\n", "\n", "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.\n", "\n", "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:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "def get_geometry(metadata: Dict[str, Any], bbox: List[float]) -> Dict[str, Any]:\n", " url = download_sidecar(\n", " metadata, metadata[\"PRODUCT_CONTENTS\"][\"FILE_NAME_ANGLE_COEFFICIENT\"]\n", " )\n", " sz = []\n", " coords = []\n", " ang_text = stac_io.read_text(pc.sign(url))\n", " if not ang_text.startswith(\"GROUP\"):\n", " raise ValueError(f\"ANG file for url {url} is incorrectly formatted\")\n", " for line in ang_text.split(\"\\n\"):\n", " if \"BAND01_NUM_L1T_LINES\" in line or \"BAND01_NUM_L1T_SAMPS\" in line:\n", " sz.append(float(line.split(\"=\")[1]))\n", " if (\n", " \"BAND01_L1T_IMAGE_CORNER_LINES\" in line\n", " or \"BAND01_L1T_IMAGE_CORNER_SAMPS\" in line\n", " ):\n", " coords.append(\n", " [float(v) for v in line.split(\"=\")[1].strip().strip(\"()\").split(\",\")]\n", " )\n", " if len(coords) == 2:\n", " break\n", " dlon = bbox[2] - bbox[0]\n", " dlat = bbox[3] - bbox[1]\n", " lons = [c / sz[1] * dlon + bbox[0] for c in coords[1]]\n", " lats = [((sz[0] - c) / sz[0]) * dlat + bbox[1] for c in coords[0]]\n", " coordinates = [\n", " [\n", " [lons[0], lats[0]],\n", " [lons[1], lats[1]],\n", " [lons[2], lats[2]],\n", " [lons[3], lats[3]],\n", " [lons[0], lats[0]],\n", " ]\n", " ]\n", "\n", " return {\"type\": \"Polygon\", \"coordinates\": coordinates}" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{\n", " \"type\": \"Polygon\",\n", " \"coordinates\": [\n", " [\n", " [\n", " -75.71075270108336,\n", " 41.3823086369878\n", " ],\n", " [\n", " -73.48924866988654,\n", " 40.980654308234485\n", " ],\n", " [\n", " -74.0425618957281,\n", " 39.25823722657151\n", " ],\n", " [\n", " -76.26093009667797,\n", " 39.66800780107756\n", " ],\n", " [\n", " -75.71075270108336,\n", " 41.3823086369878\n", " ]\n", " ]\n", " ]\n", "}\n" ] } ], "source": [ "item_geometry = get_geometry(metadata, item_bbox)\n", "print(json.dumps(item_geometry, indent=2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This would be a good time to check our work - we can print out the GeoJSON and use [geojson.io](https://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!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Create the item\n", "\n", "Now that we have the required attributes for an Item we can create it:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "item = pystac.Item(\n", " id=item_id,\n", " datetime=item_datetime,\n", " geometry=item_geometry,\n", " bbox=item_bbox,\n", " properties={},\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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:" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['https://schemas.stacspec.org/v1.0.0/item-spec/json-schema/item.json']" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "item.validate()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Add Ground Sample Distance to common metadata\n", "\n", "We'll add the Ground Sample Distance that is defined as part of the Item [Common Metadata](https://github.com/radiantearth/stac-spec/blob/master/item-spec/common-metadata.md). 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." ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "item.common_metadata.gsd = 30.0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Adding the EO extension\n", "\n", "STAC has a rich [set of extensions](https://stac-extensions.github.io/) 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](https://github.com/stac-extensions/eo), which captures fields needed for electro-optical data, like center wavelength and full-width half maximum values.\n", "\n", "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.\n", "\n", "For now, we will enable the EO extension for this item by using the `ext` property provided by the extension object:" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "eo_ext = EOExtension.ext(item, add_if_missing=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Add cloud cover\n", "\n", "Here we add cloud cover from the metadata as part of the `eo` extension." ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "def get_cloud_cover(metadata: Dict[str, Any]) -> float:\n", " return float(metadata[\"IMAGE_ATTRIBUTES\"][\"CLOUD_COVER\"])" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "43.42" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "eo_ext.cloud_cover = get_cloud_cover(metadata)\n", "eo_ext.cloud_cover" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Adding assets\n", "\n", "STAC Items contain a list of [Assets](https://github.com/radiantearth/stac-spec/blob/master/item-spec/item-spec.md#asset-object), 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.\n", "\n", "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:" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "class BandInfo(TypedDict):\n", " name: str\n", " asset_fields: Dict[str, str]\n", " extensions: List[Union[EOBand, RasterBand, List[Bitfield]]]\n", "\n", "\n", "def eo_band_info(\n", " common_name: str,\n", " href: str,\n", " name: str,\n", " description: str,\n", " center: float,\n", " fwhm: float,\n", " default_raster_band: Optional[RasterBand] = None,\n", "):\n", " raster_band = (\n", " RasterBand.create(\n", " spatial_resolution=30.0,\n", " scale=0.0000275,\n", " nodata=0,\n", " offset=-0.2,\n", " data_type=\"uint16\",\n", " )\n", " if default_raster_band is None\n", " else default_raster_band\n", " )\n", " return {\n", " \"name\": common_name,\n", " \"asset_fields\": {\n", " \"href\": href,\n", " \"media_type\": str(pystac.media_type.MediaType.COG),\n", " },\n", " \"extensions\": [\n", " EOBand.create(\n", " name=name,\n", " common_name=common_name,\n", " description=description,\n", " center_wavelength=center,\n", " full_width_half_max=fwhm,\n", " ),\n", " raster_band,\n", " ],\n", " }\n", "\n", "\n", "def plain_band_info(name: str, href: str, title: str, ext: RasterBand):\n", " return {\n", " \"name\": name,\n", " \"asset_fields\": {\n", " \"href\": href,\n", " \"media_type\": str(pystac.media_type.MediaType.COG),\n", " \"title\": title,\n", " },\n", " \"extensions\": [ext],\n", " }" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Some common raster band information definitions will also be useful." ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [], "source": [ "thermal_raster_band = RasterBand.create(\n", " spatial_resolution=30.0,\n", " scale=0.00341802,\n", " nodata=0,\n", " offset=149.0,\n", " data_type=\"uint6\",\n", " unit=\"kelvin\",\n", ")\n", "radiance_raster_band = RasterBand.create(\n", " unit=\"watt/steradian/square_meter/micrometer\",\n", " scale=1e-3,\n", " nodata=-9999,\n", " data_type=\"uint16\",\n", " spatial_resolution=30.0,\n", ")\n", "emissivity_transmission_raster_band = RasterBand.create(\n", " scale=1e-4,\n", " nodata=-9999,\n", " data_type=\"int16\",\n", " spatial_resolution=30.0,\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Several QA bands are provided that utilize bit-wise masks which we can define using the [`classification` extension](https://github.com/stac-extensions/classification). Because these definitions can be verbose, we provide some additional helper functions to minimize the length of their definition." ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [], "source": [ "def create_bitfield(\n", " offset: int,\n", " length: int,\n", " name: str,\n", " field_names_descriptions: List[Tuple[str, str]],\n", " description: Optional[str] = None,\n", ") -> Bitfield:\n", " return Bitfield.create(\n", " offset=offset,\n", " length=length,\n", " name=name,\n", " description=description,\n", " classes=[\n", " Classification.create(value=i, name=n, description=d)\n", " for (i, (n, d)) in enumerate(field_names_descriptions)\n", " ],\n", " )\n", "\n", "\n", "def create_qa_bitfield(\n", " offset: int,\n", " class_name: str,\n", " description: Optional[Union[str, Tuple[str, str]]] = None,\n", ") -> Bitfield:\n", " if description is None:\n", " descr0 = f\"{class_name.replace('_', ' ').capitalize()} confidence is not high\"\n", " descr1 = f\"High confidence {class_name.replace('_', ' ')}\"\n", " elif isinstance(description, str):\n", " descr0 = f\"{description.capitalize()} confidence is not high\"\n", " descr1 = f\"High confidence {description.lower()}\"\n", " else:\n", " descr0 = description[0]\n", " descr1 = description[1]\n", "\n", " return create_bitfield(\n", " offset, 1, class_name, [(f\"not_{class_name}\", descr0), (class_name, descr1)]\n", " )\n", "\n", "\n", "def create_confidence_bitfield(\n", " offset: int, class_name: str, use_medium: bool = False\n", ") -> Bitfield:\n", " label = class_name.replace(\"_\", \" \")\n", " return create_bitfield(\n", " offset,\n", " 2,\n", " f\"{class_name}_confidence\",\n", " [\n", " (\"not_set\", \"No confidence level set\"),\n", " (\"low\", f\"Low confidence {label}\"),\n", " (\n", " (\"medium\", f\"Medium confidence {label}\")\n", " if use_medium\n", " else (\"reserved\", \"Reserved - value not used\")\n", " ),\n", " (\"high\", f\"High confidence {label}\"),\n", " ],\n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [], "source": [ "def landsat_band_info(\n", " metadata: Dict[str, Any], downloader: Callable[[str], str]\n", ") -> List[BandInfo]:\n", " product_contents = metadata[\"PRODUCT_CONTENTS\"]\n", " return [\n", " {\n", " \"name\": \"ang\",\n", " \"asset_fields\": {\n", " \"href\": downloader(product_contents[\"FILE_NAME_ANGLE_COEFFICIENT\"]),\n", " \"media_type\": \"text/plain\",\n", " \"title\": \"Angle coefficients\",\n", " },\n", " \"extensions\": [],\n", " },\n", " {\n", " \"name\": \"mtl.txt\",\n", " \"asset_fields\": {\n", " \"href\": downloader(product_contents[\"FILE_NAME_METADATA_ODL\"]),\n", " \"media_type\": \"text/plain\",\n", " \"title\": \"Product metadata\",\n", " },\n", " \"extensions\": [],\n", " },\n", " eo_band_info(\n", " \"coastal\",\n", " downloader(product_contents[\"FILE_NAME_BAND_1\"]),\n", " \"OLI_B1\",\n", " \"Coastal/Aerosol (Operational Land Imager)\",\n", " 0.44,\n", " 0.02,\n", " ),\n", " eo_band_info(\n", " \"blue\",\n", " downloader(product_contents[\"FILE_NAME_BAND_2\"]),\n", " \"OLI_B2\",\n", " \"Visible blue (Operational Land Imager)\",\n", " 0.48,\n", " 0.06,\n", " ),\n", " eo_band_info(\n", " \"green\",\n", " downloader(product_contents[\"FILE_NAME_BAND_3\"]),\n", " \"OLI_B3\",\n", " \"Visible green (Operational Land Imager)\",\n", " 0.56,\n", " 0.06,\n", " ),\n", " eo_band_info(\n", " \"red\",\n", " downloader(product_contents[\"FILE_NAME_BAND_4\"]),\n", " \"OLI_B4\",\n", " \"Visible red (Operational Land Imager)\",\n", " 0.65,\n", " 0.04,\n", " ),\n", " eo_band_info(\n", " \"nir08\",\n", " downloader(product_contents[\"FILE_NAME_BAND_5\"]),\n", " \"OLI_B5\",\n", " \"Near infrared (Operational Land Imager)\",\n", " 0.87,\n", " 0.03,\n", " ),\n", " eo_band_info(\n", " \"swir16\",\n", " downloader(product_contents[\"FILE_NAME_BAND_6\"]),\n", " \"OLI_B6\",\n", " \"Short-wave infrared (Operational Land Imager)\",\n", " 1.61,\n", " 0.09,\n", " ),\n", " eo_band_info(\n", " \"swir22\",\n", " downloader(product_contents[\"FILE_NAME_BAND_7\"]),\n", " \"OLI_B7\",\n", " \"Short-wave infrared (Operational Land Imager)\",\n", " 2.2,\n", " 0.19,\n", " ),\n", " eo_band_info(\n", " \"lwir11\",\n", " downloader(product_contents[\"FILE_NAME_BAND_ST_B10\"]),\n", " \"TIRS_B10\",\n", " \"Long-wave infrared (Thermal InfraRed Sensor)\",\n", " 10.9,\n", " 0.59,\n", " thermal_raster_band,\n", " ),\n", " plain_band_info(\n", " \"trad\",\n", " downloader(product_contents[\"FILE_NAME_THERMAL_RADIANCE\"]),\n", " \"Thermal radiance\",\n", " radiance_raster_band,\n", " ),\n", " plain_band_info(\n", " \"urad\",\n", " downloader(product_contents[\"FILE_NAME_UPWELL_RADIANCE\"]),\n", " \"Upwelled radiance\",\n", " radiance_raster_band,\n", " ),\n", " plain_band_info(\n", " \"drad\",\n", " downloader(product_contents[\"FILE_NAME_DOWNWELL_RADIANCE\"]),\n", " \"Downwelled radiance\",\n", " radiance_raster_band,\n", " ),\n", " plain_band_info(\n", " \"emis\",\n", " downloader(product_contents[\"FILE_NAME_EMISSIVITY\"]),\n", " \"Emissivity\",\n", " emissivity_transmission_raster_band,\n", " ),\n", " plain_band_info(\n", " \"emsd\",\n", " downloader(product_contents[\"FILE_NAME_EMISSIVITY_STDEV\"]),\n", " \"Emissivity standard deviation\",\n", " emissivity_transmission_raster_band,\n", " ),\n", " plain_band_info(\n", " \"atran\",\n", " downloader(product_contents[\"FILE_NAME_ATMOSPHERIC_TRANSMITTANCE\"]),\n", " \"Atmospheric transmission\",\n", " emissivity_transmission_raster_band,\n", " ),\n", " plain_band_info(\n", " \"cdist\",\n", " downloader(product_contents[\"FILE_NAME_CLOUD_DISTANCE\"]),\n", " \"Cloud distance\",\n", " RasterBand.create(\n", " unit=\"kilometer\",\n", " scale=1e-2,\n", " nodata=-9999,\n", " data_type=\"uint16\",\n", " spatial_resolution=30.0,\n", " ),\n", " ),\n", " {\n", " \"name\": \"qa\",\n", " \"asset_fields\": {\n", " \"href\": downloader(\n", " product_contents[\"FILE_NAME_QUALITY_L2_SURFACE_TEMPERATURE\"]\n", " ),\n", " \"title\": \"Surface Temperature Quality Assessment Band\",\n", " },\n", " \"extensions\": [\n", " RasterBand.create(\n", " unit=\"kelvin\",\n", " scale=1e-2,\n", " nodata=-9999,\n", " data_type=\"int16\",\n", " spatial_resolution=30,\n", " )\n", " ],\n", " },\n", " {\n", " \"name\": \"qa_pixel\",\n", " \"asset_fields\": {\n", " \"href\": downloader(product_contents[\"FILE_NAME_QUALITY_L1_PIXEL\"]),\n", " \"media_type\": str(pystac.media_type.MediaType.COG),\n", " \"title\": \"Pixel quality assessment\",\n", " },\n", " \"extensions\": [\n", " [\n", " create_qa_bitfield(0, \"fill\", (\"Image data\", \"Fill data\")),\n", " create_qa_bitfield(\n", " 1,\n", " \"dilated_cloud\",\n", " (\"Cloud is not dilated or no cloud\", \"Dilated cloud\"),\n", " ),\n", " create_qa_bitfield(2, \"cirrus\"),\n", " create_qa_bitfield(3, \"cloud\"),\n", " create_qa_bitfield(4, \"cloud_shadow\"),\n", " create_qa_bitfield(5, \"snow\"),\n", " create_qa_bitfield(6, \"clear\"),\n", " create_qa_bitfield(7, \"water\"),\n", " create_confidence_bitfield(8, \"cloud\", True),\n", " create_confidence_bitfield(10, \"cloud_shadow\"),\n", " create_confidence_bitfield(12, \"snow\"),\n", " create_confidence_bitfield(14, \"cirrus\"),\n", " ]\n", " ],\n", " },\n", " {\n", " \"name\": \"qa_radsat\",\n", " \"asset_fields\": {\n", " \"href\": downloader(\n", " product_contents[\"FILE_NAME_QUALITY_L1_RADIOMETRIC_SATURATION\"]\n", " ),\n", " \"media_type\": str(pystac.media_type.MediaType.COG),\n", " \"description\": (\n", " \"Collection 2 Level-1 Radiometric Saturation and \"\n", " \"Terrain Occlusion Quality Assessment \"\n", " \"Band (QA_RADSAT)\"\n", " ),\n", " },\n", " \"extensions\": [\n", " [\n", " Bitfield.create(\n", " offset=i - 1,\n", " length=1,\n", " description=f\"Band {i} radiometric saturation\",\n", " classes=[\n", " Classification.create(\n", " 0, f\"Band {i} not saturated\", \"not_saturated\"\n", " ),\n", " Classification.create(\n", " 1, f\"Band {i} saturated\", \"saturated\"\n", " ),\n", " ],\n", " )\n", " for i in [1, 2, 3, 4, 5, 6, 7, 9]\n", " ]\n", " + [\n", " Bitfield.create(\n", " offset=11,\n", " length=1,\n", " description=(\n", " \"Terrain not visible from sensor due to \"\n", " \"intervening terrain\"\n", " ),\n", " classes=[\n", " Classification.create(\n", " 0, \"Terrain is not occluded\", \"not_occluded\"\n", " ),\n", " Classification.create(1, \"Terrain is occluded\", \"occluded\"),\n", " ],\n", " )\n", " ]\n", " ],\n", " },\n", " {\n", " \"name\": \"qa_aerosol\",\n", " \"asset_fields\": {\n", " \"href\": downloader(product_contents[\"FILE_NAME_QUALITY_L2_AEROSOL\"]),\n", " \"media_type\": str(pystac.media_type.MediaType.COG),\n", " \"title\": \"Aerosol Quality Assessment Band\",\n", " },\n", " \"extensions\": [\n", " [\n", " create_bitfield(\n", " 0,\n", " 1,\n", " \"fill\",\n", " [(\"not_fill\", \"Pixel is not fill\"), (\"fill\", \"Pixel is fill\")],\n", " \"Image or fill data\",\n", " ),\n", " create_bitfield(\n", " 1,\n", " 1,\n", " \"retrieval\",\n", " [\n", " (\"not_valid\", \"Pixel retrieval is not valid\"),\n", " (\"valid\", \"Pixel retrieval is valid\"),\n", " ],\n", " \"Valid aerosol retrieval\",\n", " ),\n", " create_bitfield(\n", " 2,\n", " 1,\n", " \"water\",\n", " [\n", " (\"not_water\", \"Pixel is not water\"),\n", " (\"water\", \"Pixel is water\"),\n", " ],\n", " \"Water mask\",\n", " ),\n", " create_bitfield(\n", " 5,\n", " 1,\n", " \"interpolated\",\n", " [\n", " (\"not_interpolated\", \"Pixel is not interpolated\"),\n", " (\"interpolated\", \"Pixel is interpolated\"),\n", " ],\n", " \"Aerosol interpolation\",\n", " ),\n", " create_bitfield(\n", " 6,\n", " 2,\n", " \"level\",\n", " [\n", " (\"climatology\", \"No aerosol correction applied\"),\n", " (\"low\", \"Low aerosol level\"),\n", " (\"medium\", \"Medium aerosol level\"),\n", " (\"high\", \"High aerosol level\"),\n", " ],\n", " \"Aerosol level\",\n", " ),\n", " ]\n", " ],\n", " },\n", " ]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For illustration purposes, we can look at the band info records for an example scene:" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[{'name': 'ang',\n", " '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',\n", " 'media_type': 'text/plain',\n", " 'title': 'Angle coefficients'},\n", " 'extensions': []},\n", " {'name': 'mtl.txt',\n", " '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',\n", " 'media_type': 'text/plain',\n", " 'title': 'Product metadata'},\n", " 'extensions': []},\n", " {'name': 'coastal',\n", " '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',\n", " 'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized'},\n", " 'extensions': []},\n", " {'name': 'blue',\n", " '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',\n", " 'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized'},\n", " 'extensions': []},\n", " {'name': 'green',\n", " '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',\n", " 'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized'},\n", " 'extensions': []},\n", " {'name': 'red',\n", " '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',\n", " 'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized'},\n", " 'extensions': []},\n", " {'name': 'nir08',\n", " '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',\n", " 'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized'},\n", " 'extensions': []},\n", " {'name': 'swir16',\n", " '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',\n", " 'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized'},\n", " 'extensions': []},\n", " {'name': 'swir22',\n", " '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',\n", " 'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized'},\n", " 'extensions': []},\n", " {'name': 'lwir11',\n", " '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',\n", " 'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized'},\n", " 'extensions': []},\n", " {'name': 'trad',\n", " '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',\n", " 'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized',\n", " 'title': 'Thermal radiance'},\n", " 'extensions': []},\n", " {'name': 'urad',\n", " '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',\n", " 'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized',\n", " 'title': 'Upwelled radiance'},\n", " 'extensions': []},\n", " {'name': 'drad',\n", " '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',\n", " 'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized',\n", " 'title': 'Downwelled radiance'},\n", " 'extensions': []},\n", " {'name': 'emis',\n", " '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',\n", " 'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized',\n", " 'title': 'Emissivity'},\n", " 'extensions': []},\n", " {'name': 'emsd',\n", " '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',\n", " 'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized',\n", " 'title': 'Emissivity standard deviation'},\n", " 'extensions': []},\n", " {'name': 'atran',\n", " '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',\n", " 'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized',\n", " 'title': 'Atmospheric transmission'},\n", " 'extensions': []},\n", " {'name': 'cdist',\n", " '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',\n", " 'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized',\n", " 'title': 'Cloud distance'},\n", " 'extensions': []},\n", " {'name': 'qa',\n", " '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',\n", " 'title': 'Surface Temperature Quality Assessment Band'},\n", " 'extensions': []},\n", " {'name': 'qa_pixel',\n", " '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',\n", " 'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized',\n", " 'title': 'Pixel quality assessment'},\n", " 'extensions': [[, ]>,\n", " , ]>,\n", " , ]>,\n", " , ]>,\n", " , ]>,\n", " , ]>,\n", " , ]>,\n", " , ]>,\n", " , , , ]>,\n", " , , , ]>,\n", " , , , ]>,\n", " , , , ]>]]},\n", " {'name': 'qa_radsat',\n", " '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',\n", " 'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized',\n", " 'description': 'Collection 2 Level-1 Radiometric Saturation and Terrain Occlusion Quality Assessment Band (QA_RADSAT)'},\n", " 'extensions': [[, ]>,\n", " , ]>,\n", " , ]>,\n", " , ]>,\n", " , ]>,\n", " , ]>,\n", " , ]>,\n", " , ]>,\n", " , ]>]]},\n", " {'name': 'qa_aerosol',\n", " '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',\n", " 'media_type': 'image/tiff; application=geotiff; profile=cloud-optimized',\n", " 'title': 'Aerosol Quality Assessment Band'},\n", " 'extensions': [[, ]>,\n", " , ]>,\n", " , ]>,\n", " , ]>,\n", " , , , ]>]]}]" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "landsat_band_info(metadata, download_url)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With this information we can now define a method that adds all the relevant assets for a scene to an item:" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [], "source": [ "def add_assets(item: pystac.Item, band_info: List[BandInfo]) -> None:\n", " for band in band_info:\n", " asset = pystac.Asset(**band[\"asset_fields\"])\n", " asset.set_owner(item)\n", " for ext_data in band[\"extensions\"]:\n", " if isinstance(ext_data, EOBand):\n", " EOExtension.ext(asset, add_if_missing=True).bands = [ext_data]\n", " elif isinstance(ext_data, RasterBand):\n", " RasterExtension.ext(asset, add_if_missing=True).bands = [ext_data]\n", " elif isinstance(ext_data, list):\n", " ClassificationExtension.ext(\n", " asset, add_if_missing=True\n", " ).bitfields = ext_data\n", " item.add_asset(band[\"name\"], asset)" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [], "source": [ "add_assets(item, landsat_band_info(metadata, download_url))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can examine the item to ensure that the assets appear as expected." ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'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',\n", " 'type': 'image/tiff; application=geotiff; profile=cloud-optimized',\n", " 'eo:bands': [{'name': 'OLI_B5',\n", " 'common_name': 'nir08',\n", " 'description': 'Near infrared (Operational Land Imager)',\n", " 'center_wavelength': 0.87,\n", " 'full_width_half_max': 0.03}]}" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "item.assets[\"nir08\"].to_dict()" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'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',\n", " 'type': 'image/tiff; application=geotiff; profile=cloud-optimized',\n", " 'title': 'Aerosol Quality Assessment Band',\n", " 'classification:bitfields': [{'offset': 0,\n", " 'length': 1,\n", " 'classes': [{'value': 0,\n", " 'name': 'not_fill',\n", " 'description': 'Pixel is not fill'},\n", " {'value': 1, 'name': 'fill', 'description': 'Pixel is fill'}],\n", " 'description': 'Image or fill data',\n", " 'name': 'fill'},\n", " {'offset': 1,\n", " 'length': 1,\n", " 'classes': [{'value': 0,\n", " 'name': 'not_valid',\n", " 'description': 'Pixel retrieval is not valid'},\n", " {'value': 1, 'name': 'valid', 'description': 'Pixel retrieval is valid'}],\n", " 'description': 'Valid aerosol retrieval',\n", " 'name': 'retrieval'},\n", " {'offset': 2,\n", " 'length': 1,\n", " 'classes': [{'value': 0,\n", " 'name': 'not_water',\n", " 'description': 'Pixel is not water'},\n", " {'value': 1, 'name': 'water', 'description': 'Pixel is water'}],\n", " 'description': 'Water mask',\n", " 'name': 'water'},\n", " {'offset': 5,\n", " 'length': 1,\n", " 'classes': [{'value': 0,\n", " 'name': 'not_interpolated',\n", " 'description': 'Pixel is not interpolated'},\n", " {'value': 1,\n", " 'name': 'interpolated',\n", " 'description': 'Pixel is interpolated'}],\n", " 'description': 'Aerosol interpolation',\n", " 'name': 'interpolated'},\n", " {'offset': 6,\n", " 'length': 2,\n", " 'classes': [{'value': 0,\n", " 'name': 'climatology',\n", " 'description': 'No aerosol correction applied'},\n", " {'value': 1, 'name': 'low', 'description': 'Low aerosol level'},\n", " {'value': 2, 'name': 'medium', 'description': 'Medium aerosol level'},\n", " {'value': 3, 'name': 'high', 'description': 'High aerosol level'}],\n", " 'description': 'Aerosol level',\n", " 'name': 'level'}]}" ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" } ], "source": [ "item.assets[\"qa_aerosol\"].to_dict()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Add projection information\n", "\n", "We can specify the EPSG code for the scene as part of the [projection extension](https://github.com/stac-extensions/projection). The below method, adapted from [stactools](https://github.com/stactools-packages/landsat/blob/9f595a9d5ed6b62a2e96338e79f5bb502a7d90d0/src/stactools/landsat/mtl_metadata.py#L86-L109), figures out the correct UTM Zone EPSG:" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [], "source": [ "def get_epsg(metadata: Dict[str, Any], min_lat: float, max_lat: float) -> Optional[int]:\n", " if \"UTM_ZONE\" in metadata[\"PROJECTION_ATTRIBUTES\"]:\n", " utm_zone_integer = metadata[\"PROJECTION_ATTRIBUTES\"][\"UTM_ZONE\"].zfill(2)\n", " return int(f\"326{utm_zone_integer}\")\n", " else:\n", " lat_ts = metadata[\"PROJECTION_ATTRIBUTES\"][\"TRUE_SCALE_LAT\"]\n", " if lat_ts == \"-71.00000\":\n", " # Antarctic\n", " return 3031\n", " elif lat_ts == \"71.00000\":\n", " # Arctic\n", " return 3995\n", " else:\n", " raise ValueError(\n", " f\"Unexpeced value for PROJECTION_ATTRIBUTES/TRUE_SCALE_LAT: {lat_ts} \"\n", " )" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "32618" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "proj_ext = ProjectionExtension.ext(item, add_if_missing=True)\n", "assert item.bbox is not None\n", "proj_ext.epsg = get_epsg(metadata, item.bbox[1], item.bbox[3])\n", "proj_ext.epsg" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Add view geometry information\n", "\n", "The [View Geometry](https://github.com/stac-extensions/view) 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:" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [], "source": [ "def get_view_info(metadata: Dict[str, Any]) -> Dict[str, float]:\n", " return {\n", " \"sun_azimuth\": float(metadata[\"IMAGE_ATTRIBUTES\"][\"SUN_AZIMUTH\"]),\n", " \"sun_elevation\": float(metadata[\"IMAGE_ATTRIBUTES\"][\"SUN_ELEVATION\"]),\n", " }" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'datetime': '2022-12-19T15:40:17.729916Z',\n", " 'gsd': 30.0,\n", " 'eo:cloud_cover': 43.42,\n", " 'proj:epsg': 32618,\n", " 'view:sun_azimuth': 160.86021018,\n", " 'view:sun_elevation': 23.81656674}" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "view_ext = ViewExtension.ext(item, add_if_missing=True)\n", "view_info = get_view_info(metadata)\n", "view_ext.sun_azimuth = view_info[\"sun_azimuth\"]\n", "view_ext.sun_elevation = view_info[\"sun_elevation\"]\n", "item.properties" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['https://schemas.stacspec.org/v1.0.0/item-spec/json-schema/item.json',\n", " 'https://stac-extensions.github.io/eo/v1.1.0/schema.json',\n", " 'https://stac-extensions.github.io/raster/v1.1.0/schema.json',\n", " 'https://stac-extensions.github.io/classification/v1.1.0/schema.json',\n", " 'https://stac-extensions.github.io/projection/v1.1.0/schema.json',\n", " 'https://stac-extensions.github.io/view/v1.0.0/schema.json']" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "item.validate()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Building the Collection\n", "\n", "Now that we know how to build an Item for a scene, let's build the Collection that will contain all the Items.\n", "\n", "If we look at the `__init__` method for `pystac.Collection`, we can see what properties are required:" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Help on function __init__ in module pystac.collection:\n", "\n", "__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)\n", " Initialize self. See help(type(self)) for accurate signature.\n", "\n" ] } ], "source": [ "help(pystac.Collection.__init__)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Collection `id`\n", "\n", "We'll use the location name we defined above in the ID for our Collection:" ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'philly-landsat-collection-2'" ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" } ], "source": [ "collection_id = \"{}-landsat-collection-2\".format(location_name.lower())\n", "collection_id" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Collection `title`\n", "\n", "Here we set a simple title for our collection." ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'2022 Landsat images over philly'" ] }, "execution_count": 40, "metadata": {}, "output_type": "execute_result" } ], "source": [ "collection_title = \"2022 Landsat images over {}\".format(location_name.lower())\n", "collection_title" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Collection `description`\n", "\n", "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!\n", "\n", "Notice we are using [Markdown](https://www.markdownguide.org/) 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." ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "### Philly Landsat Collection 2\n", "\n", "A collection of Landsat scenes around Philly in 2022.\n", "\n" ] } ], "source": [ "collection_description = \"\"\"### {} Landsat Collection 2\n", "\n", "A collection of Landsat scenes around {} in 2022.\n", "\"\"\".format(location_name, location_name)\n", "print(collection_description)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Collection `extent`\n", "\n", "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.\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [], "source": [ "spatial_extent = pystac.SpatialExtent([[-180.0, -90.0, 180.0, 90.0]])\n", "temporal_extent = pystac.TemporalExtent([[datetime(2013, 6, 1), None]])\n", "collection_extent = pystac.Extent(spatial_extent, temporal_extent)" ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [], "source": [ "collection = pystac.Collection(\n", " id=collection_id,\n", " title=collection_title,\n", " description=collection_description,\n", " extent=collection_extent,\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now look at our Collection as a `dict` to check our values." ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'type': 'Collection',\n", " 'id': 'philly-landsat-collection-2',\n", " 'stac_version': '1.0.0',\n", " 'description': '### Philly Landsat Collection 2\\n\\nA collection of Landsat scenes around Philly in 2022.\\n',\n", " 'links': [],\n", " 'title': '2022 Landsat images over philly',\n", " 'extent': {'spatial': {'bbox': [[-180.0, -90.0, 180.0, 90.0]]},\n", " 'temporal': {'interval': [['2013-06-01T00:00:00Z', None]]}},\n", " 'license': 'proprietary'}" ] }, "execution_count": 44, "metadata": {}, "output_type": "execute_result" } ], "source": [ "collection.to_dict()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Set the license\n", "\n", "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](https://spdx.org/licenses/) string for public domain data:" ] }, { "cell_type": "code", "execution_count": 45, "metadata": {}, "outputs": [], "source": [ "collection_license = \"PDDL-1.0\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Set the providers\n", "\n", "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:" ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [], "source": [ "collection.providers = [\n", " pystac.Provider(\n", " name=\"NASA\",\n", " roles=[pystac.ProviderRole.PRODUCER, pystac.ProviderRole.LICENSOR],\n", " url=\"https://landsat.gsfc.nasa.gov/\",\n", " ),\n", " pystac.Provider(\n", " name=\"USGS\",\n", " roles=[\n", " pystac.ProviderRole.PROCESSOR,\n", " pystac.ProviderRole.PRODUCER,\n", " pystac.ProviderRole.LICENSOR,\n", " ],\n", " url=\"https://www.usgs.gov/landsat-missions/landsat-collection-2-level-2-science-products\",\n", " ),\n", " pystac.Provider(\n", " name=\"Microsoft\",\n", " roles=[pystac.ProviderRole.HOST],\n", " url=\"https://planetarycomputer.microsoft.com\",\n", " ),\n", "]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Create items for each scene\n", "\n", "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:" ] }, { "cell_type": "code", "execution_count": 47, "metadata": {}, "outputs": [], "source": [ "def item_from_metadata(mtl_xml_url: str) -> pystac.Item:\n", " metadata = get_metadata(mtl_xml_url)\n", " download_url = partial(download_sidecar, metadata)\n", "\n", " bbox = get_bbox(metadata)\n", " item = pystac.Item(\n", " id=get_item_id(metadata),\n", " datetime=get_datetime(metadata),\n", " geometry=get_geometry(metadata, bbox),\n", " bbox=bbox,\n", " properties={},\n", " )\n", "\n", " item.common_metadata.gsd = 30.0\n", "\n", " item_eo_ext = EOExtension.ext(item, add_if_missing=True)\n", " item_eo_ext.cloud_cover = get_cloud_cover(metadata)\n", "\n", " add_assets(item, landsat_band_info(metadata, download_url))\n", "\n", " item_proj_ext = ProjectionExtension.ext(item, add_if_missing=True)\n", " assert item.bbox is not None\n", " item_proj_ext.epsg = get_epsg(metadata, item.bbox[1], item.bbox[3])\n", "\n", " item_view_ext = ViewExtension.ext(item, add_if_missing=True)\n", " view_info = get_view_info(metadata)\n", " item_view_ext.sun_azimuth = view_info[\"sun_azimuth\"]\n", " item_view_ext.sun_elevation = view_info[\"sun_elevation\"]\n", "\n", " item.validate()\n", "\n", " return item" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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:" ] }, { "cell_type": "code", "execution_count": 48, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "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\n" ] } ], "source": [ "for url in scene_mtls:\n", " try:\n", " item = item_from_metadata(url)\n", " collection.add_item(item)\n", " except Exception as e:\n", " print(e)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Reset collection extent based on items\n", "\n", "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:" ] }, { "cell_type": "code", "execution_count": 49, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'spatial': {'bbox': [[-76.29048, 39.25502, -73.48833, 41.38441]]},\n", " 'temporal': {'interval': [['2022-10-08T15:40:14.577173Z',\n", " '2022-12-19T15:40:17.729916Z']]}}" ] }, "execution_count": 49, "metadata": {}, "output_type": "execute_result" } ], "source": [ "collection.update_extent_from_items()\n", "collection.extent.to_dict()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Set the HREFs for everything in the catalog\n", "\n", "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! \n", "\n", "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](https://github.com/radiantearth/stac-spec/blob/master/best-practices.md#catalog-layout) recommendations for STAC file layout for each Catalog, Collection and Item in the STAC.\n", "\n", "Here we use that method and set the root directory to a subdirectory of our user's `home` directory:" ] }, { "cell_type": "code", "execution_count": 50, "metadata": {}, "outputs": [], "source": [ "from pathlib import Path\n", "\n", "root_path = str(Path.home() / \"{}-landsat-stac\".format(location_name))\n", "\n", "collection.normalize_hrefs(root_path)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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:" ] }, { "cell_type": "code", "execution_count": 51, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "8" ] }, "execution_count": 51, "metadata": {}, "output_type": "execute_result" } ], "source": [ "collection.validate_all()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Write the catalog locally\n", "\n", "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](https://github.com/radiantearth/stac-spec/blob/master/best-practices.md#use-of-links) documentation.\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": 52, "metadata": {}, "outputs": [], "source": [ "collection.save(pystac.CatalogType.SELF_CONTAINED)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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:" ] }, { "cell_type": "code", "execution_count": 53, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "* \n", " * \n", " * \n", " * \n", " * \n", " * \n", " * \n", " * \n", " * \n" ] } ], "source": [ "collection.describe()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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:" ] }, { "cell_type": "code", "execution_count": 54, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'type': 'Collection',\n", " 'id': 'philly-landsat-collection-2',\n", " 'stac_version': '1.0.0',\n", " 'description': '### Philly Landsat Collection 2\\n\\nA collection of Landsat scenes around Philly in 2022.\\n',\n", " 'links': [{'rel': 'root',\n", " 'href': './collection.json',\n", " 'type': 'application/json',\n", " 'title': '2022 Landsat images over philly'},\n", " {'rel': 'item',\n", " 'href': './LC80140322022353/LC80140322022353.json',\n", " 'type': 'application/json'},\n", " {'rel': 'item',\n", " 'href': './LC90140322022345/LC90140322022345.json',\n", " 'type': 'application/json'},\n", " {'rel': 'item',\n", " 'href': './LC80140322022337/LC80140322022337.json',\n", " 'type': 'application/json'},\n", " {'rel': 'item',\n", " 'href': './LC90140322022329/LC90140322022329.json',\n", " 'type': 'application/json'},\n", " {'rel': 'item',\n", " 'href': './LC80140322022321/LC80140322022321.json',\n", " 'type': 'application/json'},\n", " {'rel': 'item',\n", " 'href': './LC90140322022313/LC90140322022313.json',\n", " 'type': 'application/json'},\n", " {'rel': 'item',\n", " 'href': './LC80140322022305/LC80140322022305.json',\n", " 'type': 'application/json'},\n", " {'rel': 'item',\n", " 'href': './LC90140322022281/LC90140322022281.json',\n", " 'type': 'application/json'},\n", " {'rel': 'self',\n", " 'href': '/Users/pjh/Philly-landsat-stac/collection.json',\n", " 'type': 'application/json'}],\n", " 'title': '2022 Landsat images over philly',\n", " 'extent': {'spatial': {'bbox': [[-76.29048, 39.25502, -73.48833, 41.38441]]},\n", " 'temporal': {'interval': [['2022-10-08T15:40:14.577173Z',\n", " '2022-12-19T15:40:17.729916Z']]}},\n", " 'license': 'proprietary',\n", " 'providers': [{'name': 'NASA',\n", " 'roles': [,\n", " ],\n", " 'url': 'https://landsat.gsfc.nasa.gov/'},\n", " {'name': 'USGS',\n", " 'roles': [,\n", " ,\n", " ],\n", " 'url': 'https://www.usgs.gov/landsat-missions/landsat-collection-2-level-2-science-products'},\n", " {'name': 'Microsoft',\n", " 'roles': [],\n", " 'url': 'https://planetarycomputer.microsoft.com'}]}" ] }, "execution_count": 54, "metadata": {}, "output_type": "execute_result" } ], "source": [ "collection.to_dict()" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "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](https://github.com/radiantearth/stac-browser).\n", "\n", "We can use this simple Python server (copied from [this gist](https://gist.github.com/acdha/925e9ffc3d74ad59c3ea)) to serve our our directory at port 5555:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import os\n", "from http.server import HTTPServer, SimpleHTTPRequestHandler\n", "\n", "os.chdir(root_path)\n", "\n", "\n", "class CORSRequestHandler(SimpleHTTPRequestHandler):\n", " def end_headers(self) -> None:\n", " self.send_header(\"Access-Control-Allow-Origin\", \"*\")\n", " self.send_header(\"Access-Control-Allow-Methods\", \"GET\")\n", " self.send_header(\"Cache-Control\", \"no-store, no-cache, must-revalidate\")\n", " return super(CORSRequestHandler, self).end_headers()\n", "\n", "\n", "with HTTPServer((\"localhost\", 5555), CORSRequestHandler) as httpd:\n", " httpd.serve_forever()" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Now we can browse our STAC Collection with stac-browser in a few different ways:\n", "1. Follow the [instructions](https://github.com/radiantearth/stac-browser/blob/main/local_files.md) for starting a stac-browser instance and point it at `http://localhost:5555/collection.json` to serve out the STAC.\n", "2. If you want to avoid setting up your own stac-browser instance, you can use the [STAC Browser Demo](https://radiantearth.github.io/stac-browser/) hosted by Radiant Earth: [https://radiantearth.github.io/stac-browser/#/http://localhost:5555/collection.json](https://radiantearth.github.io/stac-browser/#/http://localhost:5555/collection.json)\n", "\n", "To quit the server, use the `Kernel` -> `Interrupt` menu option." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Acknowledgements\n", "\n", "Credit to [sat-stac-landsat](https://github.com/sat-utils/sat-stac-landsat) from which a lot of this code was based." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.3" } }, "nbformat": 4, "nbformat_minor": 2 }