Merge pull request #1717 from emilyanndavis/feature/1270-per-hectare-…

…outputs Convert per-pixel values to per-hectare values in raster outputs of multiple models
natcap · Feb 14, 2025 · 7bd4a1a · 7bd4a1a
2 parents 20e98f3 + 16d74e8
commit 7bd4a1a
Show file tree

Hide file tree

Showing 12 changed files with 633 additions and 393 deletions.
diff --git a/HISTORY.rst b/HISTORY.rst
@@ -9,7 +9,7 @@
   - Crop Pollination
   - Crop Production
   - DelineateIt
-  - Forest Carbon Edge Effects
+  - Forest Carbon Edge Effect
   - Globio
   - Habitat Quality
   - HRA
@@ -62,12 +62,24 @@ Unreleased Changes
     * Updated styling of the HTML report generated by the carbon model, for
       visual consistency with the Workbench (`InVEST #1732
       <https://github.com/natcap/invest/issues/1732>`_).
+    * Raster outputs that previously contained per-pixel values (e.g., Mg/px)
+      now contain per-hectare values (e.g., Mg/ha). (`InVEST #1270
+      <https://github.com/natcap/invest/issues/1270>`_).
 * Coastal Blue Carbon
     * The ``code`` column in the model's biophysical table input, as well as
       the ``code`` column in the preprocessor's LULC lookup table input and
       ``carbon_pool_transient_template`` output, have been renamed ``lucode``,
       for consistency with other InVEST models (`InVEST #1249
       <https://github.com/natcap/invest/issues/1249>`_).
+* Crop Production
+    * Raster outputs that previously contained per-pixel values (e.g., Mg/px)
+      now contain per-hectare values (e.g., Mg/ha). This change affects both
+      the Percentile and Regression models (`InVEST #1270
+      <https://github.com/natcap/invest/issues/1270>`_).
+* Forest Carbon Edge Effect
+    * Raster outputs that previously contained per-pixel values (e.g., Mg/px)
+      now contain per-hectare values (e.g., Mg/ha). (`InVEST #1270
+      <https://github.com/natcap/invest/issues/1270>`_).
 * Habitat Quality
     * The ``lulc`` column in the sensitivity table input, and the ``lulc_code``
       column in the rarity table outputs, have been renamed ``lucode``, for
@@ -76,6 +88,13 @@ Unreleased Changes
 * NDR
     * Align rasters to the grid of the DEM raster
       (`#1488 <https://github.com/natcap/invest/issues/1488>`_).
+    * Raster outputs that previously contained per-pixel values (e.g., kg/px)
+      now contain per-hectare values (e.g., kg/ha). (`InVEST #1270
+      <https://github.com/natcap/invest/issues/1270>`_).
+* SDR
+    * Raster outputs that previously contained per-pixel values (e.g., Mg/px)
+      now contain per-hectare values (e.g., Mg/ha). (`InVEST #1270
+      <https://github.com/natcap/invest/issues/1270>`_).
 * Urban Cooling
     * Align rasters to the grid of the LULC raster, rather than the ET0 raster
       (`#1488 <https://github.com/natcap/invest/issues/1488>`_).

diff --git a/src/natcap/invest/carbon.py b/src/natcap/invest/carbon.py
@@ -4,7 +4,6 @@
 import logging
 import os
 import time
-from functools import reduce
 
 from osgeo import gdal
 import numpy
@@ -27,7 +26,7 @@
             "scenario, mapped from the Carbon Pools table to the LULC."),
         "bands": {1: {
             "type": "number",
-            "units": u.metric_ton/u.pixel
+            "units": u.metric_ton/u.hectare
         }}
     } for pool, pool_name in [
         ('above', 'aboveground'),
@@ -154,7 +153,7 @@
             "about": gettext(
                 "The calendar year of the future scenario depicted in the "
                 "future LULC map. Required if Run Valuation model is selected."),
-            "name": f"future LULC year"
+            "name": gettext("future LULC year")
         },
         "do_valuation": {
             "type": "boolean",
@@ -201,54 +200,54 @@
             "about": "Raster showing the amount of carbon stored in each pixel for the current scenario. It is a sum of all of the carbon pools provided by the biophysical table.",
             "bands": {1: {
                 "type": "number",
-                "units": u.metric_ton/u.pixel
+                "units": u.metric_ton/u.hectare
             }}
         },
         "tot_c_fut.tif": {
             "about": "Raster showing the amount of carbon stored in each pixel for the future scenario. It is a sum of all of the carbon pools provided by the biophysical table.",
             "bands": {1: {
                 "type": "number",
-                "units": u.metric_ton/u.pixel
+                "units": u.metric_ton/u.hectare
             }},
             "created_if": "lulc_fut_path"
         },
         "tot_c_redd.tif": {
             "about": "Raster showing the amount of carbon stored in each pixel for the REDD scenario. It is a sum of all of the carbon pools provided by the biophysical table.",
             "bands": {1: {
                 "type": "number",
-                "units": u.metric_ton/u.pixel
+                "units": u.metric_ton/u.hectare
             }},
             "created_if": "lulc_redd_path"
         },
         "delta_cur_fut.tif": {
             "about": "Raster showing the difference in carbon stored between the future landscape and the current landscape. In this map some values may be negative and some positive. Positive values indicate sequestered carbon, negative values indicate carbon that was lost.",
             "bands": {1: {
                 "type": "number",
-                "units": u.metric_ton/u.pixel
+                "units": u.metric_ton/u.hectare
             }},
             "created_if": "lulc_fut_path"
         },
         "delta_cur_redd.tif": {
             "about": "Raster showing the difference in carbon stored between the REDD landscape and the current landscape. In this map some values may be negative and some positive. Positive values indicate sequestered carbon, negative values indicate carbon that was lost.",
             "bands": {1: {
                 "type": "number",
-                "units": u.metric_ton/u.pixel
+                "units": u.metric_ton/u.hectare
             }},
             "created_if": "lulc_redd_path"
         },
         "npv_fut.tif": {
             "about": "Rasters showing the economic value of carbon sequestered between the current and the future landscape dates.",
             "bands": {1: {
                 "type": "number",
-                "units": u.currency/u.pixel
+                "units": u.currency/u.hectare
             }},
             "created_if": "lulc_fut_path"
         },
         "npv_redd.tif": {
             "about": "Rasters showing the economic value of carbon sequestered between the current and the REDD landscape dates.",
             "bands": {1: {
                 "type": "number",
-                "units": u.currency/u.pixel
+                "units": u.currency/u.hectare
             }},
             "created_if": "lulc_redd_path"
         },
