Repository: natcap/invest (python)

Issue:
Habitat Quality Rarity output has confusing negative values
@newtpatrol pointed out that when a current / future LULC has more of a habitat than the baseline lulc, negative numbers are reported. Rarity is currently defined in source as:
```python
numerator = lulc_code_count_x[code] * lulc_area
denominator = lulc_code_count_b[code] * base_area
ratio = 1.0 - (numerator / denominator)
code_index[code] = ratio
```
But I wonder if denominator should really be:

`denominator = (lulc_code_count_b[code] * base_area) + numerator`

This would give values between 0 and 1, where a high numerator (not rare) would give a high ratio, and thus low rarity value. A low numerator would give a low ratio, and thus high rarity value. A current / future pixel area which is the same as the baseline would give a rarity of 0.5. So above 0.5, more rare (less abundant), less than 0.5 less rare (more abundant).

OR, we could simply cap the ratio at 0.0, so we don't see negative numbers. Then, everything above 0.0 would be more rare than the baseline and 0.0 would indicate equal abundance or more abundance.

Either way, it seems desirable that this should be a 0-1 output and better defined in the Users Guide.

- [x] Update Users Guide

Pull request description:
Update HQ rarity index calculation
# Description
This PR changes the output of habitat rarity by constraining the values to a 0-1 index. One thing to note is that 0.5 indicates no abundance change / no rarity change. Values less than 0.5 indicate more abundance and "less" rare. Values greater than 0.5 indicate less abundance and higher rarity. See issue #720 for a complete description.

Fixes #720

Users Guide update is in a PR https://github.com/natcap/invest.users-guide/pull/51 

# 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 09ebdc8df..ce97163ec 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -34,9 +34,17 @@
 .. :changelog:
 
 
-..
-  Unreleased Changes
-  ------------------
+Unreleased Changes
+------------------
+* 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
+      negative values. Now values of 0 indicate current/future LULC not
+      represented in baseline LULC; values 0 to 0.5 indicate more
+      abundance in current/future LULC and therefore less rarity; values
+      of 0.5 indicate same abundance between baseline and current/future
+      LULC; values 0.5 to 1 indicate less abundance in current/future LULC
+      and therefore higher rarity.
 
 3.9.2 (2021-10-29)
 ------------------
diff --git a/src/natcap/invest/habitat_quality.py b/src/natcap/invest/habitat_quality.py
index 82658f7ce..2c9ef84c0 100644
--- a/src/natcap/invest/habitat_quality.py
+++ b/src/natcap/invest/habitat_quality.py
@@ -755,6 +755,12 @@ def _compute_rarity_operation(
         base_lulc_path_band, lulc_path_band, new_cover_path, rarity_path):
     """Calculate habitat rarity.
 
+    Output rarity values will be an index from 0 - 1 where:
+       pixel > 0.5 - more rare
+       pixel < 0.5 - less rare
+       pixel = 0.5 - no rarity change
+       pixel = 0.0 - LULC not found in the baseline for comparison
+
     Args:
         base_lulc_path_band (tuple): a 2 tuple for the path to input base
             LULC raster of the form (path, band index).
@@ -819,7 +825,7 @@ def _compute_rarity_operation(
         if code in lulc_code_count_b:
             numerator = lulc_code_count_x[code] * lulc_area
             denominator = lulc_code_count_b[code] * base_area
-            ratio = 1.0 - (numerator / denominator)
+            ratio = 1.0 - (numerator / (denominator + numerator))
             code_index[code] = ratio
         else:
             code_index[code] = 0.0
diff --git a/src/natcap/invest/sdr/sdr.py b/src/natcap/invest/sdr/sdr.py
index 5e324e5d9..52a4d5c5e 100644
--- a/src/natcap/invest/sdr/sdr.py
+++ b/src/natcap/invest/sdr/sdr.py
@@ -789,17 +789,15 @@ def _calculate_ls_factor(
 def _calculate_rkls(
         ls_factor_path, erosivity_path, erodibility_path, stream_path,
         rkls_path):
-    """Calculate potential soil loss (tons / (pixel * year)) using RKLS.
+    """Calculate per-pixel potential soil loss using the RKLS.
 
     (revised universal soil loss equation with no C or P).
 
     Args:
         ls_factor_path (string): path to LS raster that has square pixels in
             meter units.
-        erosivity_path (string): path to erosivity raster
-            (MJ * mm / (ha * hr * yr))
+        erosivity_path (string): path to per pixel erosivity raster
         erodibility_path (string): path to erodibility raster
-            (t * ha * hr / (MJ * ha * mm))
         stream_path (string): path to drainage raster
             (1 is drainage, 0 is not)
         rkls_path (string): path to RKLS raster
@@ -817,21 +815,22 @@ def _calculate_rkls(
 
     cell_size = abs(
         pygeoprocessing.get_raster_info(ls_factor_path)['pixel_size'][0])
-    cell_area_ha = cell_size**2 / 10000.0  # hectares per pixel
+    cell_area_ha = cell_size**2 / 10000.0
 
     def rkls_function(ls_factor, erosivity, erodibility, stream):
         """Calculate the RKLS equation.
 
         Args:
