Framework project: Digital Elevation Models (DEMs)

Process chain for Copernicus DEM

Introduction

Processing a global DEM all the way from downloading tiles to create coherent landscape indices and metrics, is a typical task for which Karttur’s GeoImagine Framework was built. Once the DEM is imported and organized within the Framework, it is easy to test different algorithms and visualizations. If you have access to some ground data you can also apply the Framework for comparing the results of different algorithms and parameters against ground data, and thus select the most appropriate landscape model.

Prerequisits

You must have installed Karttur’s GeoImagine Framework. You also need a registered project that covers the geographic regions of the DEM you want to process. The latter can, however, be done as part of the project itself.

Overview

This post summarises the processes chain for Copernicus DEM in Karttur’s GeoImagine Framework. Each step in the process chain is explained in more detail in linked posts. The Copernicus DEM in 90 m spatial resolution is used as an example.

A note on hydrological corrections

Before running the DEM processing you need to decide how you want to manage hydrological errors in your DEM. If there are no errors, then you do not need to worry, but usually there are at least both pits and flat area. Both of these problems primarily affects the extraction of river basins. Different algorithms handle these problems differently. The simplest solution is to fill up all pits completely, with or without tilting or “burn in” the flow across flat surfaces. This, however, causes the, normally a-priori known, relative DEM accuracy to be replaced by an unknown error that might, or might not, have spatial or contextual biases. If the DEM is going to be solely used for flow routing this is perhaps the best method.

If there are independent vector datasets on stream networks and channels available, you can also burn a drainage pattern into the DEM. The geometric accuracy must match that of the DEM, and the vector network contiguity must be correct. For wide rivers and lakes, burning a channel also causes problems with a preferential flow line with subsequent errors when calculating e.g. flooding and evapotranspiration.

Another alternative is to only fill up minor errors and errors that can be logically identified. The latter category would then include depressed channels and negative elevations along shorelines - and this is how it is done in this manual. But then you must also apply a flow routing algorithm that can pass depressions on the fly. The GRASS GIS module r.watershed does just that. And r.watershed is used for the DEM modelling in this post. To fill up all depressions, you can instead use the GRASS GIS module r.terraflow, or any of the DEM filling algorithms of SAGA GIS. If you apply hydrological DEM corrections, you must also decide which DEM version (corrected or non-corrected) to use for extracting other indexes, like slope, curvature etc.

Process chain

The process chain is available at Kartturs GitHub account and also found under the Hide/show button below. The process chain is an ordinary text file linking to a series of json files. Empty rows and rows starting with a “#” are ignored. When called from Karttur’s GeoImagine Framework, the processes of each json file are sequentially run. In the process chain below you can click on each json file name to see the complete content. The following sections also explain in more detail the processes of each json command file.

Layout

To produce symbolised and labeled maps using the Framework, you have to define some Layout features and parameters. All map layout is done using unsigned byte range maps with values ranging from 0 to 255, Values 251 to 255 are reserved for specific purposes and not allowed for representing map values. The Framework Layout processes include:

Create Scaling

Framework process: CreateScaling

Json command file: 0001_CreatesScaling_DEM_v090.json

The CreateScaling process converts any input raster into an unsigned byte (0-255) range file. The byte file is then used for assigning a color ramp (defined in the process AddRasterPalette, see next section).

Scaling is thus not defined on the fly, instead each map composition must be associated with a predefined scaling. Once defined, this scaling can not be changed. You can, however, change the color ramping by setting another raster palette.

Add Raster Palette

Framework process: AddRasterPalette

Json command file: 0002_AddRasterPalette_DEM_v090.json

In contrast to scaling, a palette can be freely set for any map composition (but always using the same scaling). The palette to use for any map layout is defined when calling the process that generates the layout. Thus the palette must exist beforehand, created by using the process AddRasterPalette.

Create Legend

Framework process: CreateLegend

Json command file: 0003_createlegends_DEM_v090.json

