Repository: natcap/invest
Pull request #318 (2020-09-22 23:16:55)

Pour point detection function
- New function `delineateit.detect_pour_points` that takes in a flow direction raster (path, band) tuple and an output path, detects pour points from the flow direction raster, and saves them as a point vector to the output path.

- New Cython module for delineateit that does pour point detection by block. We agreed that it's not necessary right now to put all of the `detect_pour_points` functionality into the Cython module, because it's easier to read and test when most of it stays in the python module, and we are happy with the performance.

- Added UI checkbox option 'Detect pour points'. If this is checked:
    - the 'outlet features vector' input is grayed out. Detected pour points will be used instead of outlet geometries.
    - the 'skip invalid geometries' option is grayed out (detected pour point geometries should always be valid)
    - the 'snap points' box is made interactive

I'm not sure if this is the best relationship between the 'detect pour points' and 'outlet features vector' options (perhaps someone would want to use pour points in addition to their outlet geometries?) But it seems the most straightforward for now.
I'll need to update the user's guide as well. Going to wait until we're decided on how this option interacts with the other inputs.
Here's how the UI looks now with and without the box checked:

<img width="786" alt="Screen Shot 2020-09-28 at 11 11 39 AM" src="https://user-images.githubusercontent.com/43770515/94466786-36fda900-017f-11eb-95df-a1d73e511e03.png">
<img width="785" alt="Screen Shot 2020-09-28 at 11 11 51 AM" src="https://user-images.githubusercontent.com/43770515/94466844-4e3c9680-017f-11eb-8888-082c550f80c3.png">




Patch:
diff --git a/HISTORY.rst b/HISTORY.rst
index 1927b5c88..02560ebea 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -67,6 +67,8 @@ Unreleased Changes (3.9)
       'snapped' to be the same as the file name (including suffix).
     * The layer in the 'watersheds.gpkg' output has been renamed from
       'watersheds' to match the name of the vector file (including the suffix).
+    * Added pour point detection option as an alternative to providing an 
+      outlet features vector.
 * Habitat Quality:
     * Refactor of Habitat Quality that implements TaskGraph
     * Threat files are now indicated in the Threat Table csv input under
diff --git a/Makefile b/Makefile
index b2747b82e..13bb00d5e 100644
--- a/Makefile
+++ b/Makefile
@@ -2,11 +2,11 @@
 DATA_DIR := data
 GIT_SAMPLE_DATA_REPO        := https://bitbucket.org/natcap/invest-sample-data.git
 GIT_SAMPLE_DATA_REPO_PATH   := $(DATA_DIR)/invest-sample-data
-GIT_SAMPLE_DATA_REPO_REV    := 611f3d4da7562a903b3e09a826a3dc8884af2f97
+GIT_SAMPLE_DATA_REPO_REV    := 6deacc7df20c3cb7ecd77c05f5cd9277d5e75627
 
 GIT_TEST_DATA_REPO          := https://bitbucket.org/natcap/invest-test-data.git
 GIT_TEST_DATA_REPO_PATH     := $(DATA_DIR)/invest-test-data
-GIT_TEST_DATA_REPO_REV      := f1600516ab5b5704901d6478cca830bdf2fdfb20
+GIT_TEST_DATA_REPO_REV      := b76d594ee9431dcd91f88d66aded990f231e56e9
 
 GIT_UG_REPO                  := https://github.com/natcap/invest.users-guide
 GIT_UG_REPO_PATH             := doc/users-guide