@@ -297,6 +296,7 @@
 # -1.0 since carbon stocks are 0 or greater
 _CARBON_NODATA = -1.0
 
+
 def execute(args):
     """Carbon.
 
@@ -541,17 +541,15 @@ def _generate_carbon_map(
     Args:
         lulc_path (string): landcover raster with integer pixels.
         out_carbon_stock_path (string): path to output raster that will have
-            pixels with carbon storage values in them with units of Mg*C
+            pixels with carbon storage values in them with units of Mg/ha.
         carbon_pool_by_type (dict): a dictionary that maps landcover values
             to carbon storage densities per area (Mg C/Ha).
 
     Returns:
         None.
     """
-    lulc_info = pygeoprocessing.get_raster_info(lulc_path)
-    pixel_area = abs(numpy.prod(lulc_info['pixel_size']))
     carbon_stock_by_type = dict([
-        (lulcid, stock * pixel_area / 10**4)
+        (lulcid, stock)
         for lulcid, stock in carbon_pool_by_type.items()])
 
     reclass_error_details = {
@@ -704,14 +702,17 @@ def _generate_report(raster_file_set, model_args, file_registry):
             '<table><thead><tr><th>Description</th><th>Value</th><th>Units'
             '</th><th>Raw File</th></tr></thead><tbody>')
 
+        carbon_units = 'metric tons'
+
         # value lists are [sort priority, description, statistic, units]
         report = [
-            (file_registry['tot_c_cur'], 'Total cur', 'Mg of C'),
-            (file_registry['tot_c_fut'], 'Total fut', 'Mg of C'),
-            (file_registry['tot_c_redd'], 'Total redd', 'Mg of C'),
-            (file_registry['delta_cur_fut'], 'Change in C for fut', 'Mg of C'),
+            (file_registry['tot_c_cur'], 'Total cur', carbon_units),
+            (file_registry['tot_c_fut'], 'Total fut', carbon_units),
+            (file_registry['tot_c_redd'], 'Total redd', carbon_units),
+            (file_registry['delta_cur_fut'], 'Change in C for fut',
+             carbon_units),
             (file_registry['delta_cur_redd'],
-             'Change in C for redd', 'Mg of C'),
+             'Change in C for redd', carbon_units),
             (file_registry['npv_fut'],
              'Net present value from cur to fut', 'currency units'),
             (file_registry['npv_redd'],
@@ -720,11 +721,17 @@ def _generate_report(raster_file_set, model_args, file_registry):
 
         for raster_uri, description, units in report:
             if raster_uri in raster_file_set:
-                summary_stat = _accumulate_totals(raster_uri)
+                total = _accumulate_totals(raster_uri)
+                raster_info = pygeoprocessing.get_raster_info(raster_uri)
+                pixel_area = abs(numpy.prod(raster_info['pixel_size']))
+                # Since each pixel value is in Mg/ha, ``total`` is in (Mg/ha * px) = Mg•px/ha.
+                # Adjusted sum = ([total] Mg•px/ha) * ([pixel_area] m^2 / 1 px) * (1 ha / 10000 m^2) = Mg.
+                summary_stat = total * pixel_area / 10000
                 report_doc.write(
-                    '<tr><td>%s</td><td class="number">%.2f</td><td>%s</td>'
-                    '<td>%s</td></tr>' % (
-                        description, summary_stat, units, raster_uri))
+                    '<tr><td>%s</td><td class="number" data-summary-stat="%s">'
+                    '%.2f</td><td>%s</td><td>%s</td></tr>' % (
+                        description, description, summary_stat, units,
+                        raster_uri))
         report_doc.write('</tbody></table></body></html>')