diff --git a/src/natcap/invest/carbon.py b/src/natcap/invest/carbon.py index 8aa50d359..3b049cfe1 100644 --- a/src/natcap/invest/carbon.py +++ b/src/natcap/invest/carbon.py @@ -190,7 +190,41 @@ "The relative annual increase of the price of carbon. " "Required if Run Valuation model is selected."), "name": gettext("annual price change") - } + }, + "raster_output_units": { + "type": "option_string", + "options": { + "per_pixel": { + "display_name": gettext( + "metric tons (or currency units) per pixel"), + "description": gettext( + "In raster outputs reporting amounts of carbon (files " + "prefixed with 'c', 'tot_c', or 'delta'), each " + "pixel's value will represent metric tons of carbon " + "per pixel. In raster outputs reporting economic " + "value of carbon (files prefixed with 'npv'), each " + "pixel's value will represent currency per pixel.")}, + "per_hectare": { + "display_name": gettext( + "metric tons (or currency units) per hectare"), + "description": gettext( + "In raster outputs reporting amounts of carbon (files " + "prefixed with 'c', 'tot_c', or 'delta'), each " + "pixel's value will represent metric tons of carbon " + "per hectare. In raster outputs reporting economic " + "value of carbon (files prefixed with 'npv'), each " + "pixel's value will represent currency per hectare.")}, + }, + "about": gettext( + "How to report values in raster outputs. In raster outputs " + "reporting amounts of carbon (files prefixed with 'c', " + "'tot_c', or 'delta'), each pixel's value will represent " + "metric tons of carbon per pixel or per hectare. In raster " + "outputs reporting economic value of carbon (files prefixed " + "with 'npv'), each pixel's value will represent currency per " + "pixel or per hectare."), + "name": gettext("raster output units") + }, }, "outputs": { "report.html": { @@ -345,6 +379,8 @@ def execute(args): args['rate_change'] (float): Annual rate of change in price of carbon as a percentage. Used if `args['do_valuation']` is present and True. + args['raster_output_units'] (str): how to report values in raster + outputs. Options are Mg C per pixel (default) or Mg C per hectare. args['n_workers'] (int): (optional) The number of worker processes to use for processing this model. If omitted, computation will take place in the current process. @@ -421,7 +457,8 @@ def execute(args): carbon_map_task = graph.add_task( _generate_carbon_map, args=(args[lulc_key], carbon_pool_by_type, - file_registry[storage_key]), + file_registry[storage_key], + args['raster_output_units']), target_path_list=[file_registry[storage_key]], task_name='carbon_map_%s' % storage_key) storage_path_list.append(file_registry[storage_key]) @@ -499,7 +536,8 @@ def execute(args): + calculate_npv_tasks) _ = graph.add_task( _generate_report, - args=(tifs_to_summarize, args, file_registry), + args=(tifs_to_summarize, args, file_registry, + args['raster_output_units']), target_path_list=[file_registry['html_report']], dependent_task_list=tasks_to_report, task_name='generate_report') @@ -535,22 +573,32 @@ def _accumulate_totals(raster_path): def _generate_carbon_map( - lulc_path, carbon_pool_by_type, out_carbon_stock_path): + lulc_path, carbon_pool_by_type, out_carbon_stock_path, + raster_output_units): """Generate carbon stock raster by mapping LULC values to carbon pools. 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 carbon_pool_by_type (dict): a dictionary that maps landcover values - to carbon storage densities per area (Mg C/Ha). + to carbon storage densities per area (Mg C/ha). + out_carbon_stock_path (string): path to output raster that will report + carbon storage values. + raster_output_units (str): how to report values in output raster. + Options are Mg C per pixel (default) or Mg C per hectare. Returns: None. """ - carbon_stock_by_type = dict([ - (lulcid, stock) - for lulcid, stock in carbon_pool_by_type.items()]) + if (raster_output_units == 'per_hectare'): + carbon_stock_by_type = dict([ + (lulcid, stock) + for lulcid, stock in carbon_pool_by_type.items()]) + else: + 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) + for lulcid, stock in carbon_pool_by_type.items()]) reclass_error_details = { 'raster_name': 'LULC', 'column_name': 'lucode', @@ -613,13 +661,16 @@ def _calculate_npv(delta_carbon_path, valuation_constant, npv_out_path): target_path=npv_out_path) -def _generate_report(raster_file_set, model_args, file_registry): +def _generate_report(raster_file_set, model_args, file_registry, + raster_output_units): """Generate a human readable HTML report of summary stats of model run. Args: raster_file_set (set): paths to rasters that need summary stats. model_args (dict): InVEST argument dictionary. file_registry (dict): file path dictionary for InVEST workspace. + raster_output_units (str): units used to report values in output + rasters. Options are Mg C per pixel (default) or Mg C per hectare. Returns: None. @@ -718,12 +769,15 @@ def _generate_report(raster_file_set, model_args, file_registry): for raster_uri, description, units in report: if raster_uri in raster_file_set: - 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 + if raster_output_units == 'per_hectare': + 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 + else: + summary_stat = _accumulate_totals(raster_uri) report_doc.write( '%s' '%.2f%s%s' % ( diff --git a/tests/test_carbon.py b/tests/test_carbon.py index a59550447..9c75f1243 100644 --- a/tests/test_carbon.py +++ b/tests/test_carbon.py @@ -127,8 +127,8 @@ def tearDown(self): """Override tearDown function to remove temporary directory.""" shutil.rmtree(self.workspace_dir) - def test_carbon_full(self): - """Carbon: full model run.""" + def test_carbon_full_per_pixel(self): + """Carbon: full run with raster_output_units set to per_pixel.""" from natcap.invest import carbon args = { @@ -139,6 +139,78 @@ def test_carbon_full(self): 'lulc_cur_year': 2016, 'lulc_fut_year': 2030, 'discount_rate': -7.1, + 'raster_output_units': 'per_pixel', + 'n_workers': -1, + } + + # Create LULC rasters and pools csv in workspace and add them to args. + lulc_names = ['lulc_cur_path', 'lulc_fut_path', 'lulc_redd_path'] + for fill_val, lulc_name in enumerate(lulc_names, 1): + args[lulc_name] = os.path.join(args['workspace_dir'], + lulc_name + '.tif') + make_simple_raster(args[lulc_name], fill_val, -1) + + args['carbon_pools_path'] = os.path.join(args['workspace_dir'], + 'pools.csv') + make_pools_csv(args['carbon_pools_path']) + + carbon.execute(args) + + # Ensure every pixel has the correct total C value. + # Pixel size is 1 m^2 = 0.0001 ha + # Current: 15 + 10 + 60 + 1 = 86 Mg/ha * 0.0001 ha = 0.0086 Mg + assert_raster_equal_value( + os.path.join(args['workspace_dir'], 'tot_c_cur.tif'), 0.0086) + # Future: 5 + 3 + 20 + 0 = 28 Mg/ha * 0.0001 ha = 0.0028 Mg + assert_raster_equal_value( + os.path.join(args['workspace_dir'], 'tot_c_fut.tif'), 0.0028) + # REDD: 2 + 1 + 5 + 0 = 8 Mg/ha * 0.0001 ha = 0.0008 Mg + assert_raster_equal_value( + os.path.join(args['workspace_dir'], 'tot_c_redd.tif'), 0.0008) + + # Ensure deltas are correct. + assert_raster_equal_value( + os.path.join(args['workspace_dir'], 'delta_cur_fut.tif'), -0.0058) + assert_raster_equal_value( + os.path.join(args['workspace_dir'], 'delta_cur_redd.tif'), -0.0078) + + # Ensure NPV calculations are correct. + # Valuation constant based on provided args is 59.00136. + # Future: 59.00136 * -58 Mg/ha = -3422.079 Mg/ha * 0.0001 ha = -0.3422079 Mg + assert_raster_equal_value( + os.path.join(args['workspace_dir'], 'npv_fut.tif'), -0.3422079) + # REDD: 59.00136 * -78 Mg/ha = -4602.106 Mg/ha * 0.0001 ha = -0.4602106 Mg + assert_raster_equal_value( + os.path.join(args['workspace_dir'], 'npv_redd.tif'), -0.4602106) + + # Ensure aggregate results are correct. + report_path = os.path.join(args['workspace_dir'], 'report.html') + # Raster is 10 x 10; therefore, raster total is as follows: + # (x Mg / px) * (100 px) = (x * 100) Mg + for (stat, expected_value) in [ + ('Total cur', 0.86), + ('Total fut', 0.28), + ('Total redd', 0.08), + ('Change in C for fut', -0.58), + ('Change in C for redd', -0.78), + ('Net present value from cur to fut', -34.22), + ('Net present value from cur to redd', -46.02), + ]: + assert_aggregate_result_equal(report_path, stat, expected_value) + + def test_carbon_full_per_hectare(self): + """Carbon: full run with raster_output_units set to per_hectare.""" + from natcap.invest import carbon + + args = { + 'workspace_dir': self.workspace_dir, + 'do_valuation': True, + 'price_per_metric_ton_of_c': 43.0, + 'rate_change': 2.8, + 'lulc_cur_year': 2016, + 'lulc_fut_year': 2030, + 'discount_rate': -7.1, + 'raster_output_units': 'per_hectare', 'n_workers': -1, } @@ -208,6 +280,7 @@ def test_carbon_zero_rates(self): 'lulc_cur_year': 2016, 'lulc_fut_year': 2030, 'discount_rate': 0.0, + 'raster_output_units': 'per_hectare', 'n_workers': -1, } @@ -245,6 +318,7 @@ def test_carbon_without_redd(self): 'lulc_cur_year': 2016, 'lulc_fut_year': 2030, 'discount_rate': -7.1, + 'raster_output_units': 'per_hectare', 'n_workers': -1, } @@ -272,6 +346,7 @@ def test_carbon_missing_landcover_values(self): args = { 'workspace_dir': self.workspace_dir, 'do_valuation': False, + 'raster_output_units': 'per_hectare', 'n_workers': -1, } @@ -305,6 +380,7 @@ def test_carbon_full_undefined_nodata(self): 'lulc_cur_year': 2016, 'lulc_fut_year': 2030, 'discount_rate': -7.1, + 'raster_output_units': 'per_hectare', 'n_workers': -1, } @@ -341,6 +417,7 @@ def setUp(self): 'workspace_dir', 'lulc_cur_path', 'carbon_pools_path', + 'raster_output_units', ] def tearDown(self): diff --git a/workbench/src/renderer/ui_config.js b/workbench/src/renderer/ui_config.js index b9a7af3d6..c53ffdac8 100644 --- a/workbench/src/renderer/ui_config.js +++ b/workbench/src/renderer/ui_config.js @@ -79,6 +79,7 @@ const UI_SPEC = { 'discount_rate', 'rate_change', ], + ['raster_output_units'], ], enabledFunctions: { lulc_fut_path: isSufficient.bind(null, 'calc_sequestration'),