diff --git a/setup.py b/setup.py
index deb4cb935..b84cbaabe 100644
--- a/setup.py
+++ b/setup.py
@@ -49,6 +49,7 @@ setup(
         'natcap',
         'natcap.invest',
         'natcap.invest.coastal_blue_carbon',
+        'natcap.invest.delineateit',
         'natcap.invest.finfish_aquaculture',
         'natcap.invest.fisheries',
         'natcap.invest.hydropower',
@@ -88,6 +89,13 @@ setup(
         'Topic :: Scientific/Engineering :: GIS'
     ],
     ext_modules=[
+        Extension(
+            name="natcap.invest.delineateit.delineateit_core",
+            sources=['src/natcap/invest/delineateit/delineateit_core.pyx'],
+            include_dirs=[numpy.get_include()],
+            extra_compile_args=compiler_and_linker_args,
+            extra_link_args=compiler_and_linker_args,
+            language="c++"),
         Extension(
             name="natcap.invest.recreation.out_of_core_quadtree",
             sources=[
diff --git a/src/natcap/invest/delineateit/__init__.py b/src/natcap/invest/delineateit/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/src/natcap/invest/delineateit.py b/src/natcap/invest/delineateit/delineateit.py
similarity index 77%
rename from src/natcap/invest/delineateit.py
rename to src/natcap/invest/delineateit/delineateit.py
index d6e22d7e9..87df9a55d 100644
--- a/src/natcap/invest/delineateit.py
+++ b/src/natcap/invest/delineateit/delineateit.py
@@ -15,8 +15,9 @@ import pygeoprocessing
 import pygeoprocessing.routing
 import taskgraph
 
-from . import utils
-from . import validation
+from .. import utils
+from .. import validation
+from . import delineateit_core
 
 
 LOGGER = logging.getLogger(__name__)
@@ -44,9 +45,17 @@ ARGS_SPEC = {
                 "each cell."),
             "name": "Digital Elevation Model"
         },
+        "detect_pour_points": {
+            "type": "boolean",
+            "required": False,
+            "about": (
+                "If ``True``, the pour point detection algorithm will run, "
+                "creating a point vector file pour_points.gpkg."),
+            "name": "Detect pour points"
+        },
         "outlet_vector_path": {
             "type": "vector",
-            "required": True,
+            "required": "not detect_pour_points",
             "about": (
                 "This is a layer of geometries representing watershed "
                 "outlets such as municipal water intakes or lakes."),
@@ -110,6 +119,7 @@ _OUTPUT_FILES = {
     'streams': 'streams.tif',
     'snapped_outlets': 'snapped_outlets.gpkg',
     'watersheds': 'watersheds.gpkg',
+    'pour_points': 'pour_points.gpkg'
 }
 
 
@@ -145,7 +155,8 @@ def execute(args):
             'Working with the DEM' section of the InVEST User's Guide for more
             information. (required)
         args['outlet_vector_path'] (string):  This is a vector representing
-            geometries that the watersheds should be built around. (required)
+            geometries that the watersheds should be built around. Required if 
+            ``args['detect_pour_points']`` is False; not used otherwise.
         args['snap_points'] (bool): Whether to snap point geometries to the
             nearest stream pixel.  If ``True``, ``args['flow_threshold']``
             and ``args['snap_distance']`` must also be defined.
@@ -160,6 +171,9 @@ def execute(args):
             If ``False``, this tool will crash if an invalid geometry is
             found.  If ``True``, invalid geometries will be left out of
             the vector to be delineated.  Default: True
+        args['detect_pour_points'] (bool): Whether to run the pour point
+            detection algorithm. If True, detected pour points are used instead 
+            of outlet_vector_path geometries. Default: False
         args['n_workers'] (int): The number of worker processes to use with
             taskgraph. Defaults to -1 (no parallelism).
 
@@ -208,20 +222,36 @@ def execute(args):
         hash_algorithm='md5',
         copy_duplicate_artifact=True)
 
-    check_geometries_task = graph.add_task(
-        check_geometries,
-        args=(args['outlet_vector_path'],
-              file_registry['filled_dem'],
-              file_registry['preprocessed_geometries'],
-              args.get('skip_invalid_geometry', True)),
-        dependent_task_list=[fill_pits_task],
-        target_path_list=[file_registry['preprocessed_geometries']],
-        task_name='check_geometries',
-        hash_algorithm='md5',
-        copy_duplicate_artifact=True)
-    outlet_vector_path = file_registry['preprocessed_geometries']
+    if 'detect_pour_points' in args and args['detect_pour_points']:
+        # Detect pour points automatically and use them instead of
+        # user-provided geometries
+        pour_points_task = graph.add_task(
+            detect_pour_points,
+            args=((file_registry['flow_dir_d8'], 1),
+                  file_registry['pour_points']),
+            dependent_task_list=[flow_dir_task],
+            target_path_list=[file_registry['pour_points']],
+            task_name='detect_pour_points',
+            hash_algorithm='md5',
+            copy_duplicate_artifact=True)
+        outlet_vector_path = file_registry['pour_points']
+        geometry_task = pour_points_task
+    else:
+        check_geometries_task = graph.add_task(
+            check_geometries,
+            args=(args['outlet_vector_path'],
+                  file_registry['filled_dem'],
+                  file_registry['preprocessed_geometries'],
+                  args.get('skip_invalid_geometry', True)),
+            dependent_task_list=[fill_pits_task],
+            target_path_list=[file_registry['preprocessed_geometries']],
+            task_name='check_geometries',
+            hash_algorithm='md5',
+            copy_duplicate_artifact=True)
+        outlet_vector_path = file_registry['preprocessed_geometries']
+        geometry_task = check_geometries_task
 
-    delineation_dependent_tasks = [flow_dir_task, check_geometries_task]
+    delineation_dependent_tasks = [flow_dir_task, geometry_task]
     if 'snap_points' in args and args['snap_points']:
         flow_accumulation_task = graph.add_task(
             pygeoprocessing.routing.flow_accumulation_d8,
@@ -262,7 +292,7 @@ def execute(args):
                   snap_distance,
                   file_registry['snapped_outlets']),
             target_path_list=[file_registry['snapped_outlets']],
-            dependent_task_list=[streams_task, check_geometries_task],
+            dependent_task_list=[streams_task, geometry_task],
             task_name='snapped_outflow_points',
             hash_algorithm='md5',
             copy_duplicate_artifact=True)
@@ -601,6 +631,130 @@ def snap_points_to_nearest_stream(points_vector_path, stream_raster_path_band,
     points_vector = None
 
 
+def detect_pour_points(flow_dir_raster_path_band, target_vector_path):
+    """
+    Create a pour point vector from D8 flow direction raster.
+
+    A pour point is the center point of a pixel which:
+        - flows off of the raster, or
+        - flows into a nodata pixel
+
+    Args:
+        flow_dir_raster_path_band (tuple): tuple of (raster path, band index)
+            indicating the flow direction raster to use. Pixel values are D8
+            values [0 - 7] in this order:
+                321
+                4x0
+                567
+        target_vector_path (string): path to save pour point vector to.
+
+    Returns:
+        None
+    """
+    raster_info = pygeoprocessing.get_raster_info(flow_dir_raster_path_band[0])
+    pour_point_set = _find_raster_pour_points(flow_dir_raster_path_band)
+
+    # use same spatial reference as the input
+    aoi_spatial_reference = osr.SpatialReference()
+    aoi_spatial_reference.ImportFromWkt(raster_info['projection_wkt'])
+
+    gpkg_driver = ogr.GetDriverByName("GPKG")
+    target_vector = gpkg_driver.CreateDataSource(target_vector_path)
+    layer_name = os.path.splitext(
+        os.path.basename(target_vector_path))[0]
+    target_layer = target_vector.CreateLayer(
+        layer_name, aoi_spatial_reference, ogr.wkbPoint)
+    target_defn = target_layer.GetLayerDefn()
+
+    # It's important to have a user-facing unique ID field for post-processing
+    # (e.g. table-joins) that is not the FID. FIDs are not stable across file
+    # conversions that users might do.
+    target_layer.CreateField(
+        ogr.FieldDefn('point_id', ogr.OFTInteger64))
+
+    # Add a feature to the layer for each point
+    target_layer.StartTransaction()
+    for idx, (x, y) in enumerate(pour_point_set):
+        geometry = ogr.Geometry(ogr.wkbPoint)
+        geometry.AddPoint(x, y)
+        feature = ogr.Feature(target_defn)
+        feature.SetGeometry(geometry)
+        feature.SetField('point_id', idx)
+        target_layer.CreateFeature(feature)
+    target_layer.CommitTransaction()
+
+    target_layer = None
+    target_vector = None
+
+
+def _find_raster_pour_points(flow_dir_raster_path_band):
+    """
+    Memory-safe pour point calculation from a flow direction raster.
+    
+    Args:
+        flow_dir_raster_path_band (tuple): tuple of (raster path, band index)
+            indicating the flow direction raster to use.
+
+    Returns:
+        set of (x, y) coordinate tuples of pour points, in the same coordinate
+        system as the input raster.
+    """
+    flow_dir_raster_path, band_index = flow_dir_raster_path_band
+    raster_info = pygeoprocessing.get_raster_info(flow_dir_raster_path)
+    # Open the flow direction raster band
+    raster = gdal.OpenEx(flow_dir_raster_path, gdal.OF_RASTER)
+    if raster is None:
+        raise ValueError(
+            "Raster at %s could not be opened." % flow_dir_raster_path)
+    band = raster.GetRasterBand(band_index)
+    width, height = raster_info['raster_size']
+
+    pour_points = set()
+    # Read in flow direction data and find pour points one block at a time
+    for offsets in pygeoprocessing.iterblocks((flow_dir_raster_path, band_index),
+                                               offset_only=True): 
+        # Expand each block by a one-pixel-wide margin, if possible. 
+        # This way the blocks will overlap so the watershed 
+        # calculation will be continuous.
+        if offsets['xoff'] > 0:
+            offsets['xoff'] -= 1
+            offsets['win_xsize'] += 1
+        if offsets['yoff'] > 0:
+            offsets['yoff'] -= 1
+            offsets['win_ysize'] += 1
+        if offsets['xoff'] + offsets['win_xsize'] < width:
+            offsets['win_xsize'] += 1
+        if offsets['yoff'] + offsets['win_ysize'] < height:
+            offsets['win_ysize'] += 1
+
+        # Keep track of which block edges are raster edges
+        edges = numpy.empty(4, dtype=numpy.intc)
+        # edges order: top, left, bottom, right
+        edges[0] = (offsets['yoff'] == 0)
+        edges[1] = (offsets['xoff'] == 0)
+        edges[2] = (offsets['yoff'] + offsets['win_ysize'] == height)
+        edges[3] = (offsets['xoff'] + offsets['win_xsize'] == width)
+
+        flow_dir_block = band.ReadAsArray(**offsets)
+        pour_points = pour_points.union(
+            delineateit_core.calculate_pour_point_array(
+                # numpy.intc is equivalent to an int in C (normally int32 or 
+                # int64). This way it can be passed directly into a memoryview 
+                # (int[:, :]) in the Cython function.
+                flow_dir_block.astype(numpy.intc), 
+                edges, 
+                nodata=raster_info['nodata'][band_index - 1],
+                offset=(offsets['xoff'], offsets['yoff']),
+                origin=(raster_info['geotransform'][0], 
+                        raster_info['geotransform'][3]),
+                pixel_size=raster_info['pixel_size']))
+            
+    raster = None
+    band = None
+
+    return pour_points
+
+
 @validation.invest_validator
 def validate(args, limit_to=None):
     """Validate args to ensure they conform to `execute`'s contract.
diff --git a/src/natcap/invest/delineateit/delineateit_core.pyx b/src/natcap/invest/delineateit/delineateit_core.pyx
new file mode 100644
index 000000000..979423e62
--- /dev/null
+++ b/src/natcap/invest/delineateit/delineateit_core.pyx
@@ -0,0 +1,109 @@
+# cython: language_level=3
+import numpy
+import pygeoprocessing
+cimport numpy
+cimport cython
+from libcpp.set cimport set as cset
+from libcpp.pair cimport pair as cpair
+
+
+@cython.boundscheck(False)  # Deactivate bounds checking
+@cython.wraparound(False)   # Deactivate negative indexing.
+cpdef cset[cpair[double, double]] calculate_pour_point_array(
+    int[:, :] flow_dir_array, 
+    int[:] edges, 
+    int nodata,
+    cpair[int, int] offset,
+    cpair[double, double] origin,
+    cpair[double, double] pixel_size):
+    """
+    Return a binary array indicating which pixels are pour points.
+
+    Args:
+        flow_dir_array (numpy.ndarray): a 2D array of D8 flow direction 
+            values (0 - 7) and possiby the nodata value.
+        edges (list): a binary list of length 4 indicating whether or not each 
+            edge is an edge of the raster. Order: [top, left, bottom, right]
+        nodata (int): the nodata value of the flow direction array
+        offset (cpair[int, int]): the input flow_dir_array is a block taken 
+            from a larger raster. Offset is the (x, y) coordinate of the
+            upper-left corner of this block relative to the whole raster as a
+            numpy array (the raster starts at (0, 0) and each pixel is 1x1).
+        origin (cpair[double, double]): the (x, y) origin of the raster from 
+            which this block was taken, in its original coordinate system. 
+            This is equivalent to elements (0, 3) of:
+            ``pygeoprocessing.get_raster_info(flow_dir_raster)['geotransform']``.
+        pixel_size (cpair[double, double]): the (x, y) dimensions of a pixel in 
+            the raster from which this block was taken, in its original 
+            coordinate system. This is equivalent to:
+            ``pygeoprocessing.get_raster_info(flow_dir_raster)['pixel_size']``.
+
+    Returns:
+        set of (x, y) coordinates representing pour points in the coordinate
+        system of the original raster. C type cset[cpair[double, double]] is 
+        automatically cast to python type set(tuple(float, float)) when this
+        function is called from a python function.
+    """   
+
+    # flow direction: (delta rows, delta columns)
+    # for a pixel at (i, j) with flow direction x,
+    # the pixel it flows into is located at:
+    # (i + directions[x][0], j + directions[x][1])
+    cdef int *row_offsets = [0, -1, -1, -1,  0,  1, 1, 1]
+    cdef int *col_offsets = [1,  1,  0, -1, -1, -1, 0, 1]
+
+    cdef int height, width
+    cdef int row, col
+    cdef int flow_dir
+    cdef int sink_row, sink_col
+
+    # tuple unpacking doesn't work in cython
+    height, width = flow_dir_array.shape[0], flow_dir_array.shape[1]
+
+    cdef cset[cpair[double, double]] pour_points
+
+    # iterate over each pixel
+    for row in range(height):
+        for col in range(width):
+            
+            # get the flow direction [0-7] for this pixel
+            flow_dir = flow_dir_array[row, col]
+            if flow_dir != nodata:
+                # get the location of the pixel it flows into
+                sink_row = row + row_offsets[flow_dir]
+                sink_col = col + col_offsets[flow_dir]
+
+                # if the row index is -1, the pixel is flowing off of the block
+                if sink_row == -1:
+                    # if this edge of the block is an edge of the raster, 
+                    # this is a pour point. Otherwise, we can't say whether it is
+                    # because the necessary information isn't in this block.
+                    if not edges[0]:  # top edge
+                        continue
+                elif sink_col == -1:
+                    if not edges[1]:  # left edge
+                        continue
+                elif sink_row == height:
+                    if not edges[2]:  # bottom edge
+                        continue
+                elif sink_col == width:
+                    if not edges[3]:  # right edge
+                        continue
+
+                # if we get to here without having continued, the point 
+                # (sink_row, sink_col) is known to be within the bounds of 
+                # the array, so it's safe to index 
+                elif flow_dir_array[sink_row, sink_col] != nodata:
+                    continue
+                
+                # if none of the above conditions passed, this is a pour point
+                pour_points.insert(cpair[double, double](
+                    # +0.5 so that the point is centered in the pixel
+                    (col + offset.first + 0.5) * pixel_size.first + origin.first,
+                    (row + offset.second + 0.5) * pixel_size.second + origin.second))
+
+    # return set of (x, y) coordinates referenced to the same coordinate system
+    # as the original raster
+    return pour_points
+
+
diff --git a/src/natcap/invest/hydropower/hydropower_water_yield.py b/src/natcap/invest/hydropower/hydropower_water_yield.py
index 730a903ea..79af353fc 100644
--- a/src/natcap/invest/hydropower/hydropower_water_yield.py
+++ b/src/natcap/invest/hydropower/hydropower_water_yield.py
@@ -959,8 +959,8 @@ def compute_watershed_valuation(watershed_results_vector_path, val_dict):
 def compute_rsupply_volume(watershed_results_vector_path):
     """Calculate the total realized water supply volume.
 
-    And the mean realized water supply volume per hectare for the given sheds.
-    Output units in cubic meters and cubic meters per hectare respectively.
+    And the mean realized water supply volume per pixel for the given sheds.
+    Output units in cubic meters and cubic meters per pixel respectively.
 
     Args:
         watershed_results_vector_path (string): a path to a vector that
diff --git a/src/natcap/invest/ui/delineateit.py b/src/natcap/invest/ui/delineateit.py
index fb93195a7..6bc1e09c3 100644
--- a/src/natcap/invest/ui/delineateit.py
+++ b/src/natcap/invest/ui/delineateit.py
@@ -1,7 +1,7 @@
 # encoding=UTF-8
 
 from natcap.invest.ui import model, inputs
-from natcap.invest import delineateit
+from natcap.invest.delineateit import delineateit
 
 
 class Delineateit(model.InVESTModel):
@@ -21,6 +21,18 @@ class Delineateit(model.InVESTModel):
                 "for each cell."),
             validator=self.validator)
         self.add_input(self.dem_path)
+
+        self.detect_pour_points = inputs.Checkbox(
+            args_key='detect_pour_points',
+            helptext=(
+                'If this box is checked, run pour point detection algorithm. '
+                'Use detected pour points instead of outlet geometries.'),
+            label='Detect pour points'
+        )
+        self.add_input(self.detect_pour_points)
+        self.detect_pour_points.sufficiency_changed.connect(
+            self.on_detect_pour_points_changed)
+
         self.outlet_vector_path = inputs.File(
             args_key='outlet_vector_path',
             helptext=(
@@ -73,6 +85,12 @@ class Delineateit(model.InVESTModel):
             validator=self.validator)
         self.snap_points_container.add_input(self.snap_distance)
 
+    def on_detect_pour_points_changed(self, value):
+        """Change interactivity of other inputs given boolean signal value"""
+        self.outlet_vector_path.set_interactive(not value)
+        self.skip_invalid_geometry.set_interactive(not value)
+        self.snap_points_container.set_interactive(value)
+
     def _enable_point_snapping_container(self, input_valid):
         outlet_vector_path = self.outlet_vector_path.value()
         if delineateit._vector_may_contain_points(outlet_vector_path):
@@ -86,8 +104,7 @@ class Delineateit(model.InVESTModel):
             self.workspace.args_key: self.workspace.value(),
             self.suffix.args_key: self.suffix.value(),
             self.dem_path.args_key: self.dem_path.value(),
-            self.outlet_vector_path.args_key: (
-                self.outlet_vector_path.value()),
+            self.detect_pour_points.args_key: self.detect_pour_points.value(),
             self.snap_points_container.args_key: (
                 self.snap_points_container.value()),
             self.flow_threshold.args_key: self.flow_threshold.value(),
@@ -95,5 +112,11 @@ class Delineateit(model.InVESTModel):
             self.skip_invalid_geometry.args_key: (
                 self.skip_invalid_geometry.value()),
         }
+        # If the outlet_vector_path input is grayed out but still has a value,
+        # validation will still check that it overlaps the DEM.
+        # To avoid this, only include it if it's going to be used.
+        if self.detect_pour_points.value() == False:
+            args[self.outlet_vector_path.args_key] = self.outlet_vector_path.value()
+
 
         return args
diff --git a/src/natcap/invest/ui/usage.py b/src/natcap/invest/ui/usage.py
index 80c42b958..ccdc0b30d 100644
--- a/src/natcap/invest/ui/usage.py
+++ b/src/natcap/invest/ui/usage.py
@@ -10,6 +10,7 @@ import sys
 import threading
 import traceback
 import uuid
+import importlib
 
 from urllib.request import urlopen, Request
 from urllib.parse import urlencode
@@ -65,11 +66,12 @@ def log_run(module, args):
         log_exit_thread.start()
 
 
-def _calculate_args_bounding_box(args_dict):
+def _calculate_args_bounding_box(args, args_spec):
     """Calculate the bounding boxes of any GIS types found in `args_dict`.
 
     Args:
-        args_dict (dict): a string key and any value pair dictionary.
+        args (dict): a string key and any value pair dictionary.
+        args_spec (dict): the model ARGS_SPEC describing args
 
     Returns:
         bb_intersection, bb_union tuple that's either the lat/lng bounding
@@ -77,6 +79,7 @@ def _calculate_args_bounding_box(args_dict):
             in args_dict.  If no GIS types are present, this is a (None, None)
             tuple.
     """
+
     def _merge_bounding_boxes(bb1, bb2, mode):
         """Merge two bounding boxes through union or intersection.
 
@@ -106,92 +109,53 @@ def _calculate_args_bounding_box(args_dict):
         bb_out = [op(x, y) for op, x, y in zip(comparison_ops, bb1, bb2)]
         return bb_out
 
-    def _merge_local_bounding_boxes(arg, bb_intersection=None, bb_union=None):
-        """Traverse nested dictionary to merge bounding boxes of GIS types.
-
-        Args:
-            arg (dict): contains string keys and pairs that might be files to
-                gis types.  They can be any other type, including dictionaries.
-            bb_intersection (list or None): if list, has the form
-                [xmin, ymin, xmax, ymax], where coordinates are in lng, lat
-            bb_union (list or None): if list, has the form
-                [xmin, ymin, xmax, ymax], where coordinates are in lng, lat
-
-        Returns:
-            (intersection, union) bounding box tuples of all filepaths to GIS
-            data types found in the dictionary and bb_intersection and bb_union
-            inputs.  None, None if no arguments were GIS data types and input
-            bounding boxes are None.
-        """
-        def _is_spatial(arg):
-            if isinstance(arg, str) and os.path.exists(arg):
-                with utils.capture_gdal_logging():
-                    dataset = gdal.OpenEx(arg)
-                    if dataset is not None:
-                        # OGR opens CSV files.  For now, we should not
-                        # consider these to be vectors.
-                        driver_name = dataset.GetDriver().ShortName
-                        if driver_name == 'CSV':
-                            return False
-                        return True
-            return False
-
-        if isinstance(arg, dict):
-            # if dict, grab the bb's for all the members in it
-            for value in arg.values():
-                bb_intersection, bb_union = _merge_local_bounding_boxes(
-                    value, bb_intersection, bb_union)
-        elif isinstance(arg, list):
-            # if list, grab the bb's for all the members in it
-            for value in arg:
-                bb_intersection, bb_union = _merge_local_bounding_boxes(
-                    value, bb_intersection, bb_union)
+    bb_intersection = None
+    bb_union = None
+    for key, value in args.items():
+        # Using gdal.OpenEx to check if an input is spatial caused the
+        # model to hang sometimes (possible race condition), so only
+        # get the bounding box of inputs that are known to be spatial.
+        spatial_info = None
+        if args_spec['args'][key]['type'] == 'raster':
+            spatial_info = pygeoprocessing.get_raster_info(value)
+        elif args_spec['args'][key]['type'] == 'vector':
+            spatial_info = pygeoprocessing.get_vector_info(value)
+                    
+        if spatial_info:
+            local_bb = spatial_info['bounding_box']
+            projection_wkt = spatial_info['projection_wkt']
+            spatial_ref = osr.SpatialReference()
+            spatial_ref.ImportFromWkt(projection_wkt)
+
+            try:
+                # means there's a GIS type with a well defined bounding box
+                # create transform, and reproject local bounding box to
+                # lat/lng
+                lat_lng_ref = osr.SpatialReference()
+                lat_lng_ref.ImportFromEPSG(4326)  # EPSG 4326 is lat/lng
+                to_lat_trans = utils.create_coordinate_transformer(
+                    spatial_ref, lat_lng_ref)
+                for point_index in [0, 2]:
+                    local_bb[point_index], local_bb[point_index + 1], _ = (
+                        to_lat_trans.TransformPoint(
+                            local_bb[point_index],
+                            local_bb[point_index+1]))
+
+                bb_intersection = _merge_bounding_boxes(
+                    local_bb, bb_intersection, 'intersection')
+                bb_union = _merge_bounding_boxes(
+                    local_bb, bb_union, 'union')
+            except Exception as transform_error:
+                # All kinds of exceptions from bad transforms or CSV files
+                # or dbf files could get us to this point, just don't
+                # bother with the local_bb at all
+                LOGGER.exception('Error when transforming coordinates: %s',
+                                 transform_error)
         else:
-            # singular value, test if GIS type, if not, don't update bb's
-            # this is an undefined bounding box that gets returned when ogr
-            # opens a table only
-            if _is_spatial(arg):
-                with utils.capture_gdal_logging():
-                    if gdal.OpenEx(arg, gdal.OF_RASTER) is not None:
-                        spatial_info = pygeoprocessing.get_raster_info(arg)
-                    else:
-                        # If it isn't a raster, it should be a vector!
-                        spatial_info = pygeoprocessing.get_vector_info(arg)
-
-                local_bb = [0., 0., 0., 0.]
-                local_bb = spatial_info['bounding_box']
-                projection_wkt = spatial_info['projection_wkt']
-                spatial_ref = osr.SpatialReference()
-                spatial_ref.ImportFromWkt(projection_wkt)
-
-                try:
-                    # means there's a GIS type with a well defined bounding box
-                    # create transform, and reproject local bounding box to
-                    # lat/lng
-                    lat_lng_ref = osr.SpatialReference()
-                    lat_lng_ref.ImportFromEPSG(4326)  # EPSG 4326 is lat/lng
-                    to_lat_trans = utils.create_coordinate_transformer(
-                        spatial_ref, lat_lng_ref)
-                    for point_index in [0, 2]:
-                        local_bb[point_index], local_bb[point_index + 1], _ = (
-                            to_lat_trans.TransformPoint(
-                                local_bb[point_index],
-                                local_bb[point_index+1]))
-
-                    bb_intersection = _merge_bounding_boxes(
-                        local_bb, bb_intersection, 'intersection')
-                    bb_union = _merge_bounding_boxes(
-                        local_bb, bb_union, 'union')
-                except Exception as transform_error:
-                    # All kinds of exceptions from bad transforms or CSV files
-                    # or dbf files could get us to this point, just don't
-                    # bother with the local_bb at all
-                    LOGGER.exception('Error when transforming coordinates: %s',
-                                     transform_error)
-
-        return bb_intersection, bb_union
-
-    return _merge_local_bounding_boxes(args_dict)
+            LOGGER.debug(f'Arg {key} of type {args_spec["args"][key]["type"]} '
+                          'excluded from bounding box calculation')
+
+    return bb_intersection, bb_union
 
 
 def _log_exit_status(session_id, status):
@@ -252,9 +216,11 @@ def _log_model(model_name, model_args, session_id=None):
         md5.update(json.dumps(data).encode('utf-8'))
         return md5.hexdigest()
 
+    args_spec = importlib.import_module(model_name).ARGS_SPEC
+
     try:
         bounding_box_intersection, bounding_box_union = (
-            _calculate_args_bounding_box(model_args))
+            _calculate_args_bounding_box(model_args, args_spec))
 
         payload = {
             'model_name': model_name,


Test patch:
diff --git a/tests/test_delineateit.py b/tests/test_delineateit.py
index e97f05077..522eaa41d 100644
--- a/tests/test_delineateit.py
+++ b/tests/test_delineateit.py
@@ -1,5 +1,6 @@
 """Module for Testing DelineateIt."""
 import unittest
+from unittest import mock
 import tempfile
 import shutil
 import os
@@ -34,7 +35,7 @@ class DelineateItTests(unittest.TestCase):
 
     def test_delineateit_willamette(self):
         """DelineateIt: regression testing full run."""
-        from natcap.invest import delineateit
+        from natcap.invest.delineateit import delineateit
 
         args = {
             'dem_path': os.path.join(REGRESSION_DATA, 'input', 'dem.tif'),
@@ -70,17 +71,20 @@ class DelineateItTests(unittest.TestCase):
 
     def test_delineateit_validate(self):
         """DelineateIt: test validation function."""
-        from natcap.invest import delineateit
+        from natcap.invest.delineateit import delineateit
         missing_keys = {}
         validation_warnings = delineateit.validate(missing_keys)
-        self.assertEqual(len(validation_warnings), 1)
-        self.assertEqual(len(validation_warnings[0][0]), 3)
+        self.assertEqual(len(validation_warnings), 2)
+        self.assertTrue('dem_path' in validation_warnings[0][0])
+        self.assertTrue('workspace_dir' in validation_warnings[0][0])
+        self.assertEqual(validation_warnings[1][0], ['outlet_vector_path'])
 
         missing_values_args = {
             'workspace_dir': '',
             'dem_path': None,
             'outlet_vector_path': '',
             'snap_points': False,
+            'detect_pour_points': True
         }
         validation_warnings = delineateit.validate(missing_values_args)
         self.assertEqual(len(validation_warnings), 1)
@@ -125,7 +129,7 @@ class DelineateItTests(unittest.TestCase):
 
     def test_point_snapping(self):
         """DelineateIt: test point snapping."""
-        from natcap.invest import delineateit
+        from natcap.invest.delineateit import delineateit
 
         srs = osr.SpatialReference()
         srs.ImportFromEPSG(32731)  # WGS84/UTM zone 31s
@@ -201,7 +205,7 @@ class DelineateItTests(unittest.TestCase):
 
     def test_vector_may_contain_points(self):
         """DelineateIt: Check whether a layer contains points."""
-        from natcap.invest import delineateit
+        from natcap.invest.delineateit import delineateit
 
         srs = osr.SpatialReference()
         srs.ImportFromEPSG(32731)  # WGS84/UTM zone 31s
@@ -239,7 +243,7 @@ class DelineateItTests(unittest.TestCase):
 
     def test_check_geometries(self):
         """DelineateIt: Check that we can reasonably repair geometries."""
-        from natcap.invest import delineateit
+        from natcap.invest.delineateit import delineateit
         srs = osr.SpatialReference()
         srs.ImportFromEPSG(32731)  # WGS84/UTM zone 31s
         projection_wkt = srs.ExportToWkt()
@@ -345,3 +349,108 @@ class DelineateItTests(unittest.TestCase):
 
         target_layer = None
         target_vector = None
+
+    def test_detect_pour_points(self):
+        from natcap.invest.delineateit import delineateit
+ 
+        # create a flow direction raster from the sample DEM
+        flow_dir_path = os.path.join(REGRESSION_DATA, 'input/flow_dir_gura.tif')
+        output_path = os.path.join(self.workspace_dir, 'point_vector.gpkg')
+
+        delineateit.detect_pour_points((flow_dir_path, 1), output_path)
+
+        vector = gdal.OpenEx(output_path, gdal.OF_VECTOR)
+        layer = vector.GetLayer()
+
+        points = []
+        for feature in layer:
+            geom = feature.GetGeometryRef()
+            x, y, _ = geom.GetPoint()
+            points.append((x, y))
+        points.sort()
+
+        expected_points = [
+            (277713.1562500205, 9941874.499999935),
+            (277713.1562500205, 9941859.499999935)]
+
+        self.assertTrue(numpy.isclose(points[0][0], expected_points[0][0]))
+        self.assertTrue(numpy.isclose(points[0][1], expected_points[0][1]))
+        self.assertTrue(numpy.isclose(points[1][0], expected_points[1][0]))
+        self.assertTrue(numpy.isclose(points[1][1], expected_points[1][1]))
+
+    def test_calculate_pour_point_array(self):
+        from natcap.invest.delineateit import delineateit, delineateit_core
+
+        a = 100  # nodata value
+
+        # at this point the flow direction array will already have been
+        # buffered with a border of nodata
+        flow_dir_array = numpy.array([
+            [7, 7, 7, 5],
+            [6, 6, 6, 4],
+            [0, 1, a, 0],
+            [a, 2, 3, 4]
+        ], dtype=numpy.intc)
+
+        edges = numpy.array([1, 1, 0, 0], dtype=numpy.intc)  # top, left, bottom, right
+
+        expected_pour_points = {(6.25, 5.25)}
+
+        output = delineateit_core.calculate_pour_point_array(
+            flow_dir_array, 
+            edges, 
+            nodata=a,
+            offset=(0, 0),
+            origin=(5, 3),
+            pixel_size=(0.5, 1.5))
+
+        self.assertTrue(numpy.array_equal(
+            output,
+            expected_pour_points))
+
+    def test_find_pour_points_by_block(self):
+        from natcap.invest.delineateit import delineateit
+
+        a = 100  # nodata value
+        flow_dir_array = numpy.array([
+            [0, 0, 0, 0, 7, 7, 7, 1, 6, 6],
+            [2, 3, 4, 5, 6, 7, 0, 1, 1, 2],
+            [2, 2, 2, 2, 0, a, a, 3, 3, a],
+            [2, 1, 1, 1, 2, 6, 4, 1, a, a],
+            [1, 1, 0, 0, 0, 0, a, a, a, a]
+        ], dtype=numpy.int8)
+
+        srs = osr.SpatialReference()
+        srs.ImportFromEPSG(3157)
+        projection_wkt = srs.ExportToWkt()
+
+        raster_path = os.path.join(self.workspace_dir, 'small_raster.tif')
+        pygeoprocessing.numpy_array_to_raster(
+            flow_dir_array,
+            a,
+            (1, 1),
+            (0, 0),
+            projection_wkt,
+            raster_path
+        )
+
+        expected_pour_points = {(7.5, 0.5), (5.5, 1.5), (4.5, 2.5), (5.5, 4.5)}
+
+        # Mock iterblocks so that we can test with an array smaller than 128x128
+        # to make sure that the algorithm gets pour points on block edges e.g.
+        # flow_dir_array[2, 4]
+        def mock_iterblocks(*args, **kwargs):
+            xoffs = [0, 4, 8]
+            win_xsizes = [4, 4, 2]
+            for xoff, win_xsize in zip(xoffs, win_xsizes):
+                yield {
+                    'xoff': xoff, 
+                    'yoff': 0,
+                    'win_xsize': win_xsize,
+                    'win_ysize': 5}
+
+        with mock.patch(
+            'natcap.invest.delineateit.delineateit.pygeoprocessing.iterblocks', 
+            mock_iterblocks):
+            pour_points = delineateit._find_raster_pour_points((raster_path, 1))
+            self.assertEqual(pour_points, expected_pour_points)
diff --git a/ui_tests/test_model_logging.py b/ui_tests/test_model_logging.py
index 4c29103c9..74297bb02 100644
--- a/ui_tests/test_model_logging.py
+++ b/ui_tests/test_model_logging.py
@@ -61,9 +61,18 @@ class ModelLoggingTests(unittest.TestCase):
             'not_a_gis_input': 'foobar'
         }
 
+        args_spec = {
+            'args': {
+                'raster': {'type': 'raster'},
+                'vector': {'type': 'vector'},
+                'not_a_gis_input': {'type': 'freestyle_string'}
+            }
+        }
+
         output_logfile = os.path.join(self.workspace_dir, 'logfile.txt')
         with utils.log_to_file(output_logfile):
-            bb_inter, bb_union = usage._calculate_args_bounding_box(model_args)
+            bb_inter, bb_union = usage._calculate_args_bounding_box(
+                model_args, args_spec)
 
         numpy.testing.assert_allclose(
             bb_inter, [-87.234108, -85.526151, -87.233424, -85.526205])