Home

Awesome

pyDEM: A python digital elevation model analysis package

PyDEM is a package for topographic (terrain) analysis written in Python (with a bit of Cython). It takes in digital elevation model (DEM) rasters, and it outputs quantities like slope, aspect, upstream area, and topographic wetness index. PyDEM closely follows the approach taken by TauDEM to calculate these quantities. It is designed to be fast, easily extensible, and capable of handling very large datasets.

PyDEM can be used both from the command-line, or programmatically via its Python API. Examples of both usages are given below and in the examples directory. It can operate on individual elevation rasters, or on entire directories simultaneously. Most processing steps can be performed in parallel with any number of independent pyDEM processes. Note that pyDEM depends on TauDEM for certain steps (e.g., pitfilling) and it also makes extensive use of the GDAL library for working with geospatial rasters.

1. Installation

For installation notes, see install.md

2. Basic Usage

Examples and tests can be found in the pydem\examples directory.

2.1 Python Module Usage

2.1.1. Calculate quantities on a single elevation tile

Import the DEMProcessor class:

from pydem.dem_processing import DEMProcessor

Set the path to a pit filled elevation file in WGS84 coordinates:

filename_to_elevation_geotiff = 'test.tiff'

Instantiate an instance of the DEMProcessor class:

dem_proc = DEMProcessor(filename_to_elevation_geotiff)

The following three commands do not need to be called in order.

Calculate the aspect and slope magnitude:

mag, aspect = dem_proc.calc_slopes_directions()

Calculate the upstream contributing area:

uca = dem_proc.calc_uca()

Calculate the TWI:

twi = dem_proc.calc_twi()

2.1.2 Calculate TWI on a directory of elevation tiles

The ProcessManager class orchestrates multiple processes to compute TWI over multiple tiles in parallel on the same machine. Example usage is as follows:

  from pydem.processing_manager import ProcessManager
  elevation_source_path = r'/home/twi-users/elevation'
  manager = ProcessManager(
    n_workers=64, # Number of worker processes to use
    in_path=elevation_source_path,
    out_path='/home/twi-users/temporary_compute_storage/'  # Where to save intermediate data
    dem_proc_kwargs={},  # dictionary of values used to initialize DemProcessor
  )
  # Start the processing
  manager.process_twi()  # If this fails (e.g. machine goes down), can restart to pick up where it left off

You can also call individual parts of the processing:

manager.compute_grid()  # Figures out how elevation in folder are tiled
manager.process_elevation()  # Fills flats, handles pits, fixes artifacts in elevation data
manager.process_aspect_slope()
manager.process_uca() # Computes upstream contributing area (UCA) for individual tiles (embarassingly parallel)
manager.process_uca_edges() # Fixes upstream flow contribution accross tile edges -- iteratively
manager.process_twi()  # Call all the above functions internally, but skips any work already done

Finally, to export results to a single GeoTiff with overviews, use:

manager.save_non_overlap_data_geotiff(
  'float32',  # Numpy recognized filetype
  new_path=output_path,
  keys=['elev', 'uca', 'aspect', 'slope', 'twi'], # list can contain a subset of these
  overview_type='average')

Note: The name non_overlap_data comes from an implementation quirk. The temporary data is saved as a zarr file. This zarr file OVERLAPS data at the edges of tiles to make it easier to compute UCA across edges

2.1.3 DEMProcessor options

The following options are used by the DEMProcess object. They can be modified by setting the value before processing.

dem_proc = DEMProcessor(filename_to_elevation_geotiff)
dem_proc.fill_flats = False

Slopes & Directions

UCA

TWI

Other

2.2 Commandline Usage

When installing pydem using the provided setup.py file, the commandline utilities TWIDinf, AreaDinf, and DinfFlowDir are registered with the operating system.

TWIDinf :

usage: TWIDinf-script.py [-h] [--save-all]
                     Input_Pit_Filled_Elevation [Input_Number_of_Chunks]
                     [Output_D_Infinity_TWI]

Calculates a grid of topographic wetness index which is the log_e(uca / mag),
that is, the natural log of the ratio of contributing area per unit contour
length and the magnitude of the slope. Note, this function takes the elevation
as an input, and it calculates the slope, direction, and contributing area as
intermediate steps.

positional arguments:
  Input_Pit_Filled_Elevation
                    The input pit-filled elevation file in geotiff format.
  Input_Number_of_Chunks
                    The approximate number of chunks that the input file
                    will be divided into for processing (potentially on
                    multiple processors).
  Output_D_Infinity_TWI
                    Output filename for the topographic wetness index.
                    Default value = twi.tif .

optional arguments:
  -h, --help            show this help message and exit
  --save-all, --sa      If set, will save all intermediate files as well.

AreaDinf :

usage: AreaDinf-script.py [-h] [--save-all]
                      Input_Pit_Filled_Elevation [Input_Number_of_Chunks]
                      [Output_D_Infinity_Specific_Catchment_Area]

Calculates a grid of specific catchment area which is the contributing area
per unit contour length using the multiple flow direction D-infinity approach.
Note, this is different from the equivalent tauDEM function, in that it takes
the elevation (not the flow direction) as an input, and it calculates the
slope and direction as intermediate steps.

positional arguments:
  Input_Pit_Filled_Elevation
                    The input pit-filled elevation file in geotiff format.
  Input_Number_of_Chunks
                    The approximate number of chunks that the input file
                    will be divided into for processing (potentially on
                    multiple processors).
  Output_D_Infinity_Specific_Catchment_Area
                    Output filename for the flow direction. Default value = uca.tif .

optional arguments:
  -h, --help            show this help message and exit
  --save-all, --sa      If set, will save all intermediate files as well.

DinfFlowDir :

usage: DinfFlowDir-script.py [-h]
                         Input_Pit_Filled_Elevation
                         [Input_Number_of_Chunks]
                         [Output_D_Infinity_Flow_Direction]
                         [Output_D_Infinity_Slope]

Assigns a flow direction based on the D-infinity flow method using the
steepest slope of a triangular facet (Tarboton, 1997, "A New Method for the
Determination of Flow Directions and Contributing Areas in Grid Digital
Elevation Models," Water Resources Research, 33(2): 309-319).

positional arguments:
  Input_Pit_Filled_Elevation
                    The input pit-filled elevation file in geotiff format.
  Input_Number_of_Chunks
                    The approximate number of chunks that the input file
                    will be divided into for processing (potentially on
                    multiple processors).
  Output_D_Infinity_Flow_Direction
                    Output filename for the flow direction. Default value = ang.tif .
  Output_D_Infinity_Slope
                    Output filename for the flow direction. Default value = mag.tif .

optional arguments:
  -h, --help            show this help message and exit

3. Description of package Contents

4. References

Tarboton, D. G. (1997). A new method for the determination of flow directions and upslope areas in grid digital elevation models. Water resources research, 33(2), 309-319.

Ueckermann, Mattheus P., et al. (2015). "pyDEM: Global Digital Elevation Model Analysis." In K. Huff & J. Bergstra (Eds.), Scipy 2015: 14th Python in Science Conference. Paper presented at Austin, Texas, 6 - 12 July (pp. 117 - 124). http://conference.scipy.org/proceedings/scipy2015/mattheus_ueckermann.html

5. Attributions

pyDEM uses lib.pyx from the OpenPIV project for inpainting missing values in the final outputs.