STAC#

Let’s use EOReader to create SpatioTemporal Asset Catalog (STAC) items.

Note: This is experimental for now, use it at your own risk !

Warning: You will need to install pystac[validation], folium and eodag (version != 2.6.0) to run this notebook

Imports#

# Imports
import os

import pystac
import geopandas as gpd
from tempfile import TemporaryDirectory
from shapely.geometry import mapping

from eodag import setup_logging
from eodag.api.core import EODataAccessGateway

from eoreader.reader import Reader

Create logger#

# Create logger
import logging
from sertit import logs

logger = logging.getLogger("eoreader")
logs.init_logger(logger)

Linking some data#

Let’s take 3 products covering approximately the same area (over DAX city in France):

  • One Landsat-8 OLI-TIRS collection 2

  • One Landsat-5 TM collection 2

  • One Sentinel-2 L1C

prod_folder = os.path.join("/home", "data", "DATA", "PRODS")
paths = [
    # Landsat-8 OLI-TIRS collection 2
    os.path.join(prod_folder, "LANDSATS_COL2", "LC08_L1TP_200030_20201220_20210310_02_T1.tar"),
    # Landsat-5 TM collection 2    
    os.path.join(prod_folder, "LANDSATS_COL2", "LT05_L1TP_200030_20111110_20200820_02_T1.tar"),
    # Sentinel-2 L2A
    os.path.join(prod_folder, "S2", "PB 02.07+", "S2A_MSIL1C_20191215T110441_N0208_R094_T30TXP_20191215T114155.SAFE"),
]

Create STAC catalog#

Create a STAC catalog and add 3 STAC items to it.

# Create the reader
reader = Reader()

# Work in a temporary directory
tmp = TemporaryDirectory()

# Create STAC catalog
catalog_path = os.path.join(tmp.name, "catalog.json")
catalog = pystac.Catalog(
    id='SERTIT_101',
    description="SERTIT's Catalog",
    title='SERTIT Catalog',
    href=catalog_path
)

# Add all the products into the STAC catalog
for path in paths:
    logger.info(f"*** {os.path.basename(path)} ***")

    # Open the product
    prod = reader.open(path, remove_tmp=True)

    # Get item
    item = prod.stac.create_item()

    # Add item to catalogue
    catalog.add_item(item)

2023-05-31 13:19:02,675 - [INFO] - *** LC08_L1TP_200030_20201220_20210310_02_T1.tar ***
2023-05-31 13:19:05,884 - [WARNING] - No quicklook found in 20201220T104856_L8_200030_OLI_TIRS
2023-05-31 13:19:07,229 - [INFO] - *** LT05_L1TP_200030_20111110_20200820_02_T1.tar ***
2023-05-31 13:19:08,898 - [INFO] - *** S2A_MSIL1C_20191215T110441_N0208_R094_T30TXP_20191215T114155.SAFE ***

# Save catalog
catalog.describe()
catalog.normalize_and_save(tmp.name, catalog_type=pystac.CatalogType.SELF_CONTAINED)

* <Catalog id=SERTIT_101>
  * <Item id=20201220T104856_L8_200030_OLI_TIRS>
  * <Item id=20111110T103612_L5_200030_TM>
  * <Item id=20191215T110441_S2_T30TXP_L1C_114155>

list(catalog.get_items())[0]

Query the catalog#

EODAG is an opensource python library that implements STAC and allows you to query your local STAC catalog.
Look at here for a detailed tutorial.

# Create an EODAG custom STAC provider
dag = EODataAccessGateway()

# Set EODAG logging level to WARNING
setup_logging(verbose=1)

# Add the custom STAC provider, exactly like in the tutorial mentioned above
dag.update_providers_config("""
stac_http_provider:
    search:
        type: StaticStacSearch
        api_endpoint: %s
    products:
        GENERIC_PRODUCT_TYPE:
            productType: '{productType}'
    download:
        type: HTTPDownload
        base_uri: %s
        flatten_top_dirs: True
        outputs_prefix: %s
""" % (catalog_path, tmp.name, tmp.name))

# Set the custom STAC provider as preferred
dag.set_preferred_provider("stac_http_provider")

# Query every product from inside the catalog
all_products, _ = dag.search()

# Load an AOI
aoi_path = os.path.join("/home", "data", "DATA", "AOIs", "DAX.geojson")
aoi = gpd.read_file(aoi_path)
aoi_geojson = mapping(aoi.geometry.values[0])

# Query spatially with the AOI and temporally with a time period
query_args = {"start": "2020-05-01", "end": "2022-05-06", "geom": aoi.geometry.values[0]}
query_products, _ = dag.search(**query_args)

query_products[0]

EOProduct(id=20201220T104856_L8_200030_OLI_TIRS, provider=stac_http_provider)

query_products[0].assets['Blue']

{'href': 'tar+file:///home/data/DATA/PRODS/LANDSATS_COL2/LC08_L1TP_200030_20201220_20210310_02_T1.tar!LC08_L1TP_200030_20201220_20210310_02_T1_B2.TIF',
 'type': 'image/tiff; application=geotiff',
 'title': 'Blue',
 'description': 'Bathymetric mapping, distinguishing soil from vegetation and deciduous from coniferous vegetation',
 'eoreader_name': 'BLUE',
 'eo:bands': [{'name': 'Blue',
   'common_name': 'blue',
   'description': 'Bathymetric mapping, distinguishing soil from vegetation and deciduous from coniferous vegetation',
   'center_wavelength': 0.48,
   'full_width_half_max': 0.06}],
 'created': '2023-05-31T13:19:05.951933Z',
 'updated': '2023-05-31T13:19:05.951941Z',
 'start_datetime': '2020-12-20T10:48:56Z',
 'end_datetime': '2020-12-20T10:48:56Z',
 'instruments': ['OLI-TIRS'],
 'constellation': 'landsat-8',
 'gsd': 30,
 'roles': ['reflectance']}

Display the results#

We can use folium to display the results geometry over a map.

import folium

# Create a map zoomed over the search area
fmap = folium.Map((43.2, -1.05), zoom_start=7)

# Add a layer green layer for the query over the AOI
folium.GeoJson(
    data=all_products.as_geojson_object(),
    tooltip = "All products stored in the catalog",
    style_function=lambda x: {'color': 'green'}
).add_to(fmap)

# Add a layer green layer for the query over the AOI
folium.GeoJson(
    data=query_products.as_geojson_object(),
    tooltip = "Retrieved products with the query",
    style_function=lambda x: {'color': 'red'}
).add_to(fmap)

# Add a layer blue layer for the AOI
folium.GeoJson(
    data=aoi_geojson,
    tooltip = "DAX AOI",
    style_function=lambda x: {'color': 'blue'}
).add_to(fmap)

fmap
Make this Notebook Trusted to load map: File -> Trust Notebook
# Clean the tmp directory
tmp.cleanup()