Repository: natcap/invest (python)

Issue:
HQ exponential decay discrepancy
James mentioned that there is a discrepancy in how exponential decay is handled between source and the UG.

What the source code is doing: `exp(-dist / max_dist)`
What the UG says it should be doing: `exp(-(2.99 / max_dist) * dist)`

This came up from https://github.com/natcap/pygeoprocessing/issues/268#issuecomment-1299358823.

Pull request:
Decay habitat quality threats using euclidean distance transform
## Description
The Habitat Quality model has used convolutions as an implementation for decaying threat rasters over distance since InVEST 3.3.0. This approach strayed from the implementation described in the User's Guide and the User's Guide was never updated to reflect it. When thinking about updating the User's Guide to mention the convolution implementation it was not clear that the exponential decay via convolutions was providing the desired result. The degradation raster produced from the filtered threat rasters did not make intuitive sense either. Thus a decision was made to fall back to a euclidean distance transform implementation.

- Tests have been "eye" tested and updated.
- Model has been run using the Workbench UI and spot checked in GIS
- User's Guide has a PR here: https://github.com/natcap/invest.users-guide/pull/109

Fixes #1104 

## Checklist
- [x] Updated HISTORY.rst and link to any relevant issue (if these changes are user-facing)
- [x] Updated the user's guide (if needed)
- [x] Tested the affected models' UIs (if relevant)


Patch:
diff --git a/HISTORY.rst b/HISTORY.rst
index b6111f4be..1c8d8481d 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -43,6 +43,17 @@ Unreleased Changes
       (`#1151 <https://github.com/natcap/invest/issues/1151>`_)
     * The ``invest`` command-line-interface no longer opens a graphical
       interface to InVEST. (`#755 <https://github.com/natcap/invest/issues/755>`_)
