Repository: natcap/invest
Pull request #577 (2021-07-07 20:27:35)

Pr resolution 575
This PR resolves the conflict in an auto PR. The conflict was the Makefile's sampledata revision. 

https://github.com/natcap/invest/pull/575



Patch:
diff --git a/HISTORY.rst b/HISTORY.rst
index 2d0927d49..900721850 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -104,6 +104,11 @@ Unreleased Changes (3.9.1)
 * SDR
     * Fixed a bug in validation that did not warn against different coordinate
       systems (all SDR inputs must share a common coordinate system).
+    * Fixed a bug that was incorrectly using a factor of 0.0986 rather than
+      0.0896. This would have a minor effect on end-user results.
+    * Changed how SDR thresholds its L factor to allow direct thresholding
+      rather than based off of upstream area. Exposed this parameter as
+      ``l_max`` in the ``args`` input and in the user interface.
 * Wind Energy
     * Raising ValueError when AOI does not intersect Wind Data points.
 
diff --git a/Makefile b/Makefile
index 5734432a5..4b3d803b3 100644
--- a/Makefile
+++ b/Makefile
@@ -2,7 +2,7 @@
 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    := 08e499a8f21c6acf13bcc652e34d5410f9664bc7
+GIT_SAMPLE_DATA_REPO_REV    := f8c3ef11d06cee9d6f5a07f48057c05939e83028
 
 GIT_TEST_DATA_REPO          := https://bitbucket.org/natcap/invest-test-data.git
 GIT_TEST_DATA_REPO_PATH     := $(DATA_DIR)/invest-test-data
diff --git a/src/natcap/invest/sdr/sdr.py b/src/natcap/invest/sdr/sdr.py
index 39eadb47f..97fd8004c 100644
--- a/src/natcap/invest/sdr/sdr.py
+++ b/src/natcap/invest/sdr/sdr.py
@@ -141,6 +141,19 @@ ARGS_SPEC = {
             "about": "Borselli IC0 parameter.",
             "name": "Borselli IC0 Parameter"
         },
+        "l_max": {
+            "type": "number",
+            "expression": "value > 0",
+            "units": u.none,
+            "about": (
+                "Values of L (the slope length component of the LS "
+                "slope length * slope gradient factor) larger than this value "
+                "will be clamped to this value. Ranges of 122-333 (unitless) "
+                "are found in relevant literature such as "
+                "Desmet and Govers, 1996 and Renard et al., 1997 "
+                "(see user's guide)."),
+            "name": "Max L Value",
+        },
         "drainage_path": {
             "type": "raster",
             "bands": {1: {"type": "number", "units": u.none}},
@@ -390,6 +403,7 @@ def execute(args):
             f_reg['flow_accumulation_path'],
             f_reg['slope_path'],
             f_reg['weighted_avg_aspect_path'],
+            float(args['l_max']),
             f_reg['ls_path']),
         target_path_list=[f_reg['ls_path']],
         dependent_task_list=[
@@ -651,8 +665,8 @@ def execute(args):
 
 
 def _calculate_ls_factor(
-        flow_accumulation_path, slope_path, avg_aspect_path,
-        out_ls_prime_factor_path):
+        flow_accumulation_path, slope_path, avg_aspect_path, l_max,
+        target_ls_prime_factor_path):
     """Calculate LS factor.
 
     Calculates a modified LS factor as Equation 3 from "Extension and
@@ -668,7 +682,9 @@ def _calculate_ls_factor(
         slope_path (string): path to slope raster as a percent
         avg_aspect_path (string): The path to to raster of the weighted average
             of aspects based on proportional flow.
-        out_ls_prime_factor_path (string): path to output ls_prime_factor
+        l_max (float): if the calculated value of L exceeds this value
+            it is clamped to this value.
+        target_ls_prime_factor_path (string): path to output ls_prime_factor
             raster
 
     Returns:
@@ -685,13 +701,15 @@ def _calculate_ls_factor(
     cell_size = abs(flow_accumulation_info['pixel_size'][0])
     cell_area = cell_size ** 2
 
-    def ls_factor_function(percent_slope, flow_accumulation, avg_aspect):
+    def ls_factor_function(
+            percent_slope, flow_accumulation, avg_aspect, l_max):
         """Calculate the LS' factor.
 
         Args:
             percent_slope (numpy.ndarray): slope in percent
             flow_accumulation (numpy.ndarray): upstream pixels
             avg_aspect (numpy.ndarray): the weighted average aspect from MFD
+            l_max (float): max L factor, clamp to this value if L exceeds it
 
         Returns:
             ls_factor
@@ -717,7 +735,7 @@ def _calculate_ls_factor(
             16.8 * numpy.sin(slope_in_radians) - 0.5)
 
         beta = (
-            (numpy.sin(slope_in_radians) / 0.0986) /
+            (numpy.sin(slope_in_radians) / 0.0896) /
             (3 * numpy.sin(slope_in_radians)**0.8 + 0.56))
 
         # Set m value via lookup table: Table 1 in
