Repository: natcap/invest (python)

Issue:
Patch issues in SDR related to LS factor
* fix beta constant to 0.0896 (reference #532)
* threshold L factor to 122 / user defined parameter
* update test cases due to changes in L threshold and beta
* update user's guide mentioning these changes as well as how MFD `x` is calculated compared to traditional D8.

Pull request:
Bugfix/ls factor 
# Description
Fixed a bug that was incorrectly using a factor of 0.0986 rather than 0.0896 and 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.

Also changed a few stylistic PEP8 issues and one instance of a target raster being called "out".

Fixes #533 

# 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 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 8fb8aaff8..38706cfcb 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    := 598664165b911186e16a987c556f56f993cc78dd
+GIT_SAMPLE_DATA_REPO_REV    := d8c508b64bfcab8364f93535dd3fd737d45593a6
 
 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 e21967e49..6001270c1 100644
--- a/src/natcap/invest/sdr/sdr.py
+++ b/src/natcap/invest/sdr/sdr.py
@@ -152,6 +152,21 @@ ARGS_SPEC = {
             "about": "Borselli IC0 parameter.",
             "name": "Borselli IC0 Parameter"
         },
+        "l_max": {
+            "type": "number",
+            "required": True,
+            "validation_options": {
+                "expression": "value > 0",
+            },
+            "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) "
+                "and ranges are 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",
             "required": False,
@@ -400,6 +415,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=[
@@ -661,8 +677,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
@@ -678,7 +694,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:
@@ -695,13 +713,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
@@ -727,7 +747,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
@@ -745,26 +765,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)
 
 
@@ -1299,9 +1317,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 bb214f902..c1d257a38 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."""
@@ -195,24 +204,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(
@@ -260,9 +260,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(
@@ -284,9 +284,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(
@@ -309,9 +309,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(
@@ -351,17 +351,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(
@@ -377,7 +366,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(
@@ -395,10 +385,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`
         """
@@ -453,10 +443,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`
         """