Repository: natcap/invest (python)

Issue:
DelineateIt should use flow accumulation to break ties in snapping points to streams
During point snapping, if we have two pixels with the same snapping distance, we should use the Flow Accumulation raster to break the tie.  The stream pixel with the highest flow accumulation should be picked over those that are the same distance but with lower flow accumulation values.

Pull request:
DelineateIt: break ties in point snapping with flow accumulation
# Description
Fixes #60 
Also made a couple of minor code style improvements.

# Checklist
- [x] Updated HISTORY.rst (if these changes are user-facing)
- [x] Updated the user's guide (if needed)


Patch:
diff --git a/HISTORY.rst b/HISTORY.rst
index 7d73fa166..23855c8e8 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -51,6 +51,11 @@ Unreleased Changes
     * Updating the ``taskgraph`` requirement to ``0.11.0`` to resolve an issue
       where modifying a file within a roughly 2-second window would fool
       ``taskgraph`` into believing that the file had not been modified.
+* DelineateIt
+    * When snapping points to streams, if a point is equally near to more than
+      one stream pixel, it will now snap to the stream pixel with a higher
+      flow accumulation value. Before, it would snap to the stream pixel
+      encountered first in the raster (though this was not guaranteed).
 * Habitat Quality
     * Changed how Habitat Rarity outputs are calculated to be less confusing.
       Values now represent a 0 to 1 index where before there could be
diff --git a/Makefile b/Makefile
index d06f774e7..0eae4bed9 100644
--- a/Makefile
+++ b/Makefile
@@ -10,7 +10,7 @@ GIT_TEST_DATA_REPO_REV      := 586a64b2b33a124d9197bd99c06c88b926807f70
 
 GIT_UG_REPO                 := https://github.com/natcap/invest.users-guide
 GIT_UG_REPO_PATH            := doc/users-guide
-GIT_UG_REPO_REV             := 764aa04944319f3af79004860aa35c3be189b264
+GIT_UG_REPO_REV             := be19f67350e7663397a7a05408c429e763882051
 
 ENV = "./env"
 ifeq ($(OS),Windows_NT)
diff --git a/requirements-dev.txt b/requirements-dev.txt
index 6dec133c6..ab4cfd6fe 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -18,6 +18,7 @@ pypiwin32; sys_platform == 'win32'  # pip-only
 setuptools>=8.0
 Sphinx>=1.3.1,!=1.7.1
 sphinx-rtd-theme
+sphinx-reredirects
 PyInstaller>=4.1
 setuptools_scm
 requests
diff --git a/requirements-docs.yml b/requirements-docs.yml
index 102bb290c..e72a4e5a9 100644
--- a/requirements-docs.yml
+++ b/requirements-docs.yml
@@ -34,4 +34,5 @@ dependencies:
   - qtawesome  # pip-only
   - setuptools_scm
   - sphinx-rtd-theme
+  - sphinx-reredirects
 
