Python – How to Get Copernicus DEM GeoTIFFs for a Bounding Box Using Python

amazon s3amazon-web-servicesdemelevationpython

Background

According to AWS's Copernicus Digital Elevation Model page,

The Copernicus DEM is a Digital Surface Model (DSM) which represents the surface of the Earth including buildings, infrastructure and vegetation. We provide two instances of Copernicus DEM named GLO-30 Public and GLO-90. […] Data is provided as Cloud Optimized GeoTIFFs.

Aim

I would like to access these GeoTIFFs, using Python. I would like to be able to specify a bounding box, and get those GeoTIFFs which overlap or intersect with that bounding box.

My attempts

1) Access AWS directly

After installing boto3 (with pip3 install boto3), I do, relying on this answer to the question Can I use boto3 anonymously?:

import boto3
from botocore import UNSIGNED
from botocore.client import Config

s3 = boto3.client('s3', region_name='eu-central-1', config=Config(signature_version=UNSIGNED))

Then I query for list of objects in the bucket, using the second line of this answer to the question Use boto3 to download from public bucket:

objects = s3.list_objects(Bucket='copernicus-dem-30m')

I get a dict, with keys:

objects.keys()
# outputs: dict_keys(['ResponseMetadata', 'IsTruncated', 'Marker', 'Contents', 'Name', 'Prefix', 'MaxKeys', 'EncodingType'])

I access Contents, which is a list, so viewing the first element:

objects['Contents'][0]

I get:

{'ETag': '"e0c8b05f7999eedd83ba85e5c94cb670"',
 'Key': 'Copernicus_DSM_COG_10_N00_00_E006_00_DEM/Copernicus_DSM_COG_10_N00_00_E006_00_DEM.tif',
 'LastModified': datetime.datetime(2020, 11, 25, 7, 51, 3, tzinfo=tzlocal()),
 'Size': 3706379,
 'StorageClass': 'STANDARD'}

Viewing objects['Contents'][0]['Key'] gives me Copernicus_DSM_COG_10_N00_00_E006_00_DEM/Copernicus_DSM_COG_10_N00_00_E006_00_DEM.tif, which is a string. So this didn't got me much closer to accessing the actual TIFF data in a Python notebook.

2) Using opentopography.org

Googling "search in GLO-30 bucket" lead me to OpenTopography's search map. After selecting a small region on the map & filling other other details (like email), and waiting a few moments, I get to:

enter image description here

where I can download a .tar.gz file. Upon extraction, I find a .tif file.

I believe this is an example of the files I am looking for, but this method does not use Pyhton at all, and it is not suitable for automated processes.

Question

How do I access the TIF files using Python?

I am looking for a similar thing to this this:

geotiffs = client.search_geotiffs(bounding_box=[59,10,60,11])

where the returned geotiffs is a list of tiff files from the GLO-30, intersecting or overlapping with the bounding box defined by points 59N10E & 60N11E, for example.

Best Answer

Try our open source python package available through PyPI or conda-forge that merges tiles from the open AWS registry of Copernicus 30 meter DEM. (Obviously: I am one of the library developers). We are using this for InSAR processing, so there are few more basic transformations that the library can do, which you can read about on the homepage.

https://github.com/ACCESS-Cloud-Based-InSAR/dem-stitcher

bounds = [-119.085, 33.402, -118.984, 35.435]
X, p = stitch_dem(bounds,
                  dem_name='glo_30',
                  dst_ellipsoidal_height=False,
                  dst_area_or_point='Point')
# X is an m x n numpy array
# p is a dictionary (or a rasterio profile) including relevant GIS metadata

Hope this helps. Lots of information on the page. Can also utilize other global tiles such as NASADEM.

Related: You could also download one of the zipped geojsons from the library here, unzip it, load it into QGIS, and then click the URLs for tiles you are interested. That seems like that might be good too. For Copernicus 30 meter tiles, the file name is glo_30.geojson.zip.

Related Question