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
from eoreader.stac import OPTICAL_STAC_EXTENSIONS
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',
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)
2023-03-22 16:25:33,722 - [INFO] - *** LC08_L1TP_200030_20201220_20210310_02_T1.tar ***
2023-03-22 16:25:36,723 - [WARNING] - No quicklook found in 20201220T104856_L8_200030_OLI_TIRS
2023-03-22 16:25:37,565 - [INFO] - *** LT05_L1TP_200030_20111110_20200820_02_T1.tar ***
2023-03-22 16:25:39,111 - [INFO] - *** S2A_MSIL1C_20191215T110441_N0208_R094_T30TXP_20191215T114155.SAFE ***
2023-03-22 16:25:40,119 - [WARNING] - 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',
'https://stac-extensions.github.io/projection/v1.1.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:
08
nir
NARROW_NIR:
8A
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 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")
sara skipped: could not be loaded from user configuration
meteoblue skipped: could not be loaded from user configuration
cop_dataspace skipped: could not be loaded from user configuration
# 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-03-22T16:25:36.789737Z',
'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()