Skip to content

Commit

Permalink
Carbon: allow user to choose raster output units
Browse files Browse the repository at this point in the history
  • Loading branch information
@emilyanndavis
emilyanndavis committed Jan 29, 2025
1 parent acb76eb commit 3fc47c7
Show file tree
Hide file tree
Showing 3 changed files with 151 additions and 19 deletions.
88 changes: 71 additions & 17 deletions src/natcap/invest/carbon.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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": {
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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])
Expand Down Expand Up @@ -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')
Expand Down Expand Up @@ -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',
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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(
'<tr><td>%s</td><td class="number" data-summary-stat="%s">'
'%.2f</td><td>%s</td><td>%s</td></tr>' % (
Expand Down
81 changes: 79 additions & 2 deletions tests/test_carbon.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 = {
Expand All @@ -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,
}

Expand Down Expand Up @@ -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,
}

Expand Down Expand Up @@ -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,
}

Expand Down Expand Up @@ -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,
}

Expand Down Expand Up @@ -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,
}

Expand Down Expand Up @@ -341,6 +417,7 @@ def setUp(self):
'workspace_dir',
'lulc_cur_path',
'carbon_pools_path',
'raster_output_units',
]

def tearDown(self):
Expand Down
1 change: 1 addition & 0 deletions workbench/src/renderer/ui_config.js
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,6 +79,7 @@ const UI_SPEC = {
'discount_rate',
'rate_change',
],
['raster_output_units'],
],
enabledFunctions: {
lulc_fut_path: isSufficient.bind(null, 'calc_sequestration'),
Expand Down

0 comments on commit 3fc47c7

Please sign in to comment.