Qiusheng Wu
Qiusheng Wu
Assistant Professor of Geography at Binghamton University (SUNY)

A tutorial for the richdem Python package

This tutorial demonstrates the usage of the richdem Python package for depression filling.

  • HTML version: https://gishub.org/richdem-tutorial
  • GitHub repo: https://github.com/giswqs/richdem-binder
  • richdem: https://github.com/r-barnes/richdem

Launch this tutorial in an interactive cloud environment:

Binder

Binder

gif

Table of Content

Installation

The richdem Python package can be installed using the following command:

pip install richdem

If you have installed richdem Python package before and want to upgrade to the latest version, you can use the following command:

pip install richdem -U

Getting data

This section demonstrates two ways to get data into Binder so that you can test the richdem Python package in the cloud using your own data. We will use the pygis Python package to download data into Binder

Getting data from direct URLs

If you have data hosted on your own HTTP server or GitHub, you should be able to get direct URLs. With a direct URL, users can automatically download the data when the URL is clicked. For example https://github.com/giswqs/lidar/raw/master/examples/lidar-dem.zip

Import the following Python libraries and start getting data from direct URLs.

import os
import pygis

Create a folder named data and set it as the working directory.

root_dir = os.getcwd()
work_dir = os.path.join(root_dir, 'data')
if not os.path.exists(work_dir):
    os.mkdir(work_dir)
print("Working directory: {}".format(os.path.realpath(work_dir)))
Working directory: /home/jovyan/examples/data

Replace the following URL with your own direct URL hosting the data you would like to use.

url = "https://github.com/giswqs/lidar/raw/master/examples/lidar-dem.zip"

Download data the from the above URL and unzip the file if needed.

pygis.download_from_url(url, out_dir=work_dir)
Downloading lidar-dem.zip ...
Downloading done.
Unzipping lidar-dem.zip ...
Unzipping done.
Data directory: /home/jovyan/examples/data/lidar-dem

You have successfully downloaded data to Binder. Therefore, you can skip to Using richdem and start testing richdem with your own data.

Getting data from Google Drive

Alternatively, you can upload data to Google Drive and then share files publicly from Google Drive. Once the file is shared publicly, you should be able to get a shareable URL. For example, https://drive.google.com/file/d/1c6v-ep5-klb2J32Nuu1rSyqAc8kEtmdh.

Replace the following URL with your own shareable URL from Google Drive.

gfile_url = 'https://drive.google.com/file/d/1c6v-ep5-klb2J32Nuu1rSyqAc8kEtmdh'

Download the shared file from Google Drive.

pygis.download_from_gdrive(gfile_url, file_name='lidar-dem.zip', out_dir=work_dir)
Google Drive file id: 1c6v-ep5-klb2J32Nuu1rSyqAc8kEtmdh
Downloading 1c6v-ep5-klb2J32Nuu1rSyqAc8kEtmdh into /home/jovyan/examples/data/lidar-dem.zip... Done.
Unzipping...Done.

You have successfully downloaded data from Google Drive to Binder. You can now continue to Using richdem and start testing richdem with your own data.

Using richdem

Checking GDAL

The LoadGDAL function in richdem requires GDAL, so let’s first make sure GDAL has been installed in the system.

import os
print(os.popen('gdalinfo --version').read().rstrip())
GDAL 2.3.2, released 2018/09/21

Loading data

More information about loading data in richdem can be found here.

import richdem as rd
beau = rd.LoadGDAL('data/lidar-dem/dem.tif')

Depression filling

Depressions, otherwise known as pits, are areas of a landscape wherein flow ultimately terminates without reaching an ocean or the edge of a digital elevation model.

More information about depression-filling can be found on richdem.readthedocs.io

Original DEM

For reference, the original DEM appears as follows:

%matplotlib inline
beaufig = rd.rdShow(beau, ignore_colours=[0], axes=False, cmap='jet', figsize=(8,5.5))

png

Complete Filling

Depression-filling is often used to fill in all the depressions in a DEM to the level of their lowest outlet or spill-point.

The result looks as follows:

beau_filled    = rd.FillDepressions(beau, in_place=False)
beaufig_filled = rd.rdShow(beau_filled, ignore_colours=[0], axes=False, cmap='jet', vmin=beaufig['vmin'], vmax=beaufig['vmax'], figsize=(8,5.5))

png

We can visualize the difference between the two like so:

beau_diff    = beau_filled - beau
beaufig_diff = rd.rdShow(beau_diff, ignore_colours=[0], axes=False, cmap='jet', figsize=(8,5.5))

png

Epsilon Filling

A downside of complete filling is that it replaces depressions with a perfectly flat region with no local gradients. One way to deal with this is to ensure that every cell in the region is raised some small amount, ε, above cells which are closer to a depression’s spill point.

This must be done carefully. In floating-point DEMs, the value ε is non-constant and must be chosen using the !std::nextafter function. If a depression is too large, the imposed gradient may result in the interior of the depression being raised above the surrounding landscape. Using double instead of float reduces the potential for problems at a cost of twice the space used. If a problem does arise, RichDEM provides a warning.

We can visualize the difference between the epsilon-filled DEM and the original DEM like so:

beau_epsilon         = rd.FillDepressions(beau, epsilon=True, in_place=False)
beau_eps_diff        = beau_epsilon - beau
beaufig_eps_diff     = rd.rdShow(beau_eps_diff, ignore_colours=[0], axes=False, cmap='jet', figsize=(8,5.5))

png

We can visualize the difference between the epsilon-filled DEM and the completely-filled DEM as follows. Note that elevation increases with distance from the depression’s outlet: this is the effect of the epsilon.

beau_diffeps_diff    = beau_epsilon - beau_filled
beaufig_diffeps_diff = rd.rdShow(beau_diffeps_diff, ignore_colours=[0], axes=False, cmap='jet', figsize=(8,5.5))

png

comments powered by Disqus