EODAG - DestinE Data Lake Provider - DEDL Notebook Gallery

🚀 Launch in JupyterHub

Prerequisites:

To search and access DEDL data a DestinE user account is needed

References:

EODAG - Earth Observation Data Access Gateway

EODAG - DEDL

EODAG API documentation

Credit:

Code used in this context takes inspiration from the Python API User Guide produced by CS Group.

EODAG is a command line tool and a Python package for searching and downloading earth observation data via a unified API.

This tutorial demonstrates how to use EODAG to search and access DEDL data. The notebook uses the DEDL provider in EODAG to access data via Python code.

Setup: EODAG configuration to use the provider DEDL .
Search: search DEDL data, we search for Sentinel-3 data.
Filter: filter DEDL data.
Download: download DEDL data.

The complete guide on how to use EODAG Python API is available via https://eodag.readthedocs.io/en/stable/api_user_guide.html.

Note:

Please note that the two factor authentication (2FA) is still not implemented in EODAG. The users who have enabled 2FA on DESP will not be able to run this notebook.

Setup¶

In this section, we set:

The output_dir, the directory where to store downloaded products.
The DEDL credentials, you’ll be asked to enter your DEDL credentials.
The search timeout, it is of 60 seconds to avoid any unexpected errors because of long running search queries.

import os
from getpass import getpass

workspace = 'eodag_workspace'
if not os.path.isdir(workspace):
    os.mkdir(workspace)
    
os.environ["EODAG__DEDL__DOWNLOAD__OUTPUT_DIR"] = os.path.abspath(workspace)
#os.environ["EODAG__DEDL__DOWNLOAD__OUTPUTS_PREFIX"] = os.path.abspath(workspace)

os.environ["EODAG__DEDL__PRIORITY"]="10"
os.environ["EODAG__DEDL__SEARCH__TIMEOUT"]="60"


DESP_USERNAME = input("Please input your DESP username or email: ")
DESP_PASSWORD = getpass("Please input your DESP password: ")

os.environ["EODAG__DEDL__AUTH__CREDENTIALS__USERNAME"]=DESP_USERNAME
os.environ["EODAG__DEDL__AUTH__CREDENTIALS__PASSWORD"]=DESP_PASSWORD

Import EODAG and list available products on DEDL¶

We now need to import the EODataAccessGateway class. The class is going to take care of all the following operations.

We can start listing the products available using dedl as provider.


from eodag import EODataAccessGateway, setup_logging
dag = EODataAccessGateway()

setup_logging(0)

print("\033[1mID - Title\033[0m")

# Rows
for pt in dag.list_product_types("dedl"):
    print(f'{pt["ID"]} - {pt["title"]}')

Search¶

To search we use the search method passing the ID of our dataset of interest and a geo-time filter.

The search method returns a SearchResult object that stores the products obtained from a given page (default: page=1) and a given maximum number of items per page (default: items_per_page=20).

In the following cell, we change the default value of items_per_page and define the search criteria to retrieve Sentinel-2 MSI Level-2 images over Sicily, first days of July 2024. Our goal is to check whether any effects of Mount Etna’s eruptions during that period are visible in the Sentinel-2 imagery.

search_criteria = {
    "provider":"dedl",
    "productType": "EO.ESA.DAT.SENTINEL-2.MSI.L2A",
    "start": "2024-07-04T07:00:00.00Z",
    "end": "2024-07-08T07:00:00.00Z",
    "geom": {"lonmin": 12, "latmin": 37, "lonmax": 16, "latmax": 39},
    "count": True,
    "items_per_page": 50
}

products_first_page = dag.search(**search_criteria)

Results are stored in a ‘SearchResult’ object that contains the details on the single search result.

products_first_page

It is possible to list the metadata associated with a certain product, we choose the first one returned [0], and look into it.

one_product = products_first_page[0]
one_product.properties.keys()

one_product.properties['cloudCover']

Filter¶

EODAG can filter the search result. We can then refine our initial search without asking the provider again. Products can be filtered according to their properties or also with finer geometry filters.

The following example shows how to filter products to keep only those whose cloud coverage is less than 20%. And then restrict the results to products containing a smaller area over the mount Etna.

Let’s define now a smaller area around the mount Etna and a function to see the area on a map together with the results

from eodag.crunch import FilterProperty
from eodag.crunch import FilterOverlap
import shapely
import folium
from shapely.geometry import Polygon

small_geom = Polygon([[15.1, 37.7], [15.5, 37.7], [15.1, 37.75], [15.1, 37.75], [15.1, 37.7]])

smaller_area = {"lonmin": 15.1, "latmin": 37.7, "lonmax": 15.5, "latmax": 37.75}

search_geometry = shapely.geometry.box(
    smaller_area["lonmin"],
    smaller_area["latmin"],
    smaller_area["lonmax"],
    smaller_area["latmax"],
)


def create_search_result_map(search_results, extent):
    """Small utility to create an interactive map with folium
    that displays an extent in red and EO Producs in blue"""
    fmap = folium.Map([38, 14], zoom_start=7)
    folium.GeoJson(
        extent,
        style_function=lambda x: dict(color="red")
    ).add_to(fmap)
    folium.GeoJson(
        search_results
    ).add_to(fmap)
    return fmap

# Crunch the results
filtered_results = products_first_page.crunch(FilterProperty({"cloudCover": 20, "operator" : "lt"}))

print(f"Got now {len(filtered_results)} products after filtering by cloudCover.")

filtered_products = filtered_results.crunch(
    FilterOverlap(dict(contains=True)),
    geometry=small_geom
)
print(f"Got now {len(filtered_products)} products after filtering by geometry.")

Let’s use the function defined to see the area defined on a map (red) together with the initial results (blue) filtered by cloud coverage and geometry (green).

fmap = create_search_result_map(products_first_page, search_geometry)
# Create a layer that represents the filtered products in green
folium.GeoJson(
    filtered_products,
    style_function=lambda x: dict(color="green")
).add_to(fmap)
fmap

Download¶

Before downloading any product, it can be useful to have a quick look at them.

import matplotlib.pyplot as plt
import matplotlib.image as mpimg

quicklooks_dir = os.path.join(workspace, "quicklooks")
if not os.path.isdir(quicklooks_dir):
    os.mkdir(quicklooks_dir)

fig = plt.figure(figsize=(20, 40))
for i, product in enumerate(filtered_products, start=1):

    # This line takes care of downloading the quicklook
    quicklook_path = product.get_quicklook()
    
    img = mpimg.imread(quicklook_path)
    ax = fig.add_subplot(8, 2, i)
    ax.set_title(product.properties['dedl:beginningDateTime'] + "TILE: " +product.properties['dedl:tileIdentifier'])
    plt.imshow(img)
plt.tight_layout()

The quicklook shows effectively the ash plume caused by the eruptions.

EOProducts can be downloaded individually. The last image is going to be downloaded.

product_to_download = filtered_products[-1]
product_path = dag.download(product_to_download)
product_path

The location property of this product now points to a local path.

product_to_download.location