pygeoprocessing.routing.routing module

Provides PyGeprocessing Routing functionality.

Unless otherwise specified, all internal computation of rasters are done in a float64 space. The only possible loss of precision could occur when an incoming DEM type is an int64 type and values in that dem exceed 2^52 but GDAL does not support int64 rasters so no precision loss is possible with a float64.

D8 flow direction conventions encode the flow direction as:

2 1
x 0
6 7

This is slightly different from how TauDEM encodes flow direction, which is as::: 4 3 2 5 x 1 6 7 8

To convert a TauDEM flow direction raster to a pygeoprocessing-compatible flow direction raster, the following raster_map call may be used:

taudem_flow_dir_path = 'taudem_d8_flow_dir.tif'
pygeoprocessing.raster_map(
    lambda d: d+1, [taudem_flow_dir_path],
    'pygeoprocessing_d8_flow_dir.tif')

pygeoprocessing.routing.routing.calculate_subwatershed_boundary()

Calculate a linestring boundary around all subwatersheds.

Subwatersheds start where the strahler_stream_vector has a junction starting at this highest upstream to lowest and ending at the outlet of a river.

Parameters:

d8_flow_dir_raster_path_band (tuple) – raster/path band for d8 flow dir raster
strahler_stream_vector_path (str) – path to stream segment vector
target_watershed_boundary_vector_path (str) –
path to created vector of linestring for watershed boundaries. Contains the fields:
- ”stream_id”: this is the stream ID from the strahler_stream_vector_path that corresponds to this subwatershed.
- ”terminated_early”: if set to 1 this watershed generation was terminated before it could be complete. This value should always be 0 unless something is wrong as a software bug or some degenerate case of data.
- ”outlet_x”, “outlet_y”: this is the x/y coordinate in raster space of the outlet of the watershed. It can be useful when determining other properties about the watershed when indexed with underlying raster data that created the streams in strahler_stream_vector_path.
max_steps_per_watershed (int) – maximum number of steps to take when defining a watershed boundary. Useful if the DEM is large and degenerate or some other user known condition to limit long large polygons. Defaults to 1000000.
outlet_at_confluence (bool) – If True the outlet of subwatersheds starts at the confluence of streams. If False (the default) subwatersheds will start one pixel up from the confluence.

Returns:

None.

pygeoprocessing.routing.routing.detect_lowest_drain_and_sink()

Find the lowest drain and sink pixel in the DEM.

This function is used to specify conditions to DEMs that are known to have one real sink/drain, but may have several numerical sink/drains by detecting both the lowest pixel that could drain the raster on an edge and the lowest internal pixel that might sink the whole raster.

Example

raster A contains the following

pixel at (3, 4) at 10m draining to a nodata pixel
pixel at (15, 19) at 11m draining to a nodata pixel
pixel at (19, 21) at 10m draining to a nodata pixel
pit pixel at (10, 15) at 5m surrounded by non-draining pixels
pit pixel at (25, 15) at 15m surrounded by non-draining pixels
pit pixel at (2, 125) at 5m surrounded by non-draining pixels

The result is two pixels indicating the first lowest edge and first lowest sink seen:

drain_pixel = (3, 4), 10 sink_pixel = (10, 15), 5

Parameters:

dem_raster_path_band (tuple) – a raster/path band tuple to detect sinks in.

Returns:

(drain_pixel, drain_height, sink_pixel, sink_height) -: two (x, y) tuples with corresponding heights, first list is for edge drains, the second is for pit sinks. The x/y coordinate is in raster coordinate space and _height is the height of the given pixels in edge and pit respectively.

pygeoprocessing.routing.routing.detect_outlets()

Create point vector indicating flow raster outlets.

If either D8 or MFD rasters have a flow direction to the edge of the raster or to a nodata flow direction pixel the originating pixel is considered an outlet.

Parameters:

flow_dir_raster_path_band (tuple) – raster path/band tuple indicating D8 or MFD flow direction created by routing.flow_dir_d8 or routing.flow_dir_mfd.
flow_dir_type (str) – one of ‘d8’ or ‘mfd’ to indicate the flow_dir_raster_path_band is either a D8 or MFD flow direction raster.
target_outlet_vector_path (str) –
path to a vector that is created by this call that will be in the same projection units as the raster and have a point feature in the center of each pixel that is a raster outlet. Additional fields include:
- ”i” - the column raster coordinate where the outlet exists
- ”j” - the row raster coordinate where the outlet exists
- ”ID” - unique identification for the outlet.

Returns:

None.

pygeoprocessing.routing.routing.distance_to_channel_d8()

Calculate distance to channel with D8 flow.

Parameters:

flow_dir_d8_raster_path_band (tuple) –
a path/band index tuple indicating the raster that defines the D8 flow direction raster for this call. The pixel values are integers that correspond to outflow in the following configuration:
```
3 2 1
4 x 0
5 6 7
```
channel_raster_path_band (tuple) – a path/band tuple of the same dimensions and projection as flow_dir_d8_raster_path_band[0] that indicates where the channels in the problem space lie. A channel is indicated if the value of the pixel is 1. Other values are ignored.
target_distance_to_channel_raster_path (str) – path to a raster created by this call that has per-pixel distances from a given pixel to the nearest downhill channel.
weight_raster_path_band (tuple) – optional path and band number to a raster that will be used as the per-pixel flow distance weight. If None, 1 is the default distance between neighboring pixels. This raster must be the same dimensions as flow_dir_mfd_raster_path_band.
raster_driver_creation_tuple (tuple) – a tuple containing a GDAL driver name string as the first element and a GDAL creation options tuple/list as the second. Defaults to a GTiff driver tuple defined at geoprocessing.DEFAULT_GTIFF_CREATION_TUPLE_OPTIONS.

Returns:

None.

pygeoprocessing.routing.routing.distance_to_channel_mfd()

Calculate distance to channel with multiple flow direction.

MFD distance to channel for a pixel i is the average distance that water travels from i until it reaches a stream, defined as:

       ⎧ 0,                                       if i is a stream pixel
       ⎪
       ⎪ undefined,   if i is not a stream and there are 0 elements in N
d_i =  ⎨
       ⎪
       ⎪   ⎲ ((d_n + l(i, n)) * p(i, n),                       otherwise
       ⎪   ⎳
       ⎩  n ∈ N

 where
 - N is the set of immediate downstream neighbors of i that drain to
   a stream on the map
 - p(i, n) is the proportion of water on pixel i that flows onto n
 - l(i, n) is the distance between i and n. This is 1 for lateral
   neighbors, and √2 for diagonal neighbors.

Distance is measured in units of pixel lengths. Pixels must be square, and both input rasters must have the same spatial reference, geotransform, and extent, so that their pixels overlap exactly.

If a none of a pixel’s flow drains to any stream on the map, d_i is undefined, and that pixel has nodata in the output.

If only some of a pixel’s flow reaches a stream on the map, d_i is defined in terms of only that portion of flow that reaches a stream. It is the average distance traveled by the water that does reach a stream.

The algorithm begins from an arbitrary pixel and pushes downstream neighbors to the stack depth-first, until reaching a stream (the first case in the function for d_i above) or a raster edge (the second case). Non-stream pixels are only calculated once all their downstream neighbors have been calculated (the third case).

Parameters:

flow_dir_mfd_raster_path_band (tuple) – a path/band index tuple indicating the raster that defines the mfd flow accumulation raster for this call. This raster should be generated by a call to pygeoprocessing.routing.flow_dir_mfd.
channel_raster_path_band (tuple) – a path/band tuple of the same dimensions and projection as flow_dir_mfd_raster_path_band[0] that indicates where the channels in the problem space lie. A channel is indicated if the value of the pixel is 1. Other values are ignored.
target_distance_to_channel_raster_path (str) – path to a raster created by this call that has per-pixel distances from a given pixel to the nearest downhill channel.
weight_raster_path_band (tuple) – optional path and band number to a raster that will be used as the per-pixel flow distance weight. If None, 1 is the default distance between neighboring pixels. This raster must be the same dimensions as flow_dir_mfd_raster_path_band.
raster_driver_creation_tuple (tuple) – a tuple containing a GDAL driver name string as the first element and a GDAL creation options tuple/list as the second. Defaults to a GTiff driver tuple defined at geoprocessing.DEFAULT_GTIFF_CREATION_TUPLE_OPTIONS.

Returns:

None.

pygeoprocessing.routing.routing.extract_strahler_streams_d8()

Extract Strahler order stream geometry from flow accumulation.

Creates a Strahler ordered stream vector containing line segments representing each separate stream fragment. The final vector contains at least the fields:

“order” (int): an integer representing the stream order

“river_id” (int): unique ID used by all stream segments that
connect to the same outlet.

“drop_distance” (float): this is the drop distance in DEM units
from the upstream to downstream component of this stream segment.

“outlet” (int): 1 if this segment is an outlet, 0 if not.

“us_fa” (int): flow accumulation value at the upstream end of
the stream segment.

“ds_fa” (int): flow accumulation value at the downstream end of
the stream segment

“thresh_fa” (int): the final threshold flow accumulation value
used to determine the river segments.

“upstream_d8_dir” (int): a bookkeeping parameter from stream
calculations that is left in due to the overhead of deleting a field.

“ds_x” (int): the raster x coordinate for the outlet.

“ds_y” (int): the raster y coordinate for the outlet.

“ds_x_1” (int): the x raster space coordinate that is 1 pixel
upstream from the outlet.

“ds_y_1” (int): the y raster space coordinate that is 1 pixel
upstream from the outlet.

“us_x” (int): the raster x coordinate for the upstream inlet.

“us_y” (int): the raster y coordinate for the upstream inlet.

Parameters:

flow_dir_d8_raster_path_band (tuple) – a path/band representing the D8 flow direction raster.
flow_accum_raster_path_band (tuple) – a path/band representing the D8 flow accumulation raster represented by flow_dir_d8_raster_path_band.
dem_raster_path_band (tuple) – a path/band representing the DEM used to derive flow dir.
target_stream_vector_path (tuple) – a single layer line vector created by this function representing the stream segments extracted from the above arguments. Contains the fields “order” and “parent” as described above.
min_flow_accum_threshold (int) – minimum number of upstream pixels required to create a stream. If autotune_flow_accumulation is True, then the final value may be adjusted based on significant differences in 1st and 2nd order streams.
river_order (int) – what stream order to define as a river in terms of automatically determining flow accumulation threshold for that stream collection.
min_p_val (float) – minimum p_value test for significance
autotune_flow_accumulation (bool) – If true, uses a t-test to test for significant distances in order 1 and order 2 streams. If it is significant the flow accumulation parameter is adjusted upwards until the drop distances are insignificant.
osr_axis_mapping_strategy (int) – OSR axis mapping strategy for SpatialReference objects. Defaults to geoprocessing.DEFAULT_OSR_AXIS_MAPPING_STRATEGY. This parameter should not be changed unless you know what you are doing.

Returns:

None.

pygeoprocessing.routing.routing.extract_streams_mfd()

Classify a stream raster from MFD flow accumulation.

This function classifies pixels as streams that have a flow accumulation value >= flow_threshold and can trace further upstream with a fuzzy propotion if trace_threshold_proportion is set < 1.0

Parameters:

flow_accum_raster_path_band (tuple) – a string/integer tuple indicating the flow accumulation raster to use as a basis for thresholding a stream. Values in this raster that are >= flow_threshold will be classified as streams. This raster should be derived from dem_raster_path_band using pygeoprocessing.routing.flow_accumulation_mfd.
flow_dir_mfd_path_band (str) – path to multiple flow direction raster, required to join divergent streams.
flow_threshold (float) – the value in flow_accum_raster_path_band to indicate where a stream exists.
target_stream_raster_path (str) – path to the target stream raster. This raster will be the same dimensions and projection as dem_raster_path_band and will contain 1s where a stream is defined, 0 where the flow accumulation layer is defined but no stream exists, and nodata otherwise.
trace_threshold_proportion (float) – this value indicates what proportion of the flow_threshold is enough to classify a pixel as a stream after the stream has been traced from a flow_threshold drain. Setting this value < 1.0 is useful for classifying streams in regions that have highly divergent flow directions.
raster_driver_creation_tuple (tuple) – a tuple containing a GDAL driver name string as the first element and a GDAL creation options tuple/list as the second. Defaults to a GTiff driver tuple defined at geoprocessing.DEFAULT_GTIFF_CREATION_TUPLE_OPTIONS.

Returns:

None.

pygeoprocessing.routing.routing.fill_pits()

Fill the pits in a DEM.

This function defines pits as hydrologically connected regions that do not drain to the edge of the raster or a nodata pixel. After the call pits are filled to the height of the lowest pour point.

Parameters:

dem_raster_path_band (tuple) – a path, band number tuple indicating the DEM calculate flow direction.
target_filled_dem_raster_path (str) – path the pit filled dem, that’s created by a call to this function. It is functionally a single band copy of dem_raster_path_band with the pit pixels raised to the pour point. For runtime efficiency, this raster is tiled and its blocksize is set to (1<<BLOCK_BITS, 1<<BLOCK_BITS) even if dem_raster_path_band[0] was not tiled or a different block size.
working_dir (str) – If not None, indicates where temporary files should be created during this run. If this directory doesn’t exist it is created by this call. If None, a temporary directory is created by tempdir.mkdtemp which is removed after the function call completes successfully.
max_pixel_fill_count (int) – maximum number of pixels to fill a pit before leaving as a depression. Useful if there are natural large depressions. Value of -1 fills the raster with no search limit.
single_outlet_tuple (tuple) – If not None, this is an x/y tuple in raster coordinates indicating the only pixel that can be considered a drain. If None then any pixel that would drain to the edge of the raster or a nodata hole will be considered a drain.
raster_driver_creation_tuple (tuple) – a tuple containing a GDAL driver name string as the first element and a GDAL creation options tuple/list as the second. Defaults to a GTiff driver tuple defined at geoprocessing.DEFAULT_GTIFF_CREATION_TUPLE_OPTIONS.

Returns:

None.

pygeoprocessing.routing.routing.flow_accumulation_d8()

D8 flow accumulation.

Parameters:

flow_dir_raster_path_band (tuple) –
a path, band number tuple for a flow accumulation raster whose pixels indicate the flow out of a pixel in one of 8 directions in the following configuration:
```
3 2 1
4 x 0
5 6 7
```
target_flow_accum_raster_path (str) – path to flow accumulation raster created by this call. After this call, the value of each pixel will be 1 plus the number of upstream pixels that drain to that pixel. Note the target type of this raster is a 64 bit float so there is minimal risk of overflow and the possibility of handling a float dtype in weight_raster_path_band.
weight_raster_path_band (tuple) – optional path and band number to a raster that will be used as the per-pixel flow accumulation weight. If None, 1 is the default flow accumulation weight. This raster must be the same dimensions as flow_dir_mfd_raster_path_band.
raster_driver_creation_tuple (tuple) – a tuple containing a GDAL driver name string as the first element and a GDAL creation options tuple/list as the second. Defaults to a GTiff driver tuple defined at geoprocessing.DEFAULT_GTIFF_CREATION_TUPLE_OPTIONS.

Returns:

None.

pygeoprocessing.routing.routing.flow_accumulation_mfd()

Multiple flow direction accumulation.

Parameters:

flow_dir_mfd_raster_path_band (tuple) – a path, band number tuple for a multiple flow direction raster generated from a call to flow_dir_mfd. The format of this raster is described in the docstring of that function.
target_flow_accum_raster_path (str) – a path to a raster created by a call to this function that is the same dimensions and projection as flow_dir_mfd_raster_path_band[0]. The value in each pixel is 1 plus the proportional contribution of all upstream pixels that flow into it. The proportion is determined as the value of the upstream flow dir pixel in the downslope direction pointing to the current pixel divided by the sum of all the flow weights exiting that pixel. Note the target type of this raster is a 64 bit float so there is minimal risk of overflow and the possibility of handling a float dtype in weight_raster_path_band.
weight_raster_path_band (tuple) – optional path and band number to a raster that will be used as the per-pixel flow accumulation weight. If None, 1 is the default flow accumulation weight. This raster must be the same dimensions as flow_dir_mfd_raster_path_band. If a weight nodata pixel is encountered it will be treated as a weight value of 0.
raster_driver_creation_tuple (tuple) – a tuple containing a GDAL driver name string as the first element and a GDAL creation options tuple/list as the second. Defaults to a GTiff driver tuple defined at geoprocessing.DEFAULT_GTIFF_CREATION_TUPLE_OPTIONS.

Returns:

None.

pygeoprocessing.routing.routing.flow_dir_d8()

D8 flow direction.

Parameters:

dem_raster_path_band (tuple) – a path, band number tuple indicating the DEM calculate flow direction. This DEM must not have hydrological pits or else the target flow direction is undefined.
target_flow_dir_path (str) –
path to a byte raster created by this call of same dimensions as dem_raster_path_band that has a value indicating the direction of downhill flow. Values are defined as pointing to one of the eight neighbors with the following convention:
```
3 2 1
4 x 0
5 6 7
```
working_dir (str) – If not None, indicates where temporary files should be created during this run. If this directory doesn’t exist it is created by this call.
raster_driver_creation_tuple (tuple) – a tuple containing a GDAL driver name string as the first element and a GDAL creation options tuple/list as the second. Defaults to a GTiff driver tuple defined at geoprocessing.DEFAULT_GTIFF_CREATION_TUPLE_OPTIONS.

Returns:

None.

pygeoprocessing.routing.routing.flow_dir_mfd()

Multiple flow direction.

Parameters:

dem_raster_path_band (tuple) – a path, band number tuple indicating the DEM calculate flow direction. This DEM must not have hydrological pits or else the target flow direction will be undefined.
target_flow_dir_path (str) –
path to a raster created by this call of a 32 bit int raster of the same dimensions and projections as dem_raster_path_band[0]. The value of the pixel indicates the proportion of flow from that pixel to its neighbors given these indexes:
```
3 2 1
4 x 0
5 6 7
```
The pixel value is formatted as 8 separate 4 bit integers compressed into a 32 bit int. To extract the proportion of flow from a particular direction given the pixel value ‘x’ one can shift and mask as follows 0xF & (x >> (4*dir)), where dir is one of the 8 directions indicated above.
working_dir (str) – If not None, indicates where temporary files should be created during this run. If this directory doesn’t exist it is created by this call.
raster_driver_creation_tuple (tuple) – a tuple containing a GDAL driver name string as the first element and a GDAL creation options tuple/list as the second. Defaults to a GTiff driver tuple defined at geoprocessing.DEFAULT_GTIFF_CREATION_TUPLE_OPTIONS.

Returns:

None.