@@ -735,26 +753,24 @@ def _calculate_ls_factor(
             beta[big_slope_mask] / (1 + beta[big_slope_mask]))
         m_exp[~big_slope_mask] = m_table[m_indexes]
 
-        # from McCool paper: "as a final check against excessively long slope
-        # length calculations ... cap of 333m"
-        # from Rafa, this should really be the upstream area capped to
-        # "333^2 m^2" because McCool is 1D
-        contributing_area[contributing_area > 333**2] = 333**2
-
-        ls_prime_factor = (
+        l_factor = (
             ((contributing_area + cell_area)**(m_exp+1) -
              contributing_area ** (m_exp+1)) /
             ((cell_size ** (m_exp + 2)) * (avg_aspect[valid_mask]**m_exp) *
              (22.13**m_exp)))
 
-        result[valid_mask] = ls_prime_factor * slope_factor
+        # threshold L factor to l_max
+        l_factor[l_factor > l_max] = l_max
+
+        result[valid_mask] = l_factor * slope_factor
         return result
 
     # call vectorize datasets to calculate the ls_factor
     pygeoprocessing.raster_calculator(
         [(path, 1) for path in [
-            slope_path, flow_accumulation_path, avg_aspect_path]],
-        ls_factor_function, out_ls_prime_factor_path, gdal.GDT_Float32,
+            slope_path, flow_accumulation_path, avg_aspect_path]] + [
+            (l_max, 'raw')],
+        ls_factor_function, target_ls_prime_factor_path, gdal.GDT_Float32,
         _TARGET_NODATA)
 
 
@@ -1289,9 +1305,13 @@ def _generate_report(
         sed_deposition_path, watershed_results_sdr_path):
     """Create shapefile with USLE, sed export, retention, and deposition."""
     original_datasource = gdal.OpenEx(watersheds_path, gdal.OF_VECTOR)
+    if os.path.exists(watershed_results_sdr_path):
+        LOGGER.warning(f'overwriting results at {watershed_results_sdr_path}')
+        os.remove(watershed_results_sdr_path)
     driver = gdal.GetDriverByName('ESRI Shapefile')
     target_vector = driver.CreateCopy(
         watershed_results_sdr_path, original_datasource)
+
     target_layer = target_vector.GetLayer()
     target_layer.SyncToDisk()
 
diff --git a/src/natcap/invest/ui/sdr.py b/src/natcap/invest/ui/sdr.py
index ddc7f2c06..26bc44733 100644
--- a/src/natcap/invest/ui/sdr.py
+++ b/src/natcap/invest/ui/sdr.py
@@ -123,6 +123,14 @@ class SDR(model.InVESTModel):
             label='Max SDR Value',
             validator=self.validator)
         self.add_input(self.sdr_max)
