Skip to main content

STAC and GDAL

This Notebook (or rather the GDAL commands created) must be run in a terrabyte environment with direct access to the satellite data. See STAC, xarray and dask for an example of how to set up such an environment in the browser.

Example: Extracting STAC data with GDAL​

import os

In order to avoid a bug that prevents GDAL from properly interpreting local file:// URIS in the terrabyte system, make sure you have GDAL Version 3.8.0 or higher installed. (If this is not available to you, consider using the master branch installabe from the gdal-master conda channel).

!gdalinfo --version

GDAL 3.8.0dev-c95490018a-dirty, released 2023/07/21

# GDAL directly uses the search API of STAC
stac_search_url="https://stac.terrabyte.lrz.de/public/api/search"
# STAC Collection to search in
collection="cop-dem-glo-90"
# Name of the Asset in the collection above we want do download
asset="dem"

# Area of interest
bboxMunich = [11.089325,47.893787,11.923599,48.329322]
#so we can switch bounding box easily
bbox=bboxMunich

# Extra gdal options to pass
# Here we activate BigTiff support, LZW compression, and
# assign an ESRI-conform name to the band
gdalOptions=f"-co COMPRESS=LZW -co BIGTIFF=YES -mo band_1={asset}"
outputFilename="gdal_COPDEM_example.tif"

Specifying the bounding box and band name as search parameters to STAC ensures only those results are returned that intersect and match the band we want to use.

gdal_inputURl=f'STACIT:"{stac_search_url}?collections={collection}&bbox={",".join(str(b) for b in bbox)}":asset={asset}'
infoCmd=f"gdalinfo '{gdal_inputURl}'"
print(infoCmd)

gdalinfo 'STACIT:"https://stac.terrabyte.lrz.de/public/api/search?collections=cop-dem-glo-90&bbox=11.089325,47.893787,11.923599,48.329322":asset=dem'

os.system(infoCmd)

gdalinfo output

We specify the bounding box again to GDAL itself to clip the output exactly to the bounding box. Be aware that GDAL expects a different order of coordinates.

gdal_projwin=f"-projwin {' '.join(str(bbox[i]) for i in [0,3,2,1])}"
gdal_projwin

'-projwin 11.089325 48.329322 11.923599 47.893787'

translateCmd=f"gdal_translate '{gdal_inputURl}' {outputFilename} {gdal_projwin} {gdalOptions}"
print(translateCmd)

gdal_translate 'STACIT:"<https://stac.terrabyte.lrz.de/public/api/search?collections=cop-dem-glo-90&bbox=11.089325,47.893787,11.923599,48.329322>":asset=dem' gdal_COPDEM_example.tif -projwin 11.089325 48.329322 11.923599 47.893787 -co COMPRESS=LZW -co BIGTIFF=YES -mo band_1=dem

os.system(translateCmd)

Input file size is 1201, 2401

0...10...20...30...40...50...60...70...80...90...100 - done.