diff --git a/src/natcap/invest/delineateit/delineateit.py b/src/natcap/invest/delineateit/delineateit.py
index a17f7c36e..281d74da4 100644
--- a/src/natcap/invest/delineateit/delineateit.py
+++ b/src/natcap/invest/delineateit/delineateit.py
@@ -67,7 +67,10 @@ ARGS_SPEC = {
             "about": (
                 "Whether to snap point geometries to the nearest stream "
                 "pixel.  If ``True``, ``args['flow_threshold']`` and "
-                "``args['snap_distance']`` must also be defined."),
+                "``args['snap_distance']`` must also be defined. If a point "
+                "is equally near to more than one stream pixel, it will be "
+                "snapped to the stream pixel with the highest flow "
+                "accumulation value."),
             "name": "Snap points to the nearest stream"
         },
         "flow_threshold": {
@@ -146,14 +149,14 @@ def execute(args):
             selected folder does not exist, it will be created. If datasets
             already exist in the selected folder, they will be overwritten.
             (required)
-        args['results_suffix'] (string):  This text will be appended to the end of
-            output files to help separate multiple runs. (optional)
-        args['dem_path'] (string):  A GDAL-supported raster file with an elevation
-            for each cell. Make sure the DEM is corrected by filling in sinks,
-            and if necessary burning hydrographic features into the elevation
-            model (recommended when unusual streams are observed.) See the
-            'Working with the DEM' section of the InVEST User's Guide for more
-            information. (required)
+        args['results_suffix'] (string):  This text will be appended to the end
+            of output files to help separate multiple runs. (optional)
+        args['dem_path'] (string):  A GDAL-supported raster file with an
+            elevation for each cell. Make sure the DEM is corrected by filling
+            in sinks, and if necessary burning hydrographic features into the
+            elevation model (recommended when unusual streams are observed.)
+            See the '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 if
             ``args['detect_pour_points']`` is False; not used otherwise.
@@ -278,7 +281,8 @@ def execute(args):
         snapped_outflow_points_task = graph.add_task(
             snap_points_to_nearest_stream,
             args=(outlet_vector_path,
-                  (file_registry['streams'], 1),
+                  file_registry['streams'],
+                  file_registry['flow_accumulation'],
                   snap_distance,
                   file_registry['snapped_outlets']),
             target_path_list=[file_registry['snapped_outlets']],
@@ -287,7 +291,7 @@ def execute(args):
         delineation_dependent_tasks.append(snapped_outflow_points_task)
         outlet_vector_path = file_registry['snapped_outlets']
 
-    watershed_delineation_task = graph.add_task(
+    _ = graph.add_task(
         pygeoprocessing.routing.delineate_watersheds_d8,
         args=((file_registry['flow_dir_d8'], 1),
               outlet_vector_path,
@@ -460,8 +464,9 @@ def check_geometries(outlet_vector_path, dem_path, target_vector_path,
     target_vector = None
 
 
-def snap_points_to_nearest_stream(points_vector_path, stream_raster_path_band,
-                                  snap_distance, snapped_points_vector_path):
+def snap_points_to_nearest_stream(points_vector_path, stream_raster_path,
+                                  flow_accum_raster_path, snap_distance,
+                                  snapped_points_vector_path):
     """Adjust the location of points to the nearest stream pixel.
 
     The new point layer will have all fields and field values copied over from
@@ -472,9 +477,12 @@ def snap_points_to_nearest_stream(points_vector_path, stream_raster_path_band,
         points_vector_path (string): A path to a vector on disk containing
             point geometries.  Must be in the same projection as the stream
             raster.
-        stream_raster_path_band (tuple): A tuple of (path, band index), where
+        stream_raster_path (string): A path to a stream raster, where
             pixel values are ``1`` (indicating a stream pixel) or ``0``
             (indicating a non-stream pixel).
+        flow_accum_raster_path (string): A path to a flow accumulation raster
+            that is aligned with the stream raster. Used to break ties between
+            equally-near stream pixels.
         snap_distance (number): The maximum distance (in pixels) to search
             for stream pixels for each point.  This must be a positive, nonzero
             value.
@@ -494,12 +502,14 @@ def snap_points_to_nearest_stream(points_vector_path, stream_raster_path_band,
     points_vector = gdal.OpenEx(points_vector_path, gdal.OF_VECTOR)
     points_layer = points_vector.GetLayer()
 
-    stream_raster_info = pygeoprocessing.get_raster_info(
-        stream_raster_path_band[0])
+    stream_raster_info = pygeoprocessing.get_raster_info(stream_raster_path)
     geotransform = stream_raster_info['geotransform']
     n_cols, n_rows = stream_raster_info['raster_size']
-    stream_raster = gdal.OpenEx(stream_raster_path_band[0], gdal.OF_RASTER)
-    stream_band = stream_raster.GetRasterBand(stream_raster_path_band[1])
+    stream_raster = gdal.OpenEx(stream_raster_path, gdal.OF_RASTER)
+    stream_band = stream_raster.GetRasterBand(1)
+
+    flow_accum_raster = gdal.OpenEx(flow_accum_raster_path, gdal.OF_RASTER)
+    flow_accum_band = flow_accum_raster.GetRasterBand(1)
 
     driver = gdal.GetDriverByName('GPKG')
     snapped_vector = driver.Create(snapped_points_vector_path, 0, 0, 0,
@@ -560,8 +570,6 @@ def snap_points_to_nearest_stream(points_vector_path, stream_raster_path_band,
                 f'stream raster.  FID:{point_feature.GetFID()} at {point}')
             continue
 
-        x_center = x_index
-        y_center = y_index
         x_left = max(x_index - snap_distance, 0)
         y_top = max(y_index - snap_distance, 0)
         x_right = min(x_index + snap_distance, n_cols)
@@ -573,20 +581,37 @@ def snap_points_to_nearest_stream(points_vector_path, stream_raster_path_band,
             int(y_bottom - y_top))
         row_indexes, col_indexes = numpy.nonzero(
             stream_window == 1)
-        if row_indexes.size > 0:
-            # Calculate euclidean distance for sorting
+
+        # Find the closest stream pixel that meets the distance
+        # requirement. If there is a tie, snap to the stream pixel with
+        # a higher flow accumulation value.
+        if row_indexes.size > 0:  # there are streams within the snap distance
+            # Calculate euclidean distance from the point to each stream pixel
             distance_array = numpy.hypot(
-                y_center - y_top - row_indexes,
-                x_center - x_left - col_indexes,
+                # distance along y axis from the point to each stream pixel
+                y_index - y_top - row_indexes,
+                # distance along x axis from the point to each stream pixel
+                x_index - x_left - col_indexes,
                 dtype=numpy.float32)
 
-            # Find the closest stream pixel that meets the distance
-            # requirement.
-            min_index = numpy.argmin(distance_array)
-            min_row = row_indexes[min_index]
-            min_col = col_indexes[min_index]
-            offset_row = min_row - (y_center - y_top)
-            offset_col = min_col - (x_center - x_left)
+            is_nearest = distance_array == distance_array.min()
+            # if > 1 stream pixel is nearest, break tie with flow accumulation
+            if is_nearest.sum() > 1:
+                flow_accum_array = flow_accum_band.ReadAsArray(
+                    int(x_left), int(y_top), int(x_right - x_left),
+                    int(y_bottom - y_top))
+                # weight by flow accum
+                is_nearest = is_nearest * flow_accum_array[row_indexes, col_indexes]
+
+            # 1d index of max value in flattened array
+            nearest_stream_index_1d = numpy.argmax(is_nearest)
+
+            # convert 1d index back to coordinates relative to window
+            nearest_stream_row = row_indexes[nearest_stream_index_1d]
+            nearest_stream_col = col_indexes[nearest_stream_index_1d]
+
+            offset_row = nearest_stream_row - (y_index - y_top)
+            offset_col = nearest_stream_col - (x_index - x_left)
 
             y_index += offset_row
             x_index += offset_col


Test patch:
diff --git a/tests/test_delineateit.py b/tests/test_delineateit.py
index 71d9e5ced..a561ccdb8 100644
--- a/tests/test_delineateit.py
+++ b/tests/test_delineateit.py
@@ -183,9 +183,20 @@ class DelineateItTests(unittest.TestCase):
              [0, 1, 0, 0, 0, 0],
              [0, 1, 0, 0, 0, 0]], dtype=numpy.int8)
         stream_raster_path = os.path.join(self.workspace_dir, 'streams.tif')
+        flow_accum_array = numpy.array(
+            [[1, 5, 1, 1, 1, 1],
+             [1, 5, 1, 1, 1, 1],
+             [1, 5, 1, 1, 1, 1],
+             [1, 5, 1, 1, 1, 1],
+             [1, 5, 5, 5, 5, 5],
+             [1, 5, 1, 1, 1, 1],
+             [1, 5, 1, 1, 1, 1]], dtype=numpy.int8)
+        flow_accum_path = os.path.join(self.workspace_dir, 'flow_accum.tif')
         # byte datatype
         pygeoprocessing.numpy_array_to_raster(
             stream_matrix, 255, (2, -2), (2, -2), wkt, stream_raster_path)
+        pygeoprocessing.numpy_array_to_raster(
+            flow_accum_array, 255, (2, -2), (2, -2), wkt, flow_accum_path)
 
         source_points_path = os.path.join(self.workspace_dir,
                                           'source_features.geojson')
@@ -216,13 +227,13 @@ class DelineateItTests(unittest.TestCase):
         snap_distance = -1
         with self.assertRaises(ValueError) as cm:
             delineateit.snap_points_to_nearest_stream(
-                source_points_path, (stream_raster_path, 1),
+                source_points_path, stream_raster_path, flow_accum_path,
                 snap_distance, snapped_points_path)
         self.assertTrue('must be >= 0' in str(cm.exception))
 
         snap_distance = 10  # large enough to get multiple streams per point.
         delineateit.snap_points_to_nearest_stream(
-            source_points_path, (stream_raster_path, 1),
+            source_points_path, stream_raster_path, flow_accum_path,
             snap_distance, snapped_points_path)
 
         snapped_points_vector = gdal.OpenEx(snapped_points_path,
@@ -265,9 +276,20 @@ class DelineateItTests(unittest.TestCase):
              [0, 1, 0, 0, 0, 0],
              [0, 1, 0, 0, 0, 0]], dtype=numpy.int8)
         stream_raster_path = os.path.join(self.workspace_dir, 'streams.tif')
+        flow_accum_array = numpy.array(
+            [[1, 5, 1, 1, 1, 1],
+             [1, 5, 1, 1, 1, 1],
+             [1, 5, 1, 1, 1, 1],
+             [1, 5, 1, 1, 1, 1],
+             [1, 5, 5, 5, 5, 5],
+             [1, 5, 1, 1, 1, 1],
+             [1, 5, 1, 1, 1, 1]], dtype=numpy.int8)
+        flow_accum_path = os.path.join(self.workspace_dir, 'flow_accum.tif')
         # byte datatype
         pygeoprocessing.numpy_array_to_raster(
             stream_matrix, 255, (2, -2), (2, -2), wkt, stream_raster_path)
+        pygeoprocessing.numpy_array_to_raster(
+            flow_accum_array, 255, (2, -2), (2, -2), wkt, flow_accum_path)
 
         source_points_path = os.path.join(self.workspace_dir,
                                           'source_features.gpkg')
@@ -304,7 +326,7 @@ class DelineateItTests(unittest.TestCase):
                                            'snapped_points.gpkg')
         snap_distance = 10  # large enough to get multiple streams per point.
         delineateit.snap_points_to_nearest_stream(
-            source_points_path, (stream_raster_path, 1),
+            source_points_path, stream_raster_path, flow_accum_path,
             snap_distance, snapped_points_path)
 
         try:
@@ -324,6 +346,66 @@ class DelineateItTests(unittest.TestCase):
             snapped_points_layer = None
             snapped_points_vector = None
 
+    def test_point_snapping_break_ties(self):
+        """DelineateIt: distance ties are broken using flow accumulation."""
+        from natcap.invest.delineateit import delineateit
+
+        srs = osr.SpatialReference()
+        srs.ImportFromEPSG(32731)  # WGS84/UTM zone 31s
+        wkt = srs.ExportToWkt()
+
+        # need stream layer, points
+        stream_matrix = numpy.array(
+            [[0, 1, 0, 0, 0, 0],
+             [0, 1, 0, 0, 0, 0],
+             [0, 1, 0, 0, 0, 0],
+             [0, 1, 0, 0, 0, 0],
+             [0, 1, 1, 1, 1, 1],
+             [0, 1, 0, 0, 0, 0],
+             [0, 1, 0, 0, 0, 0]], dtype=numpy.int8)
+        stream_raster_path = os.path.join(self.workspace_dir, 'streams.tif')
+        flow_accum_array = numpy.array(
+            [[1, 5, 1, 1, 1, 1],
+             [1, 5, 1, 1, 1, 1],
+             [1, 5, 1, 1, 1, 1],
+             [1, 5, 1, 1, 1, 1],
+             [1, 5, 9, 9, 9, 9],
+             [1, 4, 1, 1, 1, 1],
+             [1, 4, 1, 1, 1, 1]], dtype=numpy.int8)
+        flow_accum_path = os.path.join(self.workspace_dir, 'flow_accum.tif')
+        pygeoprocessing.numpy_array_to_raster(
+            stream_matrix, 255, (2, -2), (2, -2), wkt, stream_raster_path)
+        pygeoprocessing.numpy_array_to_raster(
+            flow_accum_array, -1, (2, -2), (2, -2), wkt, flow_accum_path)
+
+        source_points_path = os.path.join(self.workspace_dir,
+                                          'source_features.geojson')
+        source_features = [Point(9, -7)]  # equidistant from two streams
+        pygeoprocessing.shapely_geometry_to_vector(
+            source_features, source_points_path, wkt, 'GeoJSON',
+            ogr_geom_type=ogr.wkbUnknown)
+
+        snapped_points_path = os.path.join(self.workspace_dir,
+                                           'snapped_points.gpkg')
+
+        snap_distance = 10  # large enough to get multiple streams per point.
+        delineateit.snap_points_to_nearest_stream(
+            source_points_path, stream_raster_path, flow_accum_path,
+            snap_distance, snapped_points_path)
+
+        snapped_points_vector = gdal.OpenEx(snapped_points_path,
+                                            gdal.OF_VECTOR)
+        snapped_points_layer = snapped_points_vector.GetLayer()
+
+        # should snap to stream point [4, 3] in the array above
+        # if not considering flow accumulation, it would snap to the
+        # nearest stream point found first in the array, at [2, 1]
+        points = [
+            shapely.wkb.loads(bytes(feature.GetGeometryRef().ExportToWkb()))
+            for feature in snapped_points_layer]
+        self.assertEqual(len(points), 1)
+        self.assertEqual((points[0].x, points[0].y), (9, -11))
+
     def test_check_geometries(self):
         """DelineateIt: Check that we can reasonably repair geometries."""
         from natcap.invest.delineateit import delineateit