With the scaling set and a palette defined, you can define a legend - that is how the combined scaling and palette will be presented in graphical layouts.

Export Legend

Framework process: ExportLegend

Json command file: 0005_exportlegend_DEM_v090.json

Once created you can export a legend. All legends are exported as png, csv and pdf. Exporting a legend requires that you have installed the application inkscape as outlined in the post Python virtual environment. The exported legends are stored under a director legends on the destination disk/volume. The image represented the raster palette is under a subfolder images (saved as a png image) and the legends with text added are under a subfolder layouts, and saved as three different versions:

  • svg
  • png
  • pdf

Create a project

0020_ManageDefRegProj_copdem-ease2n.json

Data specific processing requires that your user has a predefined project that includes a predefined tract, where the tract must be linked to a default region. All spatial processing that you process from that user and with that project will then operate on system tiles linked to parent default region of your project tract. The whole processing is set up in thin manner to minimise double processing.

Unless you have done so before you must setup the project and tract to use for processing the Copernicus DEM data. The example is for processing the DEM for the entire basins draining into Arctic waters. Thus the example assign the default region northlandease2n as the parent region. Both the user specific project and its region are called karttur-northlandease2n.

Search, download and unzip

Search and download differ for almost all datasets. Details for searching and downloading Copernicus DEM (and other DEMs) is in the post Digital Elevation Models (DEMs). Unzipping is only required if the retrieved data is in a compressed zip format.

Framework process: SearchCopernicusProducts

Json command file: 0100_SearchCopernicusProducts_CopDEM-90m.json

The search for the Copernicus DEM 90 m product is done using the special process SearchCopernicusProducts. The search requires that you installed wget, see the post on Python virtual environment, that also covers installing wget.

The process SearchCopernicusProducts results in an html file with url links to all the resources to download. For the process defined in 0100_SearchCopernicusProducts_CopDEM-90m.json, the result is stored in the file copernicusDEM90.html

Download

Framework process: DownloadCopernicus

Json command file: 0110_DownloadCopernicus_CopDEM-90m.json

Also downloading the Copernicus DEM 90 m version is defined with a customized process, DownloadCopernicus. This process reads the html file above (copernicusDEM90.html), extracts the url tag (<a>), looks for a local copy of the DEM tile, and if not found (or overwrite set to true) downloads the tile. The downloading is predefined to a path on the destination disk and can not be changed. If you want to store that data somewhere else you have to move them manually, but the system will then lose track of them. The original tiles are not registered in the database. It is thus recommended that you leave the original data in the predefined location (for the data used here the default path is /DISKPATH/DOWNLOADS/CopernicusDem90).

Unzip

Framework process: UnZipRawData

Json command file: 0120_UnZipRawData_CopDEM-90m.json

Downloaded data that is retrieved as zip files can be unzipped with the command UnZipRawData. This is a generic process, and you need to give the full path to both the source (zip) and target (expanded) datasets. The disk (volume) paths are as usual given through the parameters srcpath and dstpath, but additioally you must also give the directories on the disks for both srcdir and dstdir under parameters. If srcdir is a file, the script assumes that it is a list of the files to unzip; if srcdir is a directory the script will explode all zip files in that directory.

For the example here I created a file listing all the zip files to explode (CopernicusDem90_rawtiles.csv). In MacOS or Linux you create the csv list by piping the ls command to a file:

$ ls *.zip > CopernicusDem90_rawtiles.csv

Import

The retrieved data is not yet organized as a dataset in the Framework database. A tiled dataset, like CopDEM, must be prepared as a virtual full coverage dataset prior to importing (organization) and then retiling into one, or more, projection system(s) that are used for the actual spatial data processing in the Framework. For the CopDEM dataset, the sequence of processes thus becomes:

  • identification of raw tile bounding boxes (optional),
  • mosaic raw tiles,
  • organize (import into Framework), and
  • tiling to Framework projection system.

Bounding Boxes for raw tiles

Framework process: BBoxTiledRawData

Json command file: 0122_BBoxTiledRawData_gs-CopDEM-90m.json

The process BBoxTiledRawData creates bounding boxes for each raw tile. In the example used here, each tile to include is listed in a csv file, created in a similar way as the list of zip files in the previous section. Running this process is optional, and intended for allowing a visual inspection of the raw data source.

The process generates a single vector data source with all the included tiles.

Bounding boxes for the Copernicus DEM 90 m raw tiles.

As clearly seen in the figure above, the raw tiles are of arbitrary sizes with geographical coverage restricted to land. The next steps in the process chain converts these more arbitrary tiles to regular grids fitted to the projection system where you want to use the DEM.

Mosaic raw tiles

Framework process: MosaicAncillary

Json command file: 0125_MosaicAncillary_CopDEM-90m.json

Downloaded data that are delivered as arbitrary tiles need to be mosaicked together before importing to the Framework. (SMAP, MODIS and Sentinel data do not need mosaicking of raw tiles as these datasets either are downloaded in predefined tiling systems or as complete global datasets). The mosaicking must be done as virtual (.vrt) dataset. To actually create a true global mosaic at e.g 90 m, as in this example, would take a long time and also take up unnecessary much disk space. The raw tile mosaic will represent the raw source data imported (organized) into the Framework database. The process MosaicAncillary does not register the imported (organized in the terminology of the Framework) dataset into the database, that is done in the next step.

Virtual mosaic of the Copernicus DEM 90 m raw tiles.

Organize (import and register) ancillary dataset

Framework process: OrganizeAncillary

Json command file: 0160_OrganizeAncillary_CopDEM-90m

With the Copernicus data downloaded, unzipped and mosaicked, the data can be organized for (or imported to) the Framework. For the Copernicus DEM data, the process OrganizeAncillary as defined in the json object below only imports the virtual mosaic linked to the original (downloaded and unzipped) tiles.

Once you have run the json command, you can check in the Framework database, and you should find that the ancillary schema contains a layer with the compid dem_copdem. Apart from being defined as an individual layer, the dataset is also registered as a product (table: compprod) and a composition (table: compdef):

SELECT * FROM ancillary.layer WHERE compid = 'dem_copdem';
SELECT * FROM ancillary.compprod WHERE compid = 'dem_copdem';
SELECT * FROM ancillary.compdef WHERE compid = 'dem_copdem';

Tiling to Framework projection system

Framework process: TileAncillaryRegion

Json command file: 0180_TileAncillaryRegion_CopDem-90m.json

All basic processing in the Framework is done using predefined tiling systems. The process TileAncillaryRegion constructs tiles while also reprojecting the source (ancillary) data to the projection system defined under userproject: system. For the CopDEM data in this example, the data are reprojected and tiled for the ease2n (EASE GRID 2.0 North) projection system.

Expanded virtual mosaic tiles

Spatial processing applying seamless tiles causes edge effects for kernel and regional based processes. To overcome that the Framework can use virtually expanded tiles as input in some processes. The process for creating virtually expanded tiles is MosaicAdjacentTiles.

Mosaic adjacent tiles

Framework process: MosaicAdjacentTiles

Json command file: 0190_MosaicAdjacentTiles_CopDEM-90m.json

The process MosaicAdjacentTiles creates the expanded virtual tiles. The post Mosaics for avoiding edge effects summarises the concept and processing.

Hydrological DEM corrections

Hydrological corrections of the Copernicus DEM is done in several steps:

  • Filling of single pits and peaks adjacent to streams
  • Adjacent tile mosaic the single cell corrected DEM
  • Filling of stream associated pits
  • Adjacent tile mosaic of the dual corrected DEM

Filling of single pits and peaks adjacent to streams

Framework process: MosaicAdjacentTiles

Json command file: 0230_GrassDemFillDirTiles_CopDEM_90m.json

The first hydrological correction is the filling up of all single cell pits and the flattening of all single cell peaks adjacent to streams. The Framework process to apply is GrassDemFillDirTiles, which is a wrapper for the GRASS GIS module r.fill.dir.

Adjacent tile mosaic the single cell corrected DEM

0191_MosaicAdjacentTiles_CopDEM-90m.json

Create expanded virtual tiles of the single cell corrected DEM tiles for use in the next step. See the post on Mosaics for avoiding edge effects for details.

Fill stream associated pits

0240_GrassDemHydroDemTiles_CopDEM-90m.json

The second hydrological correction is the filling up of larger pits in streams. The Framework process is GrassDemHydroDemTiles, which is a wrapper for the GRASS GIS addon module r.hydrodem.

Adjacent tile mosaic the dual corrected DEM

0192_MosaicAdjacentTiles_CopDEM-90m.json

Create expanded virtual tiles of the single cell corrected DEM tiles for use in the next step. See the post on Mosaics for avoiding edge effects for details.

River mouth and shoreline corrections

0210_GrassOnetoManyTiles-correct-shoreline_CopDEM-90m.json

The river mouth and shoreline correction processing fills up all negative pits adjacent to the sea. No special process is used, instead the generic Framework wrapper for GRASS GIS is used by defining a sequence of GRASS commands under the process GrassOnetoManyTiles.

GRASS script

The above process generates a script file that loops over all the tiles belonging to the defined region. An example is shown under the Hide/show button below, the commented lines (starting with “#”) have been added for explanation.

Kernel DEM derivates

With Karttur’s GeoImagine Framework, kernel (filter) indexes derived from DEM data can be calculated using GDAL, GRASS GIS, SAGA GIS or using internal processes.

GDAL DEM derivatives

0301_GdalDemTiles_CopDEM-90m.json

The Framework process GdalDemTiles is a wrapper for running the Geographic Data Abstraction Library (GDAL) DEM (gdaldem) program.

Framework DEM derivatives (NumpyDemTiles)

0303_NumpyDemTiles_CopDEM-90m.json

The Framework process for internal DEM derivates is focusing of calculating Terrain Position Index [TPI] landforms and can thus only retrieve the derivates for slope and TPI plus the landform categogization. The TPI extraction and thus the landform categorisation is, however, more advanced (thetorically more correct) compared to most other solutions.

GRASS DEM derivatives (GrassOnetoManyTiles)

GRASS GIS has more options for retrieving DEM derivates than either GDAL or the internal Framework process. GRASS modules for retrieving DEM derivatives include: r.slope.aspect, r.params.scale, r.tpi, r.tri and r.geomorphon. Retrieving DEM derivates using GRASS GIS though the Framework is done through the standard GRASS wrapper GrassOnetoManyTiles.

Hillslope DEM derivates

Hillslope indexes dereived from DEM data are calcualted using either GRASS GIS or SAGA GIS.

GRASS GIS hillslope DEM derivates

0311_GrassOnetoManyTiles_hillslope-derivates_copDEM-90m.json

Hillslope DEM indexes can be created by the general GRASS GIS wrapper GrassOnetoManyTiles.

GRASS script

The above process generates a script file that loops over all the tiles belonging to the defined region. An example is shown under the Hide/show button below, the commented lines (starting with “#”) have been added for explanation.

DEM basin extraction

0321_GrassOnetoManyTilesCopDEM-basin-extract-stage2_copDEM.json

One (additional) problem with extrating river basins is the definition of outlets points, a problem unrelated to the problem of hydrological flow routing. This particualr problem is, however, of relevance only if you intend to use the basin dat for hydrolocail modeling. If you only intend to use the basin map for extracting spatial statistics, you do not need an un-ambiguosly identifed outlet point. This part is for creating a distinct, single cell, outlet point for basins. It maske use of GRASS GIS and the Framework wrapper GrassOnetoManyTiles.

GRASS script

The above process generates a script file that loops over all the tiles belonging to the defined region. An example is shown under the Hide/show button below, the commented lines (starting with “#”) have been added for explanation.