Repository: natcap/invest (python)

Issue:
HRA regression when rasterizing spatial criteria
A [user on the forums](https://community.naturalcapitalproject.org/t/habitat-risk-assessment-problems-with-risk-calculation-and-outputs/3122/15) noticed a difference in how spatial criteria were being represented between InVEST 3.9 (former HRA) and InVEST 3.12 (the new HRA implementation), specifically around how input vectors are being rasterized.

After digging into this for a while, it appears that the InVEST 3.9 implementation makes the following distinction that I missed when rewriting HRA:
* Habitat and stressor vectors are rasterized with `ALL_TOUCHED=TRUE`
* Spatial criteria vectors are rasterized with `ALL_TOUGHED=FALSE`

The current HRA does not make this distinction, but probably should, at least for backwards compatibility.

For what it's worth, I could not find this distinction in the InVEST 3.3.0 HRA.

Pull request:
Fix rasterization of spatial criteria relative to InVEST 3.9
This PR restores a bit of functionality from the InVEST 3.9 HRA implementation that I had missed in the rewrite.  In InVEST 3.9, habitat and stressor vectors are rasterized with `ALL_TOUCHED=TRUE`, while spatial criteria vectors are rasterized with `ALL_TOUCHED=FALSE`.  The HRA rewrite was rasterizing everything as `ALL_TOUCHED=TRUE`, leading to the behavior discussed in https://community.naturalcapitalproject.org/t/habitat-risk-assessment-problems-with-risk-calculation-and-outputs/3122/14.

The HRA 3.9 behavior for spatial criteria seems reasonable to me, so I've restored it here and made sure that it's being tested.

No changes in the user's guide are needed for this since we're just restoring existing functionality.

Fixes #1120 

## Checklist
- [x] Updated HISTORY.rst and link to any relevant issue (if these changes are user-facing)
~~- [ ] Updated the user's guide (if needed)~~
~~- [ ] Tested the affected models' UIs (if relevant)~~


Patch:
diff --git a/HISTORY.rst b/HISTORY.rst
index 6af58ab6e..0b26efdb2 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -45,6 +45,12 @@ Unreleased Changes
       now reprojected to the ``lulc_cur_path`` raster. This fixes a bug where
       rasters with a different SRS would appear to not intersect the
       ``lulc_cur_path`` even if they did. (https://github.com/natcap/invest/issues/1093)
+* HRA
+    * Fixed a regression relative to InVEST 3.9.0 outputs where spatial
+      criteria vectors were being rasterized with the ``ALL_TOUCHED=TRUE``
+      flag, leading to a perceived buffering of spatial criteria in certain
+      cases.  In InVEST 3.9.0, these were rasterized with ``ALL_TOUCHED=FALSE``.
+      https://github.com/natcap/invest/issues/1120
 * Workbench
     * Fixed a bug where the Workbench would become unresponsive during an
       InVEST model run if the model emitted a very high volume of log messages.
diff --git a/src/natcap/invest/hra.py b/src/natcap/invest/hra.py
index 31aed8da8..619f1ff00 100644
--- a/src/natcap/invest/hra.py
+++ b/src/natcap/invest/hra.py
@@ -361,6 +361,7 @@ def execute(args):
     alignment_dependent_tasks = []
     aligned_habitat_raster_paths = {}
     aligned_stressor_raster_paths = {}
+    habitat_stressor_vectors = set([])
     for name, attributes in itertools.chain(habitats_info.items(),
                                             stressors_info.items(),
                                             spatial_criteria_attrs.items()):
@@ -415,6 +416,10 @@ def execute(args):
         # If the input is a vector, reproject to the AOI SRS and simplify.
         # Rasterization happens in the alignment step.
         elif gis_type == pygeoprocessing.VECTOR_TYPE:
+            # Habitats and stressors are rasterized with ALL_TOUCHED=TRUE
+            if name in habitats_info or name in stressors_info:
+                habitat_stressor_vectors.add(source_filepath)
+
             # Using Shapefile here because its driver appears to not raise a
             # warning if a MultiPolygon geometry is inserted into a Polygon
             # layer, which was happening on a real-world sample dataset while
@@ -472,6 +477,7 @@ def execute(args):
             'vector_path_map': alignment_source_vector_paths,
             'target_pixel_size': (resolution, -resolution),
             'target_srs_wkt': target_srs_wkt,
+            'all_touched_vectors': habitat_stressor_vectors,
         },
         task_name='Align raster stack',
         target_path_list=list(itertools.chain(
@@ -1395,7 +1401,7 @@ def _rasterize_aoi_regions(source_aoi_vector_path, base_raster_path,
 
 
 def _align(raster_path_map, vector_path_map, target_pixel_size,
-           target_srs_wkt):
+           target_srs_wkt, all_touched_vectors=None):
     """Align a stack of rasters and/or vectors.
 
     In HRA, habitats and stressors (and optionally criteria) may be defined as
@@ -1429,6 +1435,9 @@ def _align(raster_path_map, vector_path_map, target_pixel_size,
         target_pixel_size (tuple): The pixel size of the target rasters, in
             the form (x, y), expressed in meters.
         target_srs_wkt (string): The target SRS of the aligned rasters.
+        all_touched_vectors=None (set): A set of vector paths found in
+            ``vector_path_map`` that should be rasterized with
+            ``ALL_TOUCHED=TRUE``.
 
     Returns:
         ``None``
@@ -1485,6 +1494,8 @@ def _align(raster_path_map, vector_path_map, target_pixel_size,
     # that align with the bounding box we determined earlier.
     # This approach yields more precise rasters than resampling an
     # already-rasterized vector through align_and_resize_raster_stack.
+    if all_touched_vectors is None:
+        all_touched_vectors = set([])
     if vector_path_map:
         LOGGER.info(f'Aligning {len(vector_path_map)} vectors')
         for source_vector_path, target_raster_path in vector_path_map.items():
@@ -1495,7 +1506,12 @@ def _align(raster_path_map, vector_path_map, target_pixel_size,
             raster_type = _TARGET_GDAL_TYPE_BYTE
             nodata_value = _TARGET_NODATA_BYTE
             burn_values = [1]
-            rasterize_option_list = ['ALL_TOUCHED=TRUE']
+            rasterize_option_list = []
+
+            # Only rasterize with ALL_TOUCHED=TRUE if we know it should be
+            # rasterized as such (habitats/stressors)
+            if source_vector_path in all_touched_vectors:
+                rasterize_option_list.append('ALL_TOUCHED=TRUE')
 
             for field in layer.schema:
                 fieldname = field.GetName()


Test patch:
diff --git a/tests/test_hra.py b/tests/test_hra.py
index 17fd7f845..d41c4877a 100644
--- a/tests/test_hra.py
+++ b/tests/test_hra.py
@@ -608,7 +608,8 @@ class HRAUnitTests(unittest.TestCase):
                 self.workspace_dir, 'aligned_criterion_vector.tif'),
         }
 
-        hra._align(raster_path_map, vector_path_map, (30, -30), SRS_WKT)
+        hra._align(raster_path_map, vector_path_map, (30, -30), SRS_WKT,
+                  all_touched_vectors=set([habitat_vector_path]))
 
         # Calculated by hand given the above spatial inputs and
         # (30, -30) pixels.  All rasters should share the same extents and
@@ -654,15 +655,16 @@ class HRAUnitTests(unittest.TestCase):
 
         # The aligned criterion raster should have been rasterized from the
         # rating column.
+        # This is an ALL_TOUCHED=FALSE rasterization.
         ndta = hra._TARGET_NODATA_FLOAT32
         expected_criterion_array = numpy.array([
-            [0.12, 0.12, 0.12, 0.12, 0.12, 0.12, ndta, ndta, ndta, ndta],
-            [0.12, 0.12, 0.12, 0.12, 0.12, 0.12, 0.12, ndta, ndta, ndta],
-            [0.12, 0.12, 0.12, 0.12, 0.12, 0.12, 0.12, ndta, ndta, ndta],
-            [0.12, 0.12, 0.12, 0.12, 0.12, 0.12, 0.12, ndta, ndta, ndta],
+            [ndta, ndta, 0.12, 0.12, 0.12, ndta, ndta, ndta, ndta, ndta],
+            [ndta, 0.12, 0.12, 0.12, 0.12, 0.12, ndta, ndta, ndta, ndta],
             [0.12, 0.12, 0.12, 0.12, 0.12, 0.12, 0.12, ndta, ndta, ndta],
             [0.12, 0.12, 0.12, 0.12, 0.12, 0.12, 0.12, ndta, ndta, ndta],
+            [0.12, 0.12, 0.12, 0.12, 0.12, 0.12, ndta, ndta, ndta, ndta],
             [ndta, 0.12, 0.12, 0.12, 0.12, 0.12, ndta, ndta, ndta, ndta],
+            [ndta, ndta, 0.12, 0.12, ndta, ndta, ndta, ndta, ndta, ndta],
             [ndta, ndta, ndta, ndta, ndta, ndta, ndta, ndta, ndta, ndta],
             [ndta, ndta, ndta, ndta, ndta, ndta, ndta, ndta, ndta, ndta]],
             dtype=numpy.float32)