PyDEM is a package for topographic (terrain) analysis written in Python 2.7 (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.
For installation notes, see install.md
Examples and tests can be found in the
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
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()
from pydem.processing_manager import ProcessManager
Set the path to the directory of pit-filled WGS84 geotiff elevation files. Important: The files in this directory needs to follow the naming convention used. See
examples.process_manager_directory.py for an example of using
utils.rename_files to rename elevation files.
elevation_source_path = r'C:\test_directory'
Set the target location (where files will be saved):
# This is the default location if unspecified. elevation_target_path = r'C:\test_directory\processed_data'
Instantiate an instance of the
pm = ProcessManager(elevation_source_path, elevation_target_path)
Process the elevation files in that directory:
Topographic wetness index calculated for each elevation file in the supplied directory
will be located in
C:\test_directory\processed_data\twi. These TWI files will not have
edge effects on edges interior to the data set.
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
chunk_size_slp_dir: Chunk size for slopes_directions calculation. Default
chunk_overlap_slp_dir: Overlap to use for resolving slopes and directions at chunk edges. Default
fill_flats: Fill/interpolate the elevation for flat regions before calculating slopes and directions. The direction cannot be calculated in regions where the slope is 0 because the nominal elevation is all the same. This can happen in very gradual terrain or in lake and river beds, particularly when the input elevation is composed of integers. When
True, the elevation is interpolated in those regions so that a reasonable slope and direction can be calculated. Default
fill_flats_below_sea: Interpolate the elevation for flat regions that are below sea level. Water will never flow out of these "pits", so in many cases you can ignore these regions and achieve faster processing times. Default
fill_flats_source_tol: When filling flats, the algorithm finds adjacent "source" pixels and "drain" pixels for each flat region and interpolates the elevation using these data points. This sets the tolerance for the elevation of source pixels above the flat region (i.e. shallow sources are used as sources but not steep cliffs). Default
fill_flats_peaks: Interpolate the elevation for flat regions that are "peaks" (local maxima). These regions have a higher elevation than all adjacent pixels, so there are no "source" pixels to use for interpolation. When
True, a single pixel is selected approximately in the center of the flat region as the "peak"/"source". Default
fill_flats_pits: Interpolate the elevation for flat regions that are "pits" (local minima). These regions have a lower elevation than all adjacent pixels, so there are no "drain" pixels to use for interpolation. When
True, a single pixel is selected approximately in the center of the flat region as the "pit"/"drain". Default
resolve_edges: Ensure edge UCA is continuous across chunks. Default
chunk_size_uca: Chunk size for uca calculation. Default
chunk_overlap_uca: Overlap to use for resolving uca at chunk edges. Default
drain_pits: Drain from "pits" to nearby but non-adjacent pixels. Pits have no lower adjacent pixels to drain to directly. Note that with
fill_flats_pits off, this setting will still drain each pixel in large flat regions, but it may be slower and produces less reasonable results. Default
drain_pits_max_iter: Maximum number of iterations to look for drain pixels for pits. Generally, "nearby drains" for a pit/flat region are found by expanding the region upward/outward iteratively. Default
drain_pits_max_dist: Maximum distance in coordnate-space to (non-adjacent) drains for pits. Pits that are too far from another pixel with a lower elevation will not drain. Default
drain_pits_max_dist_XY: Maximum distance in real-space to (non-adjacent) drains for pits. Pits that are too far from another pixel with a lower elevation will not drain. This filter is applied after
drain_pits_max_dist; if the X and Y resolution are similar, this filter is generally unnecessary. Default
drain_flats: [Deprecated, replaced by
drain_pits] Drains flat regions and pits by draining all pixels in the region to an arbitrary pixel in the region and then draining that pixel to the border of the flat region. Ignored if
apply_uca_limit_edges: Mark edges as completed if the maximum UCA is reached when resolving drainage across edges. Default
False. If True, it may speed up large calculations.
apply_twi_limits: When calculating TWI, limit TWI to max value. Default
apply_twi_limits_on_uca: When calculating TWI, limit UCA to max value. Default
Import and Instantiate a
from pydem.processing_manager import ProcessManager
pm = ProcessManager(r'C:\test_directory')
Define a custom command. For example, to calculate hill-shading using gdal:
def command(source_elevation_file, target_file_save_location): cmd = ['gdaldem', 'hillshade', '-s', '111120', '-compute_edges', '-co', 'BIGTIFF=YES', '-of', 'GTiff', '-co', 'compress=lzw', '-co', 'TILED=YES', esfile, fn] print '<'*8, ' '.join(cmd), '>'*8 status = subprocess.call(cmd) return status
Process directory of files using custom command: pm.process_command(command, save_name='hillshade_files')
The results of this operation will be found in
When installing pydem using the provided setup.py file, the commandline utilities
DinfFlowDir are registered with the operating system.
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.
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.
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
commandline_utils.py: Contains the functions that wrap the python modules into command line utilities.
dem_processing.py: Contains the main algorithms.
processing_manager.py: Implements a class that manages the calculation of TWI for a directory of files.
test_pydem.py: A few helper utilities that create analytic test-cases used to develop/test pyDEM.
utils.py: A few helper utility functions.
gdalfunctions for reading and writing geotiff files.
cyfuncs: Directory containing cythonized versions of python functions in
cyfuncs.cyutils.pyx: Computationally efficient implementations of algorithms used to calculate upstream contributing area.
examples: Directory containing a few examples, along with an end-to-end test of the cross-tile calculations.
examples.compare_tile_to_chunk.py: Compares the calculation of the upstream contributing area over a full tile compared to multiple chunks in a file. This tests that the upstream contributing area calculation correctly drains across tile edges.
examples.compare_to_taudem.py: This compares the calculation of magnitude and aspect to taudem's algorithms. This validates that the algorithms are correctly implemented from Tarboton (1997), and also shows the differences when taking the change in coordinates into account.
examples.cross-tile_process_manager_test.py: End-to-end test to make sure that the cross-tile calculations are correctly performed.
examples.process_manager_directory.py: This shows how to use the
ProcessingManagerto calculate all of the elevation files within a directory.
reader: Directory containing python code used to deal with opening and closing geotiff files. Essentially wraps
gdalto provide simpler usage.
reader.gdal_reader.py: Contains the `GDALReader class used to read and write geotiff files.
reader.inpaint.pyx: Cython function used to fill no-data values in geotiffs.
reader.my_types.py: Defines classes used to deal with different grid coordinate systems.
taudem: Directory containing a copy of taudem for convenience.
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