-            ls_factor (numpy.ndarray): length/slope factor. unitless.
-            erosivity (numpy.ndarray): related to peak rainfall events. units:
-                MJ * mm / (ha * hr * yr)
-            erodibility (numpy.ndarray): related to the potential for soil to
-                erode. units: t * ha * hr / (MJ * ha * mm)
-            stream (numpy.ndarray): stream mask (1 stream, 0 no stream)
+            ls_factor (numpy.ndarray): length/slope factor
+        erosivity (numpy.ndarray): related to peak rainfall events
+        erodibility (numpy.ndarray): related to the potential for soil to
+            erode
+        stream (numpy.ndarray): stream mask (1 stream, 0 no stream)
 
         Returns:
-            numpy.ndarray of RKLS values in tons / (pixel * year))
+            ls_factor * erosivity * erodibility * usle_c_p * avg_aspect or
+            nodata if any values are nodata themselves.
+
         """
         rkls = numpy.empty(ls_factor.shape, dtype=numpy.float32)
         nodata_mask = (
@@ -844,11 +843,9 @@ def _calculate_rkls(
         valid_mask = nodata_mask & (stream == 0)
         rkls[:] = _TARGET_NODATA
 
-        rkls[valid_mask] = (           # rkls units are tons / (pixel * year)
-            ls_factor[valid_mask] *    # unitless
-            erosivity[valid_mask] *    # MJ * mm / (ha * hr * yr)
-            erodibility[valid_mask] *  # t * ha * hr / (MJ * ha * mm)
-            cell_area_ha)              # ha / pixel
+        rkls[valid_mask] = (
+            ls_factor[valid_mask] * erosivity[valid_mask] *
+            erodibility[valid_mask] * cell_area_ha)
 
         # rkls is 1 on the stream
         rkls[nodata_mask & (stream == 1)] = 1


Test patch:
diff --git a/tests/test_habitat_quality.py b/tests/test_habitat_quality.py
index 55afd7849..bbbdeb94b 100644
--- a/tests/test_habitat_quality.py
+++ b/tests/test_habitat_quality.py
@@ -233,8 +233,8 @@ class HabitatQualityTests(unittest.TestCase):
                 'deg_sum_f_regression.tif': 33.931896,
                 'quality_c_regression.tif': 7499.983,
                 'quality_f_regression.tif': 4999.9893,
-                'rarity_c_regression.tif': 2500.0000000,
-                'rarity_f_regression.tif': 2500.0000000}.items():
+                'rarity_c_regression.tif': 3333.3335,
+                'rarity_f_regression.tif': 3333.3335}.items():
             raster_path = os.path.join(args['workspace_dir'], output_filename)
             # Check that the raster's computed values are what we expect.
             # In this case, the LULC and threat rasters should have been
@@ -297,8 +297,8 @@ class HabitatQualityTests(unittest.TestCase):
                 'deg_sum_f_regression.tif': 46.279358,
                 'quality_c_regression.tif': 7499.9414,
                 'quality_f_regression.tif': 4999.955,
-                'rarity_c_regression.tif': 2500.0000000,
-                'rarity_f_regression.tif': 2500.0000000}.items():
+                'rarity_c_regression.tif': 3333.3335,
+                'rarity_f_regression.tif': 3333.3335}.items():
             raster_path = os.path.join(args['workspace_dir'], output_filename)
             # Check that the raster's computed values are what we expect.
             # In this case, the LULC and threat rasters should have been
@@ -358,8 +358,8 @@ class HabitatQualityTests(unittest.TestCase):
                 'deg_sum_f_regression.tif': 46.279358,
                 'quality_c_regression.tif': 7499.9414,
                 'quality_f_regression.tif': 4999.955,
-                'rarity_c_regression.tif': 2500.0000000,
-                'rarity_f_regression.tif': 2500.0000000}.items():
+                'rarity_c_regression.tif': 3333.3335,
+                'rarity_f_regression.tif': 3333.3335}.items():
             raster_path = os.path.join(args['workspace_dir'], output_filename)
             # Check that the raster's computed values are what we expect.
             # In this case, the LULC and threat rasters should have been
@@ -420,8 +420,8 @@ class HabitatQualityTests(unittest.TestCase):
                 'deg_sum_f_regression.tif': 46.279358,
                 'quality_c_regression.tif': 7499.9414,
                 'quality_f_regression.tif': 4999.955,
-                'rarity_c_regression.tif': 2500.0000000,
-                'rarity_f_regression.tif': 2500.0000000}.items():
+                'rarity_c_regression.tif': 3333.3335,
+                'rarity_f_regression.tif': 3333.3335}.items():
             raster_path = os.path.join(args['workspace_dir'], output_filename)
             # Check that the raster's computed values are what we expect.
             # In this case, the LULC and threat rasters should have been
@@ -563,8 +563,8 @@ class HabitatQualityTests(unittest.TestCase):
                 'deg_sum_f_regression.tif': 46.279358,
                 'quality_c_regression.tif': 7499.9414,
                 'quality_f_regression.tif': 4999.955,
-                'rarity_c_regression.tif': 2500.0000000,
-                'rarity_f_regression.tif': 2500.0000000
+                'rarity_c_regression.tif': 3333.3335,
+                'rarity_f_regression.tif': 3333.3335
         }.items():
             assert_array_sum(os.path.join(
                 args['workspace_dir'], output_filename), assert_value)
@@ -1078,7 +1078,7 @@ class HabitatQualityTests(unittest.TestCase):
             7499.524)
         assert_array_sum(
             os.path.join(args['workspace_dir'], 'rarity_c.tif'),
-            2500.000000)
+            3333.3335)
 
     def test_habitat_quality_missing_lucodes_in_table(self):
         """Habitat Quality: on missing lucodes in the sensitivity table."""