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 logging

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
from eoreader.stac import OPTICAL_STAC_EXTENSIONS

Create logger#

logger = logging.getLogger("eoreader")
logger.setLevel(logging.INFO)

# create console handler and set level to debug
ch = logging.StreamHandler()
ch.setLevel(logging.INFO)

# create formatter
formatter = logging.Formatter('%(message)s')

# add formatter to ch
ch.setFormatter(formatter)

# add ch to logger
logger.addHandler(ch)

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',
    stac_extensions=OPTICAL_STAC_EXTENSIONS,
    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)
*** LC08_L1TP_200030_20201220_20210310_02_T1.tar ***
No quicklook found in 20201220T104856_L8_200030_OLI_TIRS
*** LT05_L1TP_200030_20111110_20200820_02_T1.tar ***
*** S2A_MSIL1C_20191215T110441_N0208_R094_T30TXP_20191215T114155.SAFE ***
No quicklook found in 20191215T110441_S2_T30TXP_L1C_114155
# See of what the STAC object of the last product is composed
prod.stac
STAC Item attributes:
	id: 20191215T110441_S2_T30TXP_L1C_114155
	constellation: sentinel-2
	gsd: 10.0
	datetime: 2019-12-15 11:04:41
	geometry:
		{'coordinates': (((-1.7679672999632903, 43.25815531863286, 0.0),
		                  (-1.747519127782274, 44.24655093448751, 0.0),
		                  (-0.3733540263, 44.22320973997534, 0.0),
		                  (-0.416191720429816, 43.23560095969715, 0.0),
		                  (-1.7679672999632903, 43.25815531863286, 0.0)),),
		 'type': 'Polygon'}
	bbox: 
		[-1.7679672999632903, 43.23560095969715, -0.3733540263, 44.24655093448751]
	stac_extensions: 
		['https://stac-extensions.github.io/eo/v1.0.0/schema.json',
		 'https://stac-extensions.github.io/projection/v1.0.0/schema.json',
		 'https://stac-extensions.github.io/view/v1.0.0/schema.json']
	properties - tilename: T30TXP

	Electro-Optical STAC Extension attributes:
		eo:cloud_cover: 0.6668
		eo:bands:
			COASTAL_AEROSOL:
				01
				coastal
			BLUE:
				02
				blue
			GREEN:
				03
				green
			RED:
				04
				red
			VEGETATION_RED_EDGE_1:
				05
				rededge
			VEGETATION_RED_EDGE_2:
				06
				rededge
			VEGETATION_RED_EDGE_3:
				07
				rededge
			NIR:
				8A
				nir
			NARROW_NIR:
				08
				nir08
			WATER_VAPOUR:
				09
				nir09
			CIRRUS:
				10
				cirrus
			SWIR_1:
				11
				swir16
			SWIR_2:
				12
				swir22
	Projection STAC Extension attributes:
		proj:epsg: 32630
		proj:geometry: 
			{'coordinates': (((600000.0, 4790220.0, 0.0),
			                  (600000.0, 4900020.0, 0.0),
			                  (709800.0, 4900020.0, 0.0),
			                  (709800.0, 4790220.0, 0.0),
			                  (600000.0, 4790220.0, 0.0)),),
			 'type': 'Polygon'}
		proj:bbox: 
			[600000.0, 4790220.0, 709800.0, 4900020.0]
		proj:centroid: 
			{'lat': 43.74293327815072, 'lon': -1.0762579548348727}
		proj:shape: [10980, 10980]
		proj:transform: 
			Affine(10.0, 0.0, 600000.0,
			       0.0, -10.0, 4900020.0)
	View STAC Extension attributes:
		view:sun_azimuth: 167.543494280222
		view:sun_elevation: 67.991496343423
# 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>

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 INFO
setup_logging(verbose=2)

# 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")
2022-10-10 14:13:02,388 eodag.config                     [INFO    ] stac_http_provider: unknown provider found in user conf, trying to use provided configuration
# Query every product from inside the catalog
all_products, _ = dag.search()
2022-10-10 14:13:02,496 eodag.core                       [INFO    ] No product type could be guessed with provided arguments
2022-10-10 14:13:03,172 eodag.core                       [INFO    ] Searching product type 'None' on provider: stac_http_provider
2022-10-10 14:13:03,280 eodag.plugins.crunch.filter_overlap [WARNING ] geometry not found in cruncher arguments, filtering disabled.
2022-10-10 14:13:03,285 eodag.core                       [INFO    ] Found 3 result(s) on provider 'stac_http_provider'
# 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)
2022-10-10 14:13:03,301 eodag.core                       [INFO    ] No product type could be guessed with provided arguments
2022-10-10 14:13:03,302 eodag.core                       [INFO    ] Searching product type 'None' on provider: stac_http_provider
2022-10-10 14:13:03,413 eodag.plugins.crunch.filter_date [INFO    ] Finished filtering products. 1 resulting products
2022-10-10 14:13:03,414 eodag.plugins.crunch.filter_overlap [INFO    ] Finished filtering products. 1 resulting products
2022-10-10 14:13:03,417 eodag.core                       [INFO    ] Found 1 result(s) on provider 'stac_http_provider'
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': '2022-10-10T14:12:56.703962Z',
 '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()