Repository: natcap/invest
Pull request #492 (2021-02-12 21:24:24)

Fix carbon valuation when rates are 0
# Description
Fixes a bug in the carbon model valuation that came up in [this forum post](https://community.naturalcapitalproject.org/t/economic-value-raster-not-assigning-correct-values-to-pixels/1699). If rate change and discount rate are both 0, the calculated valuation constant became a rate per year instead of a total value, giving results that are off by a factor of `fut_year - cur_year`.

# Checklist
- [x] Updated HISTORY.rst (if these changes are user-facing)

- [ ] Updated the user's guide (if needed)



Patch:
diff --git a/HISTORY.rst b/HISTORY.rst
index 20bbe32d6..a0be0802b 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -54,6 +54,10 @@ Unreleased Changes (3.9.1)
       requirement version to ``0.10.3``
     * Provide a better validation error message when an overview '.ovr' file
       is input instead of a valid raster.
+* Carbon
+    * Fixed a bug where, if rate change and discount rate were set to 0, the
+      valuation results were in $/year rather than $, too small by a factor of 
+      ``lulc_fut_year - lulc_cur_year``.
 * Fisheries Habitat Scenario Tool
     * Fixed divide-by-zero bug that was causing a RuntimeWarning in the logs.
       This bug did not affect the output.
diff --git a/src/natcap/invest/carbon.py b/src/natcap/invest/carbon.py
index 95e385d65..118eeb709 100644
--- a/src/natcap/invest/carbon.py
+++ b/src/natcap/invest/carbon.py
@@ -497,15 +497,21 @@ def _calculate_valuation_constant(
         a floating point number that can be used to multiply a delta carbon
         storage value by to calculate NPV.
     """
-    n_years = int(lulc_fut_year) - int(lulc_cur_year) - 1
+    n_years = lulc_fut_year - lulc_cur_year
     ratio = (
-        1.0 / ((1 + float(discount_rate) / 100.0) *
-               (1 + float(rate_change) / 100.0)))
-    valuation_constant = (
-        float(price_per_metric_ton_of_c) /
-        (float(lulc_fut_year) - float(lulc_cur_year)))
-    if ratio != 1.0:
-        valuation_constant *= (1.0 - ratio ** (n_years + 1)) / (1.0 - ratio)
+        1 / ((1 + discount_rate / 100) * 
+             (1 + rate_change / 100)))
+    valuation_constant = (price_per_metric_ton_of_c / n_years)
+    # note: the valuation formula in the user's guide uses sum notation.
+    # here it's been simplified to remove the sum using the general rule
+    # sum(r^k) from k=0 to N  ==  (r^(N+1) - 1) / (r - 1)
+    # where N = n_years-1 and r = ratio
+    if ratio == 1:
+        # if ratio == 1, we would divide by zero in the equation below
+        # so use the limit as ratio goes to 1, which is n_years
+        valuation_constant *= n_years
+    else:
+        valuation_constant *= (1 - ratio ** n_years) / (1 - ratio)
     return valuation_constant
 
 


Test patch:
diff --git a/tests/test_carbon.py b/tests/test_carbon.py
index f78716e0b..65f3a2589 100644
--- a/tests/test_carbon.py
+++ b/tests/test_carbon.py
@@ -80,7 +80,9 @@ def make_pools_csv(pools_csv_path):
     with open(pools_csv_path, 'w') as open_table:
         open_table.write('C_above,C_below,C_soil,C_dead,lucode,LULC_Name\n')
         open_table.write('15,10,60,1,1,"lulc code 1"\n')
+        # total change from 1 -> 2: -58 metric tons per hectare
         open_table.write('5,3,20,0,2,"lulc code 2"\n')
+        # total change from 1 -> 3: -78 metric tons per hectare
         open_table.write('2,1,5,0,3,"lulc code 3"\n')
 
 
@@ -161,11 +163,14 @@ class CarbonTests(unittest.TestCase):
         carbon.execute(args)
 
         # Add assertions for npv for future and REDD scenarios.
-        # The npv was calculated based on _calculate_npv in carbon.py.
+        # carbon change from cur to fut: 
+        # -58 Mg/ha * .0001 ha/pixel * 43 $/Mg = -0.2494 $/pixel
         assert_raster_equal_value(
-            os.path.join(args['workspace_dir'], 'npv_fut.tif'), -0.0178143)
+            os.path.join(args['workspace_dir'], 'npv_fut.tif'), -0.2494)
+        # carbon change from cur to redd: 
+        # -78 Mg/ha * .0001 ha/pixel * 43 $/Mg = -0.3354 $/pixel
         assert_raster_equal_value(
-            os.path.join(args['workspace_dir'], 'npv_redd.tif'), -0.0239571)
+            os.path.join(args['workspace_dir'], 'npv_redd.tif'), -0.3354)
 
     def test_carbon_future(self):
         """Carbon: regression testing future scenario."""