+        self.l_max = inputs.Text(
+            args_key='l_max',
+            helptext=(
+                'L will not exceed this value. Ranges of 122-333 (unitless) '
+                'are found in relevant literature.'),
+            label='Max L Value',
+            validator=self.validator)
+        self.add_input(self.l_max)
 
     def assemble_args(self):
         args = {


Test patch:
diff --git a/tests/test_sdr.py b/tests/test_sdr.py
index a48521f08..baca238ab 100644
--- a/tests/test_sdr.py
+++ b/tests/test_sdr.py
@@ -26,11 +26,19 @@ def assert_expected_results_in_vector(expected_results, vector_path):
     watershed_results_vector = None
     watershed_results_layer = None
     watershed_results_feature = None
+    incorrect_vals = {}
     for key in expected_results:
         # Using relative tolerance here because results with different
         # orders of magnitude are tested
-        numpy.testing.assert_allclose(
-            actual_results[key], expected_results[key], rtol=0.000001, atol=0)
+        try:
+            numpy.testing.assert_allclose(
+                actual_results[key], expected_results[key],
+                rtol=0.00001, atol=0)
+        except AssertionError:
+            incorrect_vals[key] = (actual_results[key], expected_results[key])
+    if incorrect_vals:
+        raise AssertionError(
+            f'these key (actual/expected) errors occured: {incorrect_vals}')
 
 
 class SDRTests(unittest.TestCase):
@@ -58,6 +66,7 @@ class SDRTests(unittest.TestCase):
             'k_param': '2',
             'lulc_path': os.path.join(SAMPLE_DATA, 'landuse_90.tif'),
             'sdr_max': '0.8',
+            'l_max': '122',
             'threshold_flow_accumulation': '1000',
             'watersheds_path': os.path.join(SAMPLE_DATA, 'watersheds.shp'),
             'workspace_dir': workspace_dir,
@@ -103,7 +112,7 @@ class SDRTests(unittest.TestCase):
         # use predefined directory so test can clean up files during teardown
         args = {}
         validation_warnings = sdr.validate(args, limit_to=None)
-        self.assertEqual(len(validation_warnings[0][0]), 11)
+        self.assertEqual(len(validation_warnings[0][0]), 12)
 
     def test_sdr_validation_key_no_value(self):
         """SDR test validation that's missing a value on a key."""
@@ -197,24 +206,15 @@ class SDRTests(unittest.TestCase):
         from natcap.invest.sdr import sdr
 
         # use predefined directory so test can clean up files during teardown
-        args = SDRTests.generate_base_args(self.workspace_dir)
+        args = SDRTests.generate_base_args('sdr_test_workspace') #self.workspace_dir)
         # make args explicit that this is a base run of SWY
 
-        gpkg_driver = ogr.GetDriverByName('GPKG')
-        base_vector = ogr.Open(args['watersheds_path'])
-        target_watersheds_path = os.path.join(
-            args['workspace_dir'], 'input_watersheds.gpkg')
-        target_vector = gpkg_driver.CopyDataSource(
-            base_vector, target_watersheds_path)
-        base_vector = None
-        target_vector = None
-        args['watersheds_path'] = target_watersheds_path
         sdr.execute(args)
         expected_results = {
-            'usle_tot': 12.04494380951,
-            'sed_retent': 367660.25,
-            'sed_export': 0.71140885353,
-            'sed_dep': 7.84880876541,
+            'usle_tot': 14.25030517578,
+            'sed_retent': 443994.1875,
+            'sed_export': 0.87300693989,
+            'sed_dep': 9.32623577118,
         }
 
         vector_path = os.path.join(
@@ -262,9 +262,9 @@ class SDRTests(unittest.TestCase):
 
         sdr.execute(args)
         expected_results = {
-            'sed_retent': 367660.25,
-            'sed_export': 0.71140885353,
-            'usle_tot': 12.04494380951,
+            'sed_retent': 443994.1875,
+            'sed_export': 0.87300693989,
+            'usle_tot': 14.25030517578,
         }
 
         vector_path = os.path.join(
@@ -286,9 +286,9 @@ class SDRTests(unittest.TestCase):
         sdr.execute(args)
 
         expected_results = {
-            'sed_retent': 335578.875,
-            'sed_export': 0.60814452171,
-            'usle_tot': 11.43409824371,
+            'sed_retent': 376752.75,
+            'sed_export': 0.6872266531,
+            'usle_tot': 12.6965303421,
         }
 
         vector_path = os.path.join(
@@ -311,9 +311,9 @@ class SDRTests(unittest.TestCase):
         sdr.execute(args)
 
         expected_results = {
-            'sed_retent': 414067.65625,
-            'sed_export': 0.88248616457,
-            'usle_tot': 11.37281990051,
+            'sed_retent': 479059.4375,
+            'sed_export': 1.03590250015,
+            'usle_tot': 12.97211265564,
         }
 
         vector_path = os.path.join(
@@ -353,17 +353,6 @@ class SDRTests(unittest.TestCase):
 
         # use predefined directory so test can clean up files during teardown
         args = SDRTests.generate_base_args(self.workspace_dir)
-        # make args explicit that this is a base run of SWY
-
-        gpkg_driver = ogr.GetDriverByName('GPKG')
-        base_vector = ogr.Open(args['watersheds_path'])
-        target_watersheds_path = os.path.join(
-            args['workspace_dir'], 'input_watersheds.gpkg')
-        target_vector = gpkg_driver.CopyDataSource(
-            base_vector, target_watersheds_path)
-        base_vector = None
-        target_vector = None
-        args['watersheds_path'] = target_watersheds_path
 
         # remove a row from the biophysical table so that lulc value is missing
         bad_biophysical_path = os.path.join(
@@ -379,7 +368,8 @@ class SDRTests(unittest.TestCase):
             sdr.execute(args)
         self.assertTrue(
             "The missing values found in the LULC raster but not the table"
-            " are: [2.]" in str(context.exception))
+            " are: [2.]" in str(context.exception),
+            f'error does not match {str(context.exception)}')
 
     @staticmethod
     def _assert_regression_results_equal(
@@ -397,10 +387,10 @@ class SDRTests(unittest.TestCase):
                 expected aggregated_results.shp table in the form of
                 fid,vri_sum,qb_val per line
 
-        Returns:
-            None
+        Return:
+            ``None``
 
-        Raises:
+        Raise:
             AssertionError if any files are missing or results are out of
             range by `tolerance_places`
         """
@@ -455,10 +445,10 @@ class SDRTests(unittest.TestCase):
             directory_path (string): a path to a directory whose contents will
                 be checked against the files listed in `base_list_file`
 
-        Returns:
-            None
+        Return:
+            ``None``
 
-        Raises:
+        Raise:
             AssertionError when there are files listed in `base_list_file`
                 that don't exist in the directory indicated by `path`
         """