+* Habitat Quality
+    * The model now uses an euclidean distance implementation for decaying
+      threat rasters both linearly and exponentially. Since InVEST 3.3.0 a
+      convolution implementation has been used, which reflected how
+      the density of a threat or surrounding threat pixels could have an
+      even greater, cumulative impact and degradation over space. However, this
+      was never properly documented in the User's Guide and is not the approach
+      taken in the publication. The convolution implementation also produced
+      degradation and quality outputs that were difficult to interpret.
+    * There should be a noticeable runtime improvement from calculating
+      euclidean distances vs convolutions.
 * 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>`_)
diff --git a/src/natcap/invest/habitat_quality.py b/src/natcap/invest/habitat_quality.py
index 0eb644bde..d7aa15118 100644
--- a/src/natcap/invest/habitat_quality.py
+++ b/src/natcap/invest/habitat_quality.py
@@ -5,10 +5,9 @@ import logging
 import os
 
 import numpy
+from osgeo import gdal
 import pygeoprocessing
 import taskgraph
-from osgeo import gdal
-from osgeo import osr
 
 from .model_metadata import MODEL_METADATA
 from . import spec_utils
@@ -265,8 +264,7 @@ def execute(args):
     output_dir = args['workspace_dir']
     intermediate_output_dir = os.path.join(
         args['workspace_dir'], 'intermediate')
-    kernel_dir = os.path.join(intermediate_output_dir, 'kernels')
-    utils.make_directories([intermediate_output_dir, output_dir, kernel_dir])
+    utils.make_directories([intermediate_output_dir, output_dir])
 
     taskgraph_working_dir = os.path.join(
         intermediate_output_dir, '_taskgraph_working_dir')
@@ -483,7 +481,7 @@ def execute(args):
     for lulc_key, lulc_path in lulc_path_dict.items():
         LOGGER.info(f'Calculating habitat quality for landuse: {lulc_path}')
 
-        threat_convolve_task_list = []
+        threat_decay_task_list = []
         sensitivity_task_list = []
 
         # Create raster of habitat based on habitat field
@@ -535,36 +533,32 @@ def execute(args):
                     f"The max distance for threat: '{threat}' is less than"
                     " or equal to 0. MAX_DIST should be a positive value.")
 
-            kernel_path = os.path.join(
-                kernel_dir, f'kernel_{threat}{lulc_key}{file_suffix}.tif')
-
-            decay_type = threat_data['decay']
+            distance_raster_path = os.path.join(
+                intermediate_output_dir,
+                f'{threat}_distance_transform{lulc_key}{file_suffix}.tif')
 
-            create_kernel_task = task_graph.add_task(
-                func=_create_decay_kernel,
-                args=((threat_raster_path, 1), kernel_path, decay_type,
-                      threat_data['max_dist']),
-                target_path_list=[kernel_path],
+            dist_edt_task = task_graph.add_task(
+                func=pygeoprocessing.distance_transform_edt,
+                args=((threat_raster_path, 1), distance_raster_path),
+                target_path_list=[distance_raster_path],
                 dependent_task_list=[align_task],
-                task_name=f'decay_kernel_{decay_type}{lulc_key}_{threat}')
+                task_name=f'distance edt {lulc_key} {threat}')
+
+            decay_type = threat_data['decay']
 
             filtered_threat_raster_path = os.path.join(
                 intermediate_output_dir,
-                f'filtered_{threat}{lulc_key}{file_suffix}.tif')
-
-            convolve_task = task_graph.add_task(
-                func=pygeoprocessing.convolve_2d,
-                args=((threat_raster_path, 1), (kernel_path, 1),
-                      filtered_threat_raster_path),
-                kwargs={
-                    'target_nodata': _OUT_NODATA,
-                    'ignore_nodata_and_edges': False,
-                    'mask_nodata': False
-                },
+                f'filtered_{decay_type}_{threat}{lulc_key}{file_suffix}.tif')
+
+            dist_decay_task = task_graph.add_task(
+                func=_decay_distance,
+                args=(
+                    distance_raster_path, threat_data['max_dist'],
+                    decay_type, filtered_threat_raster_path),
                 target_path_list=[filtered_threat_raster_path],
-                dependent_task_list=[create_kernel_task],
-                task_name=f'convolve_{decay_type}{lulc_key}_{threat}')
-            threat_convolve_task_list.append(convolve_task)
+                dependent_task_list=[dist_edt_task],
+                task_name=f'distance decay {lulc_key} {threat}')
+            threat_decay_task_list.append(dist_decay_task)
 
             # create sensitivity raster based on threat
             sens_raster_path = os.path.join(
@@ -622,7 +616,7 @@ def execute(args):
             args=(deg_raster_list, deg_sum_raster_path, weight_list),
             target_path_list=[deg_sum_raster_path],
             dependent_task_list=[
-                *threat_convolve_task_list, *sensitivity_task_list,
+                *threat_decay_task_list, *sensitivity_task_list,
                 *access_task_list],
             task_name=f'tot_degradation_{decay_type}{lulc_key}_{threat}')
 
@@ -862,48 +856,6 @@ def _compute_rarity_operation(
                 f' {os.path.basename(lulc_path_band[0])} land cover.')
 
 
-def _create_decay_kernel(raster_path_band, kernel_path, decay_type, max_dist):
-    """Create a decay kernel as a raster.
-
-    Args:
-        raster_path_band (tuple): a 2 tuple of the form
-            (filepath to raster, band index) for raster to decay.
-        kernel_path (string): path to output kernel raster.
-        decay_type (string): type of decay kernel to create, either
-            'linear' | 'exponentional'.
-        max_dist (float): max distance of threat in KM.
-
-    Returns:
-        None
-    """
-    # need the pixel size for the raster so we can create an appropriate
-    # kernel for convolution
-    threat_pixel_size = pygeoprocessing.get_raster_info(
-        raster_path_band[0])['pixel_size']
-
-    # convert max distance (given in KM) to meters
-    max_dist_m = max_dist * 1000.0
-
-    # convert max distance from meters to the number of pixels that
-    # represents on the raster
-    max_dist_pixel = max_dist_m / abs(threat_pixel_size[0])
-    LOGGER.debug(f'Max distance in pixels: {max_dist_pixel}')
-
-    # blur the raster based on the decay type
-    if decay_type == 'linear':
-        decay_func = _make_linear_decay_kernel_path
-    elif decay_type == 'exponential':
-        decay_func = utils.exponential_decay_kernel_raster
-    else:
-        raise ValueError(
-            "Unknown type of decay in biophysical table, should be"
-            f" either 'linear' or 'exponential'. Input was {decay_type}"
-            " for threat"
-            f" {os.path.splitext(os.path.basename(raster_path_band[0]))[0]}")
-
-    decay_func(max_dist_pixel, kernel_path)
-
-
 def _raster_pixel_count(raster_path_band):
     """Count unique pixel values in single band raster.
 
@@ -928,60 +880,6 @@ def _raster_pixel_count(raster_path_band):
     return counts
 
 
-def _make_linear_decay_kernel_path(max_distance, kernel_path):
-    """Create a linear decay kernel as a raster.
-
-    Pixels in raster are equal to d / max_distance where d is the distance to
-    the center of the raster in number of pixels.
-
-    Args:
-        max_distance (int): number of pixels out until the decay is 0.
-        kernel_path (string): path to output raster whose values are in (0,1)
-            representing distance to edge.
-            Size is (``max_distance`` * 2 + 1)
-
-    Returns:
-        None
-    """
-    kernel_size = int(numpy.round(max_distance * 2 + 1))
-
-    driver = gdal.GetDriverByName('GTiff')
-    kernel_dataset = driver.Create(
-        kernel_path.encode('utf-8'), kernel_size, kernel_size, 1,
-        gdal.GDT_Float32, options=[
-            'BIGTIFF=IF_SAFER', 'TILED=YES', 'BLOCKXSIZE=256',
-            'BLOCKYSIZE=256'])
-
-    # Make some kind of geotransform, it doesn't matter what but
-    # will make GIS libraries behave better if it's all defined
-    kernel_dataset.SetGeoTransform([444720, 30, 0, 3751320, 0, -30])
-    srs = osr.SpatialReference()
-    srs.SetUTM(11, 1)
-    srs.SetWellKnownGeogCS('NAD27')
-    kernel_dataset.SetProjection(srs.ExportToWkt())
-
-    kernel_band = kernel_dataset.GetRasterBand(1)
-    kernel_band.SetNoDataValue(-9999)
-
-    col_index = numpy.array(range(kernel_size))
-    integration = 0.0
-    for row_index in range(kernel_size):
-        distance_kernel_row = numpy.sqrt(
-            (row_index - max_distance) ** 2 +
-            (col_index - max_distance) ** 2).reshape(1, kernel_size)
-        kernel = numpy.where(
-            distance_kernel_row > max_distance, 0.0,
-            (max_distance - distance_kernel_row) / max_distance)
-        integration += numpy.sum(kernel)
-        kernel_band.WriteArray(kernel, xoff=0, yoff=row_index)
-
-    for row_index in range(kernel_size):
-        kernel_row = kernel_band.ReadAsArray(
-            xoff=0, yoff=row_index, win_xsize=kernel_size, win_ysize=1)
-        kernel_row /= integration
-        kernel_band.WriteArray(kernel_row, 0, row_index)
-
-
 def _raster_values_in_bounds(raster_path_band, lower_bound, upper_bound):
     """Check raster values are between ``lower_bound`` and ``upper_bound``.
 
@@ -1013,6 +911,78 @@ def _raster_values_in_bounds(raster_path_band, lower_bound, upper_bound):
     return values_valid
 
 
+def _decay_distance(dist_raster_path, max_dist, decay_type, target_path):
+    """Apply an exponential or linear decay to a distance transform raster.
+
+    The function will set pixels greater than ``max_dist`` to 0.
+
+    Args:
+        dist_raster_path (string): a filepath for the raster to decay.
+            The raster is expected to be a euclidean distance transform with
+            values measuring distance in pixels.
+        max_dist (float): max distance of threat in KM.
+        decay_type (string): a string defining which decay method to use.
+            Options include: 'linear' | 'exponential'.
+        target_path (string): a filepath for a float output raster.
+
+    Returns:
+        None
+    """
+    # get raster pixel size to determine how many pixels max_dist covers
+    threat_pixel_size = pygeoprocessing.get_raster_info(
+        dist_raster_path)['pixel_size']
+
+    # convert max distance (given in KM) to meters
+    max_dist_m = max_dist * 1000.0
+
+    # convert max distance from meters to the number of pixels that
+    # represents on the raster
+    max_dist_pixel = max_dist_m / abs(threat_pixel_size[0])
+    LOGGER.debug(f'Max distance in pixels: {max_dist_pixel}')
+
+    def linear_op(dist):
+        """Linear decay operation."""
+        valid_mask = ~utils.array_equals_nodata(dist, _OUT_NODATA)
+        result = numpy.empty(dist.shape, dtype=numpy.float32)
+        result[:] = _OUT_NODATA
+
+        result[valid_mask] = numpy.where(
+            dist[valid_mask] > max_dist_pixel, 0.0,
+            (max_dist_pixel - dist[valid_mask]) / max_dist_pixel)
+        return result
+
+    def exp_op(dist):
+        """Exponential decay operation."""
+        valid_mask = ~utils.array_equals_nodata(dist, _OUT_NODATA)
+        result = numpy.empty(dist.shape, dtype=numpy.float32)
+        result[:] = _OUT_NODATA
+
+        # Some background on where the 2.99 constant comes from:
+        # With the constant of 2.99, the impact of the threat is reduced by
+        # 95% (to 5%) at the specified max threat distance. So I suspect it's
+        # based on the traditional 95% cutoff that is used in statistics. We
+        # could tweak this cutoff (e.g., 99% decay at max distance), if we
+        # wanted. - Lisa Mandle
+        result[valid_mask] = numpy.where(
+            dist[valid_mask] > max_dist_pixel, 0.0,
+            numpy.exp((-dist[valid_mask] * 2.99) / max_dist_pixel))
+        return result
+
+    if decay_type == 'linear':
+        decay_op = linear_op
+    elif decay_type == 'exponential':
+        decay_op = exp_op
+    else:
+        raise ValueError(
+            "Unknown type of decay in threat table, should be"
+            f" either 'linear' or 'exponential'. Input was '{decay_type}' for"
+            f" output raster path : '{target_path}'")
+
+    pygeoprocessing.raster_calculator(
+        [(dist_raster_path, 1)], decay_op, target_path, gdal.GDT_Float32,
+        _OUT_NODATA)
+
+
 def _validate_threat_path(threat_path, lulc_key):
     """Check ``threat_path`` is a valid raster file against ``lulc_key``.
 


Test patch:
diff --git a/tests/test_habitat_quality.py b/tests/test_habitat_quality.py
index 1fbe47bf5..2c4d9a4d6 100644
--- a/tests/test_habitat_quality.py
+++ b/tests/test_habitat_quality.py
@@ -232,10 +232,10 @@ class HabitatQualityTests(unittest.TestCase):
 
         # Assert values were obtained by summing each output raster.
         for output_filename, assert_value in {
-                'deg_sum_c_regression.tif': 27.153614,
-                'deg_sum_f_regression.tif': 46.279358,
-                'quality_c_regression.tif': 7499.9414,
-                'quality_f_regression.tif': 4999.955,
+                'deg_sum_c_regression.tif': 1422.2422,
+                'deg_sum_f_regression.tif': 1438.791,
+                'quality_c_regression.tif': 6961.4536,
+                'quality_f_regression.tif': 4920.1206,
                 'rarity_c_regression.tif': 3333.3335,
                 'rarity_f_regression.tif': 3333.3335}.items():
             raster_path = os.path.join(args['workspace_dir'], output_filename)
@@ -258,7 +258,7 @@ class HabitatQualityTests(unittest.TestCase):
 
         access_vector_path = os.path.join(
             args['workspace_dir'], 'access_samp.shp')
-        make_access_shp('access_vector_path')
+        make_access_shp(access_vector_path)
         # reproject access vector to WGS84 EPSG 4326
         access_wgs84_path = os.path.join(
             args['workspace_dir'], 'access_samp_wgs84.shp')
@@ -266,7 +266,7 @@ class HabitatQualityTests(unittest.TestCase):
         wgs84_srs.ImportFromEPSG(4326)
         target_access_wkt = wgs84_srs.ExportToWkt()
         pygeoprocessing.reproject_vector(
-            'access_vector_path', target_access_wkt, access_wgs84_path)
+            access_vector_path, target_access_wkt, access_wgs84_path)
         args['access_vector_path'] = access_wgs84_path
 
         scenarios = ['_bas_', '_cur_', '_fut_']
@@ -313,10 +313,10 @@ class HabitatQualityTests(unittest.TestCase):
 
         # Assert values were obtained by summing each output raster.
         for output_filename, assert_value in {
-                'deg_sum_c_regression.tif': 18.91135,
-                'deg_sum_f_regression.tif': 33.931896,
-                'quality_c_regression.tif': 7499.983,
-                'quality_f_regression.tif': 4999.9893,
+                'deg_sum_c_regression.tif': 1422.2422,
+                'deg_sum_f_regression.tif': 1438.791,
+                'quality_c_regression.tif': 6961.4536,
+                'quality_f_regression.tif': 4920.1206,
                 'rarity_c_regression.tif': 3333.3335,
                 'rarity_f_regression.tif': 3333.3335}.items():
             raster_path = os.path.join(args['workspace_dir'], output_filename)
@@ -374,10 +374,10 @@ class HabitatQualityTests(unittest.TestCase):
 
         # Assert values were obtained by summing each output raster.
         for output_filename, assert_value in {
-                'deg_sum_c_regression.tif': 27.153614,
-                'deg_sum_f_regression.tif': 46.279358,
-                'quality_c_regression.tif': 7499.9414,
-                'quality_f_regression.tif': 4999.955,
+                'deg_sum_c_regression.tif': 1422.2422,
+                'deg_sum_f_regression.tif': 1438.791,
+                'quality_c_regression.tif': 6961.4536,
+                'quality_f_regression.tif': 4920.1206,
                 'rarity_c_regression.tif': 3333.3335,
                 'rarity_f_regression.tif': 3333.3335}.items():
             raster_path = os.path.join(args['workspace_dir'], output_filename)
@@ -436,10 +436,10 @@ class HabitatQualityTests(unittest.TestCase):
 
         # Assert values were obtained by summing each output raster.
         for output_filename, assert_value in {
-                'deg_sum_c_regression.tif': 27.153614,
-                'deg_sum_f_regression.tif': 46.279358,
-                'quality_c_regression.tif': 7499.9414,
-                'quality_f_regression.tif': 4999.955,
+                'deg_sum_c_regression.tif': 1422.2422,
+                'deg_sum_f_regression.tif': 1438.791,
+                'quality_c_regression.tif': 6961.4536,
+                'quality_f_regression.tif': 4920.1206,
                 'rarity_c_regression.tif': 3333.3335,
                 'rarity_f_regression.tif': 3333.3335}.items():
             raster_path = os.path.join(args['workspace_dir'], output_filename)
@@ -579,10 +579,10 @@ class HabitatQualityTests(unittest.TestCase):
 
         # Assert values were obtained by summing each output raster.
         for output_filename, assert_value in {
-                'deg_sum_c_regression.tif': 27.153614,
-                'deg_sum_f_regression.tif': 46.279358,
-                'quality_c_regression.tif': 7499.9414,
-                'quality_f_regression.tif': 4999.955,
+                'deg_sum_c_regression.tif': 1422.2422,
+                'deg_sum_f_regression.tif': 1438.791,
+                'quality_c_regression.tif': 6961.4536,
+                'quality_f_regression.tif': 4920.1206,
                 'rarity_c_regression.tif': 3333.3335,
                 'rarity_f_regression.tif': 3333.3335
         }.items():
@@ -1002,7 +1002,7 @@ class HabitatQualityTests(unittest.TestCase):
         # Reasonable to just check quality out in this case
         assert_array_sum(
             os.path.join(args['workspace_dir'], 'quality_c.tif'),
-            7499.524)
+            5852.088)
 
     def test_habitat_quality_case_insensitivty(self):
         """Habitat Quality: with table columns that have camel case."""
@@ -1048,9 +1048,10 @@ class HabitatQualityTests(unittest.TestCase):
         habitat_quality.execute(args)
 
         # Reasonable to just check quality out in this case
+        #assert_array_sum(
         assert_array_sum(
             os.path.join(args['workspace_dir'], 'quality_c.tif'),
-            7499.524)
+            5852.088)
 
     def test_habitat_quality_lulc_baseline_current(self):
         """Habitat Quality: on missing future LULC raster."""
@@ -1095,7 +1096,7 @@ class HabitatQualityTests(unittest.TestCase):
         # Reasonable to just check quality out in this case
         assert_array_sum(
             os.path.join(args['workspace_dir'], 'quality_c.tif'),
-            7499.524)
+            5852.088)
         assert_array_sum(
             os.path.join(args['workspace_dir'], 'rarity_c.tif'),
             3333.3335)