Hint

You can run this notebook in a live session with .

Get Copernicus DEM using EODAG

This tutorial will help you to search and download Copernicus DEM using EODAG

import os
import boto3
import json
import zipfile

from botocore.handlers import disable_signing

from eodag import EODataAccessGateway
from eodag.utils import ProgressCallback

Create a workspace directory eodag_workspace_copdem where all our files will live:

workspace = 'eodag_workspace_copdem'
if not os.path.isdir(workspace):
    os.mkdir(workspace)

Download grid index from AWS:

grid_zipfilename = "grid.zip"
progress_callback = ProgressCallback()
path_to_zipfile = os.path.join(os.path.abspath(workspace), grid_zipfilename)

s3_resource = boto3.resource("s3")
s3_resource.meta.client.meta.events.register(
    "choose-signer.s3.*", disable_signing
)
s3_bucket = s3_resource.Bucket("copernicus-dem-30m")

# file size
progress_callback.reset(total=sum([
    p.size for p in s3_bucket.objects.filter(Prefix=grid_zipfilename)
]))
# download
objects = s3_bucket.download_file(
    grid_zipfilename,
    path_to_zipfile,
    Callback=progress_callback,
)
progress_callback.close()

# unzip
with zipfile.ZipFile(path_to_zipfile, 'r') as zip_ref:
    zip_ref.extractall(workspace)

Init EODAG and add a custom provider for this local Copernicus DEM index:

dag = EODataAccessGateway()
dag.update_providers_config("""
cop_dem_aws:
    search:
        type: StacSearch
        api_endpoint: https://fake-endpoint
        metadata_mapping:
            id: '$.properties.id'
            title: '$.properties.id'
            downloadLink: 's3://copernicus-dem-30m/{id}'
    products:
        GENERIC_PRODUCT_TYPE:
            productType: '{productType}'
    download:
        type: AwsDownload
        base_uri:  https://fake-endpoint
        flatten_top_dirs: True
    auth:
        type: AwsAuth
        credentials:
""")

Convert downloaded geojson to EODAG serialized search result

dem_in = os.path.join(os.path.abspath(workspace), "dem30mGrid.json")
dem_out = os.path.join(os.path.abspath(workspace), "dem30mGrid_eodag.json")
with open(dem_in) as fin, open(dem_out, "w") as fout:
    data = json.load(fin)
    data["features"] = [dict(feat, **{
            "id": feat["properties"]["id"],
            "properties" : {
                "title": feat["properties"]["id"],
                "eodag_product_type": "COP_DEM",
                "eodag_provider": "cop_dem_aws",
                "eodag_search_intersection": {"type": "Polygon","coordinates": []},
                "downloadLink": "s3://copernicus-dem-30m/%s" % feat["properties"]["id"],
            }
        }) for feat in data["features"]
    ]
    json.dump(data, fout)
print("%s items written to %s" % (len(data["features"]), dem_out))

items = dag.deserialize_and_register(dem_out)
print("%s items loaded as SearchResult" % len(items))

64800 items written to /home/sylvain/workspace/eodag/examples/eodag_workspace_copdem/dem30mGrid_eodag.json
64800 items loaded as SearchResult

Filter on a specific area

Draw geometry on the map:

import ipyleaflet as ipyl
import ipywidgets as widgets

center = [43.53, 1.33] # Toulouse
zoom = 6
m = ipyl.Map(center=center, zoom=zoom)

def clear_m():
    global drawn_polygons
    drawn_polygons = set()

clear_m()

myDrawControl = ipyl.DrawControl(polygon={'shapeOptions':{'color':'#00F'}}, circlemarker={}, polyline={})

def handle_draw(self, action, geo_json):
    global drawn_polygons, m
    polygon=[]
    for coords in geo_json['geometry']['coordinates'][0][:-1][:]:
        polygon.append(tuple(coords))
    polygon = tuple(polygon)
    if action == 'created':
        drawn_polygons.add(polygon)
    elif action == 'deleted':
        drawn_polygons.discard(polygon)
        # clear found tiles
        for l in m.layers[1:]:
            m.remove_layer(l)

myDrawControl.on_draw(handle_draw)
m.add_control(myDrawControl)
m

Filter using drawn geometry and show on map:

from shapely.geometry import Polygon

filtered_by_geom = []
if drawn_polygons:
    drawn_polygon = Polygon(list(next(iter(drawn_polygons))))

    # Filter products
    filtered_by_geom = items.filter_overlap(
        intersects=True,
        geometry=drawn_polygon
    )
    print("Found %s products intersecting given geometry" % len(filtered_by_geom))

    # show on map
    filtered_by_geom_layer = ipyl.GeoJSON(data=filtered_by_geom.as_geojson_object(), style=dict(color='green'))
    m.add_layer(filtered_by_geom_layer)

Found 3 products intersecting given geometry

Download the filtered DEM files to workspace

paths = dag.download_all(
    filtered_by_geom,
    outputs_prefix=workspace,
)
paths

['/home/sylvain/workspace/eodag/examples/eodag_workspace_copdem/Copernicus_DSM_COG_10_N42_00_E001_00_DEM',
 '/home/sylvain/workspace/eodag/examples/eodag_workspace_copdem/Copernicus_DSM_COG_10_N43_00_E001_00_DEM',
 '/home/sylvain/workspace/eodag/examples/eodag_workspace_copdem/Copernicus_DSM_COG_10_N42_00_E002_00_DEM']