Get Copernicus DEM using EODAG
Contents
Get Copernicus DEM using EODAG#
This tutorial will help you to search and download Copernicus DEM using EODAG
[1]:
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:
[2]:
workspace = 'eodag_workspace_copdem'
if not os.path.isdir(workspace):
os.mkdir(workspace)
Download grid index from AWS:
[3]:
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:
[4]:
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#
[5]:
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:
[6]:
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:
[9]:
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#
[10]:
paths = dag.download_all(
filtered_by_geom,
outputs_prefix=workspace,
)
paths
[10]:
['/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']