Use D8 stream extraction

HISTORY.rst

@@ -45,11 +45,14 @@ Unreleased Changes
       interface to InVEST. (`#755 <https://github.com/natcap/invest/issues/755>`_)
     * The classic InVEST user-interface has been removed in favor of the Workbench.
     * Replace the ``ARGS_SPEC`` with ``MODEL_SPEC`` which describes all model
-      outputs as well as inputs in a structured format 
+      outputs as well as inputs in a structured format
       (`#596 <https://github.com/natcap/invest/issues/596>`_)
 * Workbench
     * Added tooltips to the model tabs so that they can be identified even when
       several tabs are open (`#1071 <https://github.com/natcap/invest/issues/1088>`_)
+* DelineateIt
+    * DelineateIt now uses ``pygeoprocessing.routing.extract_streams_d8`` for D8
+      stream thresholding. https://github.com/natcap/invest/issues/1143
     * Represent boolean inputs with a toggle switch rather than radio buttons.
 * Habitat Quality
     * The model now uses an euclidean distance implementation for decaying
@@ -69,6 +72,12 @@ Unreleased Changes
       trailing whitespace is also now removed from all string fields in the
       criteria table, which should also help reduce the chance of errors.
       https://github.com/natcap/invest/issues/1191
+* GLOBIO
+    * Deprecated the GLOBIO model
+      (`#1131 <https://github.com/natcap/invest/issues/1131>`_)
+* RouteDEM
+    * RouteDEM now uses ``pygeoprocessing.routing.extract_streams_d8`` for D8
+      stream thresholding. https://github.com/natcap/invest/issues/1143
 * Scenic Quality
     * Any points over nodata (and therefore excluded from the viewshed
       analysis) will now correctly have their FID reported in the logging.
@@ -109,9 +118,6 @@ Unreleased Changes
     * Fixed a bug where uncaught exceptions in the React tree would result in
       a blank browser window.
       (`#1119 <https://github.com/natcap/invest/issues/1119>`_)
-* GLOBIO
-    * Deprecated the GLOBIO model
-      (`#1131 <https://github.com/natcap/invest/issues/1131>`_)
 * Habitat Quality
     * All spatial inputs including the access vector and threat rasters are
       now reprojected to the ``lulc_cur_path`` raster. This fixes a bug where

requirements.txt

@@ -17,7 +17,7 @@ numpy>=1.11.0,!=1.16.0
 Rtree>=0.8.2,!=0.9.1
 shapely>=1.7.1,<1.8.2  # https://github.com/shapely/shapely/issues/1385
 scipy>=1.9.0
-pygeoprocessing>=2.3.5  # pip-only
+pygeoprocessing>=2.4.0  # pip-only
 taskgraph[niced_processes]>=0.11.0  # pip-only
 psutil>=5.6.6
 chardet>=3.0.4

src/natcap/invest/delineateit/delineateit.py

@@ -290,23 +290,17 @@ def execute(args):
         snap_distance = int(args['snap_distance'])
         flow_threshold = int(args['flow_threshold'])
 
-        out_nodata = 255
-        flow_accumulation_task.join()  # wait so we can read the nodata value
-        flow_accumulation_nodata = pygeoprocessing.get_raster_info(
-            file_registry['flow_accumulation'])['nodata']
         streams_task = graph.add_task(
-            pygeoprocessing.raster_calculator,
-            args=([(file_registry['flow_accumulation'], 1),
-                   (flow_accumulation_nodata, 'raw'),
-                   (out_nodata, 'raw'),
-                   (flow_threshold, 'raw')],
-                  _threshold_streams,
-                  file_registry['streams'],
-                  gdal.GDT_Byte,
-                  out_nodata),
+            pygeoprocessing.routing.extract_streams_d8,
+            kwargs={
+                'flow_accum_raster_path_band':
+                    (file_registry['flow_accumulation'], 1),
+                'flow_threshold': flow_threshold,
+                'target_stream_raster_path': file_registry['streams'],
+            },
             target_path_list=[file_registry['streams']],
             dependent_task_list=[flow_accumulation_task],
-            task_name='threshold_streams')
+            task_name='extract streams')
 
         snapped_outflow_points_task = graph.add_task(
             snap_points_to_nearest_stream,

src/natcap/invest/routedem.py

@@ -119,7 +119,7 @@
     'D8': {
         'flow_accumulation': pygeoprocessing.routing.flow_accumulation_d8,
         'flow_direction': pygeoprocessing.routing.flow_dir_d8,
-        'threshold_flow': None,  # Defined in source code as _threshold_flow
+        'threshold_flow': pygeoprocessing.routing.extract_streams_d8,
         'distance_to_channel': pygeoprocessing.routing.distance_to_channel_d8,
     },
     'MFD': {
@@ -131,36 +131,6 @@
 }
 
 
-def _threshold_flow(flow_accum_pixels, threshold, in_nodata, out_nodata):
-    """Raster_calculator local_op to threshold D8 stream flow.
-
-    Args:
-        flow_accum_pixels (numpy.ndarray): Array representing the number of
-            pixels upslope of a given pixel.
-        threshold (int or float): The threshold above which we have a stream.
-        in_nodata (int or float): The nodata value of the flow accumulation
-            raster.
-        out_nodata (int): The nodata value of the target stream mask raster.
-
-    Returns:
-        A numpy.ndarray (dtype is numpy.uint8) with pixel values of 1 where
-        flow accumulation > threshold, 0 where flow accumulation < threshold
-        and out_nodata where flow accumulation is equal to in_nodata.
-
-    """
-    out_matrix = numpy.empty(flow_accum_pixels.shape, dtype=numpy.uint8)
-    out_matrix[:] = out_nodata
-    stream_mask = (flow_accum_pixels > threshold)
-
-    valid_mask = slice(None)
-    if in_nodata is not None:
-        valid_mask = ~utils.array_equals_nodata(flow_accum_pixels, in_nodata)
-
-    out_matrix[valid_mask & stream_mask] = 1
-    out_matrix[valid_mask & ~stream_mask] = 0
-    return out_matrix
-
-
 def execute(args):
     """RouteDEM: Hydrological routing.
 
@@ -292,8 +262,7 @@ def execute(args):
             flow_accum_task = graph.add_task(
                 routing_funcs['flow_accumulation'],
                 args=((flow_dir_path, 1),
-                      flow_accumulation_path
-                ),
+                      flow_accumulation_path),
                 target_path_list=[flow_accumulation_path],
                 task_name='flow_accumulation_%s' % algorithm,
                 dependent_task_list=[flow_direction_task])
@@ -302,24 +271,20 @@ def execute(args):
                     bool(args['calculate_stream_threshold'])):
                 stream_mask_path = os.path.join(
                         args['workspace_dir'],
-                    _STREAM_MASK_FILE_PATTERN % file_suffix)
+                        _STREAM_MASK_FILE_PATTERN % file_suffix)
                 if algorithm == 'D8':
-                    flow_accum_task.join()
-                    flow_accum_info = pygeoprocessing.get_raster_info(
-                        flow_accumulation_path)
                     stream_threshold_task = graph.add_task(
-                        pygeoprocessing.raster_calculator,
-                        args=(((flow_accumulation_path, 1),
-                               (float(args['threshold_flow_accumulation']), 'raw'),
-                               (flow_accum_info['nodata'][0], 'raw'),
-                               (255, 'raw')),
-                              _threshold_flow,
-                              stream_mask_path,
-                              gdal.GDT_Byte,
-                              255),
+                        pygeoprocessing.routing.extract_streams_d8,
+                        kwargs={
+                            'flow_accum_raster_path_band':
+                                (flow_accumulation_path, 1),
+                            'flow_threshold': float(
+                                args['threshold_flow_accumulation']),
+                            'target_stream_raster_path': stream_mask_path,
+                        },
                         target_path_list=[stream_mask_path],
-                        task_name='stream_thresholding_D8',
-                        dependent_task_list=[flow_accum_task])
+                        dependent_task_list=[flow_accum_task],
+                        task_name='stream_thresholding_d8')
                 else:  # MFD
                     stream_threshold_task = graph.add_task(
                         routing_funcs['threshold_flow'],
@@ -344,6 +309,7 @@ def execute(args):
                         target_path_list=[distance_path],
                         task_name='downslope_distance_%s' % algorithm,
                         dependent_task_list=[stream_threshold_task])
+    graph.close()
     graph.join()