diff --git a/HISTORY.rst b/HISTORY.rst index db68e1b914..bf685cdb85 100644 --- 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 `_). + * Raster outputs that previously contained per-pixel values (e.g., Mg/px) + now contain per-hectare values (e.g., Mg/ha). (`InVEST #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 `_). +* 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 + `_). +* 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 + `_). * 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 `_). + * Raster outputs that previously contained per-pixel values (e.g., kg/px) + now contain per-hectare values (e.g., kg/ha). (`InVEST #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 + `_). * Urban Cooling * Align rasters to the grid of the LULC raster, rather than the ET0 raster (`#1488 `_). diff --git a/src/natcap/invest/carbon.py b/src/natcap/invest/carbon.py index 9a899558ab..944cda9770 100644 --- 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,14 +200,14 @@ "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" }, @@ -216,7 +215,7 @@ "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" }, @@ -224,7 +223,7 @@ "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" }, @@ -232,7 +231,7 @@ "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" }, @@ -240,7 +239,7 @@ "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" }, @@ -248,7 +247,7 @@ "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): '') + 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( - '' - '' % ( - description, summary_stat, units, raster_uri)) + '' % ( + description, description, summary_stat, units, + raster_uri)) report_doc.write('
DescriptionValueUnits' 'Raw File
%s%.2f%s%s
%s' + '%.2f%s%s
') diff --git a/src/natcap/invest/crop_production_percentile.py b/src/natcap/invest/crop_production_percentile.py index ea4dc057f2..fefc3de14f 100644 --- a/src/natcap/invest/crop_production_percentile.py +++ b/src/natcap/invest/crop_production_percentile.py @@ -761,7 +761,6 @@ def execute(args): args['landcover_raster_path'], interpolated_yield_percentile_raster_path, crop_lucode, - pixel_area_ha, percentile_crop_production_raster_path), target_path_list=[percentile_crop_production_raster_path], dependent_task_list=[ @@ -836,7 +835,7 @@ def execute(args): args=([(args['landcover_raster_path'], 1), (interpolated_observed_yield_raster_path, 1), (observed_yield_nodata, 'raw'), (landcover_nodata, 'raw'), - (crop_lucode, 'raw'), (pixel_area_ha, 'raw')], + (crop_lucode, 'raw')], _mask_observed_yield_op, observed_production_raster_path, gdal.GDT_Float32, observed_yield_nodata), target_path_list=[observed_production_raster_path], @@ -853,7 +852,7 @@ def execute(args): output_dir, 'result_table%s.csv' % file_suffix) crop_names = crop_to_landcover_df.index.to_list() - tabulate_results_task = task_graph.add_task( + _ = task_graph.add_task( func=tabulate_results, args=(nutrient_df, yield_percentile_headers, crop_names, pixel_area_ha, @@ -870,14 +869,14 @@ def execute(args): output_dir, _AGGREGATE_VECTOR_FILE_PATTERN % (file_suffix)) aggregate_results_table_path = os.path.join( output_dir, _AGGREGATE_TABLE_FILE_PATTERN % file_suffix) - aggregate_results_task = task_graph.add_task( + _ = task_graph.add_task( func=aggregate_to_polygons, args=(args['aggregate_polygon_path'], target_aggregate_vector_path, landcover_raster_info['projection_wkt'], crop_names, nutrient_df, - yield_percentile_headers, output_dir, file_suffix, - aggregate_results_table_path), + yield_percentile_headers, pixel_area_ha, output_dir, + file_suffix, aggregate_results_table_path), target_path_list=[target_aggregate_vector_path, aggregate_results_table_path], dependent_task_list=dependent_task_list, @@ -888,14 +887,14 @@ def execute(args): def calculate_crop_production(lulc_path, yield_path, crop_lucode, - pixel_area_ha, target_path): + target_path): """Calculate crop production for a particular crop. The resulting production value is: - nodata, where either the LULC or yield input has nodata - 0, where the LULC does not match the given LULC code - - yield * pixel area, where the given LULC code exists + - yield (in Mg/ha), where the given LULC code exists Args: lulc_path (str): path to a raster of LULC codes @@ -903,7 +902,6 @@ def calculate_crop_production(lulc_path, yield_path, crop_lucode, by ``crop_lucode``, in units per hectare crop_lucode (int): LULC code that identifies the crop of interest in the ``lulc_path`` raster. - pixel_area_ha (number): Pixel area in hectares for both input rasters target_path (str): Path to write the output crop production raster Returns: @@ -911,7 +909,7 @@ def calculate_crop_production(lulc_path, yield_path, crop_lucode, """ pygeoprocessing.raster_map( op=lambda lulc, _yield: numpy.where( - lulc == crop_lucode, _yield * pixel_area_ha, 0), + lulc == crop_lucode, _yield, 0), rasters=[lulc_path, yield_path], target_path=target_path, target_nodata=_NODATA_YIELD) @@ -941,7 +939,7 @@ def _zero_observed_yield_op(observed_yield_array, observed_yield_nodata): def _mask_observed_yield_op( lulc_array, observed_yield_array, observed_yield_nodata, - landcover_nodata, crop_lucode, pixel_area_ha): + landcover_nodata, crop_lucode): """Mask total observed yield to crop lulc type. Args: @@ -950,7 +948,6 @@ def _mask_observed_yield_op( observed_yield_nodata (float): yield raster nodata value landcover_nodata (float): landcover raster nodata value crop_lucode (int): code used to mask in the current crop - pixel_area_ha (float): area of lulc raster cells (hectares) Returns: numpy.ndarray with float values of yields masked to crop_lucode @@ -959,13 +956,13 @@ def _mask_observed_yield_op( result = numpy.empty(lulc_array.shape, dtype=numpy.float32) if landcover_nodata is not None: result[:] = observed_yield_nodata - valid_mask = ~pygeoprocessing.array_equals_nodata(lulc_array, landcover_nodata) + valid_mask = ~pygeoprocessing.array_equals_nodata(lulc_array, + landcover_nodata) result[valid_mask] = 0 else: result[:] = 0 lulc_mask = lulc_array == crop_lucode - result[lulc_mask] = ( - observed_yield_array[lulc_mask] * pixel_area_ha) + result[lulc_mask] = observed_yield_array[lulc_mask] return result @@ -986,7 +983,7 @@ def tabulate_results( landcover_raster_path (string): path to landcover raster landcover_nodata (float): landcover raster nodata value output_dir (string): the file path to the output workspace. - file_suffix (string): string to appened to any output filenames. + file_suffix (string): string to append to any output filenames. target_table_path (string): path to 'result_table.csv' in the output workspace @@ -1007,6 +1004,11 @@ def tabulate_results( for nutrient_id in _EXPECTED_NUTRIENT_TABLE_HEADERS for yield_percentile_id in sorted(yield_percentile_headers) + [ 'yield_observed']] + + # Since pixel values in observed and percentile rasters are Mg/(ha•yr), + # raster sums are (Mg•px)/(ha•yr). Before recording sums in + # production_lookup dictionary, convert to Mg/yr by multiplying by ha/px. + with open(target_table_path, 'w') as result_table: result_table.write( 'crop,area (ha),' + 'production_observed,' + @@ -1038,6 +1040,7 @@ def tabulate_results( production_pixel_count += numpy.count_nonzero( valid_mask & (yield_block > 0)) yield_sum += numpy.sum(yield_block[valid_mask]) + yield_sum *= pixel_area_ha production_area = production_pixel_count * pixel_area_ha production_lookup['observed'] = yield_sum result_table.write(',%f' % production_area) @@ -1055,6 +1058,7 @@ def tabulate_results( yield_sum += numpy.sum( yield_block[~pygeoprocessing.array_equals_nodata( yield_block, _NODATA_YIELD)]) + yield_sum *= pixel_area_ha production_lookup[yield_percentile_id] = yield_sum result_table.write(",%f" % yield_sum) @@ -1080,7 +1084,8 @@ def tabulate_results( (landcover_raster_path, 1)): if landcover_nodata is not None: total_area += numpy.count_nonzero( - ~pygeoprocessing.array_equals_nodata(band_values, landcover_nodata)) + ~pygeoprocessing.array_equals_nodata(band_values, + landcover_nodata)) else: total_area += band_values.size result_table.write( @@ -1090,8 +1095,8 @@ def tabulate_results( def aggregate_to_polygons( base_aggregate_vector_path, target_aggregate_vector_path, - landcover_raster_projection, crop_names, - nutrient_df, yield_percentile_headers, output_dir, file_suffix, + landcover_raster_projection, crop_names, nutrient_df, + yield_percentile_headers, pixel_area_ha, output_dir, file_suffix, target_aggregate_table_path): """Write table with aggregate results of yield and nutrient values. @@ -1108,8 +1113,9 @@ def aggregate_to_polygons( nutrient_df (pandas.DataFrame): a table of nutrient values by crop yield_percentile_headers (list): list of strings indicating percentiles at which yield was calculated. + pixel_area_ha (float): area of lulc raster cells (hectares) output_dir (string): the file path to the output workspace. - file_suffix (string): string to appened to any output filenames. + file_suffix (string): string to append to any output filenames. target_aggregate_table_path (string): path to 'aggregate_results.csv' in the output workspace @@ -1124,6 +1130,10 @@ def aggregate_to_polygons( target_aggregate_vector_path, driver_name='ESRI Shapefile') + # Since pixel values are Mg/(ha•yr), zonal stats sum is (Mg•px)/(ha•yr). + # Before writing sum to results tables or when using sum to calculate + # nutrient yields, convert to Mg/yr by multiplying by ha/px. + # loop over every crop and query with pgp function total_yield_lookup = {} total_nutrient_table = collections.defaultdict( @@ -1153,11 +1163,12 @@ def aggregate_to_polygons( crop_name, yield_percentile_id)]: total_nutrient_table[nutrient_id][ yield_percentile_id][id_index] += ( - nutrient_factor * - total_yield_lookup['%s_%s' % ( - crop_name, yield_percentile_id)][ - id_index]['sum'] * - nutrient_df[nutrient_id][crop_name]) + nutrient_factor + * total_yield_lookup[ + '%s_%s' % (crop_name, yield_percentile_id) + ][id_index]['sum'] + * pixel_area_ha + * nutrient_df[nutrient_id][crop_name]) # process observed observed_yield_path = os.path.join( @@ -1171,10 +1182,11 @@ def aggregate_to_polygons( for id_index in total_yield_lookup[f'{crop_name}_observed']: total_nutrient_table[ nutrient_id]['observed'][id_index] += ( - nutrient_factor * - total_yield_lookup[ - f'{crop_name}_observed'][id_index]['sum'] * - nutrient_df[nutrient_id][crop_name]) + nutrient_factor + * total_yield_lookup[ + f'{crop_name}_observed'][id_index]['sum'] + * pixel_area_ha + * nutrient_df[nutrient_id][crop_name]) # report everything to a table with open(target_aggregate_table_path, 'w') as aggregate_table: @@ -1193,7 +1205,8 @@ def aggregate_to_polygons( for id_index in list(total_yield_lookup.values())[0]: aggregate_table.write('%s,' % id_index) aggregate_table.write(','.join([ - str(total_yield_lookup[yield_header][id_index]['sum']) + str(total_yield_lookup[yield_header][id_index]['sum'] + * pixel_area_ha) for yield_header in sorted(total_yield_lookup)])) for nutrient_id in _EXPECTED_NUTRIENT_TABLE_HEADERS: diff --git a/src/natcap/invest/crop_production_regression.py b/src/natcap/invest/crop_production_regression.py index b795ee0d8b..3c18c221fe 100644 --- a/src/natcap/invest/crop_production_regression.py +++ b/src/natcap/invest/crop_production_regression.py @@ -1,4 +1,4 @@ -"""InVEST Crop Production Percentile Model.""" +"""InVEST Crop Production Regression Model.""" import collections import logging import os @@ -218,8 +218,8 @@ "about": f"{x} {name} production within the polygon", "type": "number", "units": units - } for nutrient, name, units in NUTRIENTS - for x in ["modeled", "observed"] + } for (nutrient, name, units) in NUTRIENTS + for x in ["modeled", "observed"] } } }, @@ -251,8 +251,8 @@ "about": f"{x} {name} production from the crop", "type": "number", "units": units - } for nutrient, name, units in NUTRIENTS - for x in ["modeled", "observed"] + } for (nutrient, name, units) in NUTRIENTS + for x in ["modeled", "observed"] } } }, @@ -660,8 +660,7 @@ def execute(args): (regression_parameter_raster_path_lookup['c_n'], 1), (args['landcover_raster_path'], 1), (crop_to_fertilization_rate_df['nitrogen_rate'][crop_name], - 'raw'), - (crop_lucode, 'raw'), (pixel_area_ha, 'raw')], + 'raw'), (crop_lucode, 'raw')], _x_yield_op, nitrogen_yield_raster_path, gdal.GDT_Float32, _NODATA_YIELD), target_path_list=[nitrogen_yield_raster_path], @@ -679,8 +678,7 @@ def execute(args): (regression_parameter_raster_path_lookup['c_p2o5'], 1), (args['landcover_raster_path'], 1), (crop_to_fertilization_rate_df['phosphorus_rate'][crop_name], - 'raw'), - (crop_lucode, 'raw'), (pixel_area_ha, 'raw')], + 'raw'), (crop_lucode, 'raw')], _x_yield_op, phosphorus_yield_raster_path, gdal.GDT_Float32, _NODATA_YIELD), target_path_list=[phosphorus_yield_raster_path], @@ -698,8 +696,7 @@ def execute(args): (regression_parameter_raster_path_lookup['c_k2o'], 1), (args['landcover_raster_path'], 1), (crop_to_fertilization_rate_df['potassium_rate'][crop_name], - 'raw'), - (crop_lucode, 'raw'), (pixel_area_ha, 'raw')], + 'raw'), (crop_lucode, 'raw')], _x_yield_op, potassium_yield_raster_path, gdal.GDT_Float32, _NODATA_YIELD), target_path_list=[potassium_yield_raster_path], @@ -794,7 +791,7 @@ def execute(args): args=([(args['landcover_raster_path'], 1), (interpolated_observed_yield_raster_path, 1), (observed_yield_nodata, 'raw'), (landcover_nodata, 'raw'), - (crop_lucode, 'raw'), (pixel_area_ha, 'raw')], + (crop_lucode, 'raw')], _mask_observed_yield_op, observed_production_raster_path, gdal.GDT_Float32, observed_yield_nodata), target_path_list=[observed_production_raster_path], @@ -834,10 +831,10 @@ def execute(args): func=aggregate_regression_results_to_polygons, args=(args['aggregate_polygon_path'], target_aggregate_vector_path, + aggregate_results_table_path, landcover_raster_info['projection_wkt'], - crop_names, nutrient_df, - output_dir, file_suffix, - aggregate_results_table_path), + crop_names, nutrient_df, pixel_area_ha, + output_dir, file_suffix), target_path_list=[target_aggregate_vector_path, aggregate_results_table_path], dependent_task_list=dependent_task_list, @@ -848,12 +845,12 @@ def execute(args): def _x_yield_op( - y_max, b_x, c_x, lulc_array, fert_rate, crop_lucode, pixel_area_ha): + y_max, b_x, c_x, lulc_array, fert_rate, crop_lucode): """Calc generalized yield op, Ymax*(1-b_NP*exp(-cN * N_GC)). The regression model has identical mathematical equations for - the nitrogen, phosphorus, and potassium. The only difference is - the scalars in the equation (fertilization rate and pixel area). + the nitrogen, phosphorus, and potassium. The only difference is + the scalar in the equation (fertilization rate). """ result = numpy.empty(b_x.shape, dtype=numpy.float32) result[:] = _NODATA_YIELD @@ -862,7 +859,7 @@ def _x_yield_op( ~pygeoprocessing.array_equals_nodata(b_x, _NODATA_YIELD) & ~pygeoprocessing.array_equals_nodata(c_x, _NODATA_YIELD) & (lulc_array == crop_lucode)) - result[valid_mask] = pixel_area_ha * y_max[valid_mask] * ( + result[valid_mask] = y_max[valid_mask] * ( 1 - b_x[valid_mask] * numpy.exp( -c_x[valid_mask] * fert_rate)) @@ -897,7 +894,7 @@ def _zero_observed_yield_op(observed_yield_array, observed_yield_nodata): def _mask_observed_yield_op( lulc_array, observed_yield_array, observed_yield_nodata, - landcover_nodata, crop_lucode, pixel_area_ha): + landcover_nodata, crop_lucode): """Mask total observed yield to crop lulc type. Args: @@ -906,7 +903,6 @@ def _mask_observed_yield_op( observed_yield_nodata (float): yield raster nodata value landcover_nodata (float): landcover raster nodata value crop_lucode (int): code used to mask in the current crop - pixel_area_ha (float): area of lulc raster cells (hectares) Returns: numpy.ndarray with float values of yields masked to crop_lucode @@ -915,13 +911,13 @@ def _mask_observed_yield_op( result = numpy.empty(lulc_array.shape, dtype=numpy.float32) if landcover_nodata is not None: result[:] = observed_yield_nodata - valid_mask = ~pygeoprocessing.array_equals_nodata(lulc_array, landcover_nodata) + valid_mask = ~pygeoprocessing.array_equals_nodata( + lulc_array, landcover_nodata) result[valid_mask] = 0 else: result[:] = 0 lulc_mask = lulc_array == crop_lucode - result[lulc_mask] = ( - observed_yield_array[lulc_mask] * pixel_area_ha) + result[lulc_mask] = observed_yield_array[lulc_mask] return result @@ -952,6 +948,11 @@ def tabulate_regression_results( nutrient_id + '_' + mode for nutrient_id in _EXPECTED_NUTRIENT_TABLE_HEADERS for mode in ['modeled', 'observed']] + + # Since pixel values in observed and percentile rasters are Mg/(ha•yr), + # raster sums are (Mg•px)/(ha•yr). Before recording sums in + # production_lookup dictionary, convert to Mg/yr by multiplying by ha/px. + with open(target_table_path, 'w') as result_table: result_table.write( 'crop,area (ha),' + 'production_observed,production_modeled,' + @@ -982,6 +983,7 @@ def tabulate_regression_results( production_pixel_count += numpy.count_nonzero( valid_mask & (yield_block > 0.0)) yield_sum += numpy.sum(yield_block[valid_mask]) + yield_sum *= pixel_area_ha production_area = production_pixel_count * pixel_area_ha production_lookup['observed'] = yield_sum result_table.write(',%f' % production_area) @@ -997,6 +999,7 @@ def tabulate_regression_results( # _NODATA_YIELD will always have a value (defined above) yield_block[~pygeoprocessing.array_equals_nodata( yield_block, _NODATA_YIELD)]) + yield_sum *= pixel_area_ha production_lookup['modeled'] = yield_sum result_table.write(",%f" % yield_sum) @@ -1031,9 +1034,8 @@ def tabulate_regression_results( def aggregate_regression_results_to_polygons( base_aggregate_vector_path, target_aggregate_vector_path, - landcover_raster_projection, crop_names, - nutrient_df, output_dir, file_suffix, - target_aggregate_table_path): + aggregate_results_table_path, landcover_raster_projection, + crop_names, nutrient_df, pixel_area_ha, output_dir, file_suffix): """Write table with aggregate results of yield and nutrient values. Use zonal statistics to summarize total observed and interpolated @@ -1042,15 +1044,16 @@ def aggregate_regression_results_to_polygons( Args: base_aggregate_vector_path (string): path to polygon vector - target_aggregate_vector_path (string): - path to re-projected copy of polygon vector + target_aggregate_vector_path (string): path to re-projected copy of + polygon vector + aggregate_results_table_path (string): path to CSV file where aggregate + results will be reported. landcover_raster_projection (string): a WKT projection string crop_names (list): list of crop names nutrient_df (pandas.DataFrame): a table of nutrient values by crop + pixel_area_ha (float): area of lulc raster cells (hectares) output_dir (string): the file path to the output workspace. file_suffix (string): string to append to any output filenames. - target_aggregate_table_path (string): path to 'aggregate_results.csv' - in the output workspace Returns: None @@ -1062,6 +1065,10 @@ def aggregate_regression_results_to_polygons( target_aggregate_vector_path, driver_name='ESRI Shapefile') + # Since pixel values are Mg/(ha•yr), zonal stats sum is (Mg•px)/(ha•yr). + # Before writing sum to results tables or when using sum to calculate + # nutrient yields, convert to Mg/yr by multiplying by ha/px. + # loop over every crop and query with pgp function total_yield_lookup = {} total_nutrient_table = collections.defaultdict( @@ -1085,10 +1092,11 @@ def aggregate_regression_results_to_polygons( for fid_index in total_yield_lookup['%s_modeled' % crop_name]: total_nutrient_table[nutrient_id][ 'modeled'][fid_index] += ( - nutrient_factor * - total_yield_lookup['%s_modeled' % crop_name][ - fid_index]['sum'] * - nutrient_df[nutrient_id][crop_name]) + nutrient_factor + * total_yield_lookup['%s_modeled' % crop_name][ + fid_index]['sum'] + * pixel_area_ha + * nutrient_df[nutrient_id][crop_name]) # process observed observed_yield_path = os.path.join( @@ -1103,15 +1111,14 @@ def aggregate_regression_results_to_polygons( '%s_observed' % crop_name]: total_nutrient_table[ nutrient_id]['observed'][fid_index] += ( - nutrient_factor * # percent crop used * 1000 [100g per Mg] - total_yield_lookup[ - '%s_observed' % crop_name][fid_index]['sum'] * - nutrient_df[nutrient_id][crop_name]) # nutrient unit per 100g crop + nutrient_factor # percent crop used * 1000 [100g per Mg] + * total_yield_lookup[ + '%s_observed' % crop_name][fid_index]['sum'] + * pixel_area_ha + * nutrient_df[nutrient_id][crop_name]) # nutrient unit per 100g crop # report everything to a table - aggregate_table_path = os.path.join( - output_dir, _AGGREGATE_TABLE_FILE_PATTERN % file_suffix) - with open(aggregate_table_path, 'w') as aggregate_table: + with open(aggregate_results_table_path, 'w') as aggregate_table: # write header aggregate_table.write('FID,') aggregate_table.write(','.join(sorted(total_yield_lookup)) + ',') @@ -1127,7 +1134,8 @@ def aggregate_regression_results_to_polygons( for id_index in list(total_yield_lookup.values())[0]: aggregate_table.write('%s,' % id_index) aggregate_table.write(','.join([ - str(total_yield_lookup[yield_header][id_index]['sum']) + str(total_yield_lookup[yield_header][id_index]['sum'] + * pixel_area_ha) for yield_header in sorted(total_yield_lookup)])) for nutrient_id in _EXPECTED_NUTRIENT_TABLE_HEADERS: diff --git a/src/natcap/invest/forest_carbon_edge_effect.py b/src/natcap/invest/forest_carbon_edge_effect.py index 9d18eb26af..d506d83420 100644 --- a/src/natcap/invest/forest_carbon_edge_effect.py +++ b/src/natcap/invest/forest_carbon_edge_effect.py @@ -7,7 +7,6 @@ import os import pickle import time -import uuid import numpy import pandas @@ -187,13 +186,13 @@ "outputs": { "carbon_map.tif": { "about": ( - "A map of carbon stock per pixel, with the amount in forest derived from the regression based on " - "distance to forest edge, and the amount in non-forest classes according to the biophysical table. " - "Note that because the map displays carbon per pixel, coarser resolution maps should have higher " - "values for carbon, because the pixel areas are larger."), + "A map of carbon stock per hectare, with the amount in forest " + "derived from the regression based on distance to forest " + "edge, and the amount in non-forest classes according to the " + "biophysical table. "), "bands": {1: { "type": "number", - "units": u.metric_ton/u.pixel + "units": u.metric_ton/u.hectare }} }, "aggregated_carbon_stocks.shp": { @@ -205,7 +204,7 @@ "units": u.metric_ton, "about": "Total carbon in the area." }, - "c_ha_mean":{ + "c_ha_mean": { "type": "number", "units": u.metric_ton/u.hectare, "about": "Mean carbon density in the area." @@ -580,17 +579,24 @@ def _aggregate_carbon_map( target_aggregate_layer.ResetReading() target_aggregate_layer.StartTransaction() + # Since pixel values are Mg/ha, raster sum is (Mg•px)/ha. + # To convert to Mg, multiply by ha/px. + raster_info = pygeoprocessing.get_raster_info(carbon_map_path) + pixel_area = abs(numpy.prod(raster_info['pixel_size'])) + ha_per_px = pixel_area / 10000 + for poly_feat in target_aggregate_layer: poly_fid = poly_feat.GetFID() poly_feat.SetField( - 'c_sum', float(serviceshed_stats[poly_fid]['sum'])) + 'c_sum', float(serviceshed_stats[poly_fid]['sum'] * ha_per_px)) # calculates mean pixel value per ha in for each feature in AOI poly_geom = poly_feat.GetGeometryRef() poly_area_ha = poly_geom.GetArea() / 1e4 # converts m^2 to hectare poly_geom = None poly_feat.SetField( 'c_ha_mean', - float(serviceshed_stats[poly_fid]['sum'] / poly_area_ha)) + float(serviceshed_stats[poly_fid]['sum'] / poly_area_ha + * ha_per_px)) target_aggregate_layer.SetFeature(poly_feat) target_aggregate_layer.CommitTransaction() @@ -629,9 +635,6 @@ def _calculate_lulc_carbon_map( biophysical_table_path, **MODEL_SPEC['args']['biophysical_table_path']) lucode_to_per_cell_carbon = {} - cell_size = pygeoprocessing.get_raster_info( - lulc_raster_path)['pixel_size'] # in meters - cell_area_ha = abs(cell_size[0]) * abs(cell_size[1]) / 10000 # Build a lookup table for lucode, row in biophysical_df.iterrows(): @@ -648,8 +651,7 @@ def _calculate_lulc_carbon_map( "Could not interpret carbon pool value as a number. " f"lucode: {lucode}, pool_type: {carbon_pool_type}, " f"value: {row[carbon_pool_type]}") - lucode_to_per_cell_carbon[lucode] = row[carbon_pool_type] * cell_area_ha - + lucode_to_per_cell_carbon[lucode] = row[carbon_pool_type] # map aboveground carbon from table to lulc that is not forest reclass_error_details = { @@ -873,7 +875,6 @@ def _calculate_tropical_forest_edge_carbon_map( cell_xsize, cell_ysize = pygeoprocessing.get_raster_info( edge_distance_path)['pixel_size'] cell_size_km = (abs(cell_xsize) + abs(cell_ysize))/2 / 1000 - cell_area_ha = (abs(cell_xsize) * abs(cell_ysize)) / 10000 # Loop memory block by memory block, calculating the forest edge carbon # for every forest pixel. @@ -957,19 +958,19 @@ def _calculate_tropical_forest_edge_carbon_map( biomass[mask_1] = ( thetas[mask_1][:, 0] - thetas[mask_1][:, 1] * numpy.exp( -thetas[mask_1][:, 2] * valid_edge_distances_km[mask_1]) - ) * cell_area_ha + ) # logarithmic model # biomass_2 = t1 + t2 * numpy.log(edge_dist_km) biomass[mask_2] = ( thetas[mask_2][:, 0] + thetas[mask_2][:, 1] * numpy.log( - valid_edge_distances_km[mask_2])) * cell_area_ha + valid_edge_distances_km[mask_2])) # linear regression # biomass_3 = t1 + t2 * edge_dist_km biomass[mask_3] = ( thetas[mask_3][:, 0] + thetas[mask_3][:, 1] * - valid_edge_distances_km[mask_3]) * cell_area_ha + valid_edge_distances_km[mask_3]) # reshape the array so that each set of points is in a separate # dimension, here distances are distances to each valid model diff --git a/src/natcap/invest/ndr/ndr.py b/src/natcap/invest/ndr/ndr.py index 4c39f31e89..b552c4ad14 100644 --- a/src/natcap/invest/ndr/ndr.py +++ b/src/natcap/invest/ndr/ndr.py @@ -1,6 +1,5 @@ """InVEST Nutrient Delivery Ratio (NDR) module.""" import copy -import itertools import logging import os import pickle @@ -94,16 +93,13 @@ "about": gettext( "The distance after which it is assumed that this " "LULC type retains the nutrient at its maximum " - "capacity. If nutrients travel a shorter distance " - "that this, the retention " - "efficiency will be less than the maximum value " - "eff_x, following an exponential decay.")}, + "capacity.")}, "proportion_subsurface_n": { "type": "ratio", "required": "calc_n", "about": gettext( "The proportion of the total amount of nitrogen that " - "are dissolved into the subsurface. By default, this " + "is dissolved into the subsurface. By default, this " "value should be set to 0, indicating that all " "nutrients are delivered via surface flow. There is " "no equivalent of this for phosphorus.")} @@ -210,28 +206,28 @@ "about": "A pixel level map showing how much phosphorus from each pixel eventually reaches the stream by surface flow.", "bands": {1: { "type": "number", - "units": u.kilogram/u.pixel + "units": u.kilogram/u.hectare }} }, "n_surface_export.tif": { "about": "A pixel level map showing how much nitrogen from each pixel eventually reaches the stream by surface flow.", "bands": {1: { "type": "number", - "units": u.kilogram/u.pixel + "units": u.kilogram/u.hectare }} }, "n_subsurface_export.tif": { "about": "A pixel level map showing how much nitrogen from each pixel eventually reaches the stream by subsurface flow.", "bands": {1: { "type": "number", - "units": u.kilogram/u.pixel + "units": u.kilogram/u.hectare }} }, "n_total_export.tif": { "about": "A pixel level map showing how much nitrogen from each pixel eventually reaches the stream by either flow.", "bands": {1: { "type": "number", - "units": u.kilogram/u.pixel + "units": u.kilogram/u.hectare }} }, "intermediate_outputs": { @@ -1032,12 +1028,15 @@ def execute(args): # raster_map equation: Multiply a series of arrays element-wise def _mult_op(*array_list): return numpy.prod(numpy.stack(array_list), axis=0) + # raster_map equation: Sum a list of arrays element-wise def _sum_op(*array_list): return numpy.sum(array_list, axis=0) + # raster_map equation: calculate inverse of S factor def _inverse_op(base_val): return numpy.where(base_val == 0, 0, 1 / base_val) + # raster_map equation: rescale and threshold slope between 0.005 and 1 def _slope_proportion_and_threshold_op(slope): slope_fraction = slope / 100 @@ -1138,8 +1137,13 @@ def _add_fields_to_shapefile(field_pickle_map, target_vector_path): for feature in target_layer: fid = feature.GetFID() for field_name in field_pickle_map: + # Since pixel values are kg/(ha•yr), raster sum is (kg•px)/(ha•yr). + # To convert to kg/yr, multiply by ha/px. + pixel_area = field_summaries[field_name]['pixel_area'] + ha_per_px = pixel_area / 10000 feature.SetField( - field_name, float(field_summaries[field_name][fid]['sum'])) + field_name, float( + field_summaries[field_name][fid]['sum']) * ha_per_px) # Save back to datasource target_layer.SetFeature(feature) target_layer = None @@ -1207,8 +1211,8 @@ def _normalize_raster(base_raster_path_band, target_normalized_raster_path): """ value_sum, value_count = pygeoprocessing.raster_reduce( - function=lambda sum_count, block: # calculate both in one pass - (sum_count[0] + numpy.sum(block), sum_count[1] + block.size), + function=lambda sum_count, block: ( # calculate both in one pass + sum_count[0] + numpy.sum(block), sum_count[1] + block.size), raster_path_band=base_raster_path_band, initializer=(0, 0)) @@ -1224,29 +1228,25 @@ def _normalize_raster(base_raster_path_band, target_normalized_raster_path): def _calculate_load(lulc_raster_path, lucode_to_load, target_load_raster): - """Calculate load raster by mapping landcover and multiplying by area. + """Calculate load raster by mapping landcover. Args: lulc_raster_path (string): path to integer landcover raster. lucode_to_load (dict): a mapping of landcover IDs to per-area nutrient load. target_load_raster (string): path to target raster that will have - total load per pixel. + load values (kg/ha) mapped to pixels based on LULC. Returns: None. """ - cell_area_ha = abs(numpy.prod(pygeoprocessing.get_raster_info( - lulc_raster_path)['pixel_size'])) * 0.0001 - def _map_load_op(lucode_array): """Convert unit load to total load & handle nodata.""" result = numpy.empty(lucode_array.shape) for lucode in numpy.unique(lucode_array): try: - result[lucode_array == lucode] = ( - lucode_to_load[lucode] * cell_area_ha) + result[lucode_array == lucode] = (lucode_to_load[lucode]) except KeyError: raise KeyError( 'lucode: %d is present in the landuse raster but ' @@ -1442,6 +1442,12 @@ def _aggregate_and_pickle_total( base_raster_path_band, aggregate_vector_path, working_dir=os.path.dirname(target_pickle_path)) + # Write pixel area to pickle file so that _add_fields_to_shapefile + # can adjust totals as needed. + raster_info = pygeoprocessing.get_raster_info(base_raster_path_band[0]) + pixel_area = abs(numpy.prod(raster_info['pixel_size'])) + result['pixel_area'] = pixel_area + with open(target_pickle_path, 'wb') as target_pickle_file: pickle.dump(result, target_pickle_file) diff --git a/src/natcap/invest/sdr/sdr.py b/src/natcap/invest/sdr/sdr.py index f1540982b7..a6fb96c01e 100644 --- a/src/natcap/invest/sdr/sdr.py +++ b/src/natcap/invest/sdr/sdr.py @@ -147,20 +147,20 @@ "about": "The contribution of vegetation to keeping soil from eroding from each pixel. (Eq. (82))", "bands": {1: { "type": "number", - "units": u.metric_ton/u.pixel + "units": u.metric_ton/u.hectare }} }, "avoided_export.tif": { "about": "The contribution of vegetation to keeping erosion from entering a stream. This combines local/on-pixel sediment retention with trapping of erosion from upslope of the pixel. (Eq. (83))", "bands": {1: { "type": "number", - "units": u.metric_ton/u.pixel + "units": u.metric_ton/u.hectare }} }, "rkls.tif": { "bands": {1: { "type": "number", - "units": u.metric_ton/u.pixel + "units": u.metric_ton/u.hectare }}, "about": "Total potential soil loss per pixel in the original land cover from the RKLS equation. Equivalent to the soil loss for bare soil. (Eq. (68), without applying the C or P factors)." }, @@ -168,14 +168,14 @@ "about": "The total amount of sediment deposited on the pixel from upslope sources as a result of trapping. (Eq. (80))", "bands": {1: { "type": "number", - "units": u.metric_ton/u.pixel + "units": u.metric_ton/u.hectare }} }, "sed_export.tif": { "about": "The total amount of sediment exported from each pixel that reaches the stream. (Eq. (76))", "bands": {1: { "type": "number", - "units": u.metric_ton/u.pixel + "units": u.metric_ton/u.hectare }} }, "stream.tif": spec_utils.STREAM, @@ -185,10 +185,10 @@ "bands": {1: {"type": "integer"}} }, "usle.tif": { - "about": "Total potential soil loss per pixel in the original land cover calculated from the USLE equation. (Eq. (68))", + "about": "Total potential soil loss per hectare in the original land cover calculated from the USLE equation. (Eq. (68))", "bands": {1: { "type": "number", - "units": u.metric_ton/u.pixel + "units": u.metric_ton/u.hectare }} }, "watershed_results_sdr.shp": { @@ -960,6 +960,7 @@ def _avoided_export_op(avoided_erosion, sdr, sed_deposition): # known as a modified USLE) return avoided_erosion * sdr + sed_deposition + def add_drainage_op(stream, drainage): """raster_map equation: add drainage mask to stream layer. @@ -974,9 +975,11 @@ def add_drainage_op(stream, drainage): """ return numpy.where(drainage == 1, 1, stream) + # raster_map equation: calculate USLE def usle_op(rkls, cp_factor): return rkls * cp_factor + # raster_map equation: calculate the inverse ws factor def inverse_ws_op(w_factor, s_factor): return 1 / (w_factor * s_factor) @@ -1185,7 +1188,7 @@ def ls_factor_function(percent_slope, flow_accumulation): def _calculate_rkls( ls_factor_path, erosivity_path, erodibility_path, stream_path, rkls_path): - """Calculate potential soil loss (tons / (pixel * year)) using RKLS. + """Calculate potential soil loss (tons / (ha * year)) using RKLS. (revised universal soil loss equation with no C or P). @@ -1211,10 +1214,6 @@ def _calculate_rkls( stream_nodata = pygeoprocessing.get_raster_info( stream_path)['nodata'][0] - cell_size = abs( - pygeoprocessing.get_raster_info(ls_factor_path)['pixel_size'][0]) - cell_area_ha = cell_size**2 / 10000.0 # hectares per pixel - def rkls_function(ls_factor, erosivity, erodibility, stream): """Calculate the RKLS equation. @@ -1227,7 +1226,7 @@ def rkls_function(ls_factor, erosivity, erodibility, stream): stream (numpy.ndarray): stream mask (1 stream, 0 no stream) Returns: - numpy.ndarray of RKLS values in tons / (pixel * year)) + numpy.ndarray of RKLS values in tons / (ha * year)) """ rkls = numpy.empty(ls_factor.shape, dtype=numpy.float32) nodata_mask = ( @@ -1243,11 +1242,10 @@ def rkls_function(ls_factor, erosivity, erodibility, stream): valid_mask = nodata_mask & (stream == 0) rkls[:] = _TARGET_NODATA - rkls[valid_mask] = ( # rkls units are tons / (pixel * year) + rkls[valid_mask] = ( # rkls units are tons / (ha * 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 + erodibility[valid_mask]) # t * ha * hr / (MJ * ha * mm) return rkls # aligning with index 3 that's the stream and the most likely to be @@ -1386,18 +1384,6 @@ def _calculate_d_up( target_path=out_d_up_path) -def _calculate_d_up_bare( - s_bar_path, flow_accumulation_path, out_d_up_bare_path): - """Calculate s_bar * sqrt(flow accumulation * cell area).""" - cell_area = abs( - pygeoprocessing.get_raster_info(s_bar_path)['pixel_size'][0])**2 - pygeoprocessing.raster_map( - op=lambda s_bar, flow_accum: ( - numpy.sqrt(flow_accum * cell_area) * s_bar), - rasters=[s_bar_path, flow_accumulation_path], - target_path=out_d_up_bare_path) - - def _calculate_ic(d_up_path, d_dn_path, out_ic_factor_path): """Calculate log10(d_up/d_dn).""" # ic can be positive or negative, so float.min is a reasonable nodata value @@ -1521,13 +1507,20 @@ def _generate_report( field_def.SetPrecision(11) target_layer.CreateField(field_def) + # Since pixel values are t/(ha•yr), raster sum is (t•px)/(ha•yr). + # To convert to t/yr, multiply by ha/px. + raster_info = pygeoprocessing.get_raster_info(usle_path) + pixel_area = abs(numpy.prod(raster_info['pixel_size'])) + ha_per_px = pixel_area / 10000 + target_layer.ResetReading() for feature in target_layer: feature_id = feature.GetFID() for field_name in field_summaries: feature.SetField( field_name, - float(field_summaries[field_name][feature_id]['sum'])) + float(field_summaries[field_name][feature_id]['sum'] + * ha_per_px)) target_layer.SetFeature(feature) target_vector = None target_layer = None diff --git a/tests/test_carbon.py b/tests/test_carbon.py index ed0e111094..bee15ff1c8 100644 --- a/tests/test_carbon.py +++ b/tests/test_carbon.py @@ -4,16 +4,19 @@ import tempfile import shutil import os +import re from osgeo import gdal from osgeo import osr import numpy import numpy.random import numpy.testing + import pygeoprocessing gdal.UseExceptions() + def make_simple_raster(base_raster_path, fill_val, nodata_val): """Create a 10x10 raster on designated path with fill value. @@ -68,7 +71,8 @@ def assert_raster_equal_value(base_raster_path, val_to_compare): array_to_compare = numpy.empty(base_array.shape) array_to_compare.fill(val_to_compare) - numpy.testing.assert_allclose(base_array, array_to_compare, rtol=0, atol=1e-6) + numpy.testing.assert_allclose(base_array, array_to_compare, + rtol=0, atol=1e-3) def make_pools_csv(pools_csv_path): @@ -90,6 +94,28 @@ def make_pools_csv(pools_csv_path): open_table.write('2,1,5,0,3,"lulc code 3"\n') +def assert_aggregate_result_equal(html_report_path, stat_name, val_to_compare): + """Assert that the given stat in the HTML report has a specific value. + + Args: + html_report_path (str): path to the HTML report generated by the model. + stat_name (str): name of the stat to find. Must match the name listed + in the HTML. + val_to_compare (float): the value to check against. + + Returns: + None. + """ + with open(html_report_path) as file: + report = file.read() + pattern = (r'data-summary-stat="' + + stat_name + + r'">(\-?\d+\.\d{2})') + match = re.search(pattern, report) + stat_str = match.groups()[0] + assert float(stat_str) == val_to_compare + + class CarbonTests(unittest.TestCase): """Tests for the Carbon Model.""" @@ -131,12 +157,46 @@ def test_carbon_full(self): carbon.execute(args) - # Add assertions for npv for future and REDD scenarios. - # The npv was calculated based on _calculate_npv in carbon.py. + # Ensure every pixel has the correct total C value. + # Current: 15 + 10 + 60 + 1 = 86 Mg/ha + assert_raster_equal_value( + os.path.join(args['workspace_dir'], 'tot_c_cur.tif'), 86) + # Future: 5 + 3 + 20 + 0 = 28 Mg/ha + assert_raster_equal_value( + os.path.join(args['workspace_dir'], 'tot_c_fut.tif'), 28) + # REDD: 2 + 1 + 5 + 0 = 8 Mg/ha + assert_raster_equal_value( + os.path.join(args['workspace_dir'], 'tot_c_redd.tif'), 8) + + # Ensure deltas are correct. assert_raster_equal_value( - os.path.join(args['workspace_dir'], 'npv_fut.tif'), -0.3422078) + os.path.join(args['workspace_dir'], 'delta_cur_fut.tif'), -58) assert_raster_equal_value( - os.path.join(args['workspace_dir'], 'npv_redd.tif'), -0.4602106) + os.path.join(args['workspace_dir'], 'delta_cur_redd.tif'), -78) + + # Ensure NPV calculations are correct. + # Valuation constant based on provided args is 59.00136. + # Future: 59.00136 * -58 Mg/ha = -3422.079 Mg/ha + assert_raster_equal_value( + os.path.join(args['workspace_dir'], 'npv_fut.tif'), -3422.079) + # REDD: 59.00136 * -78 Mg/ha = -4602.106 Mg/ha + assert_raster_equal_value( + os.path.join(args['workspace_dir'], 'npv_redd.tif'), -4602.106) + + # Ensure aggregate results are correct. + report_path = os.path.join(args['workspace_dir'], 'report.html') + # Raster size is 100 m^2; therefore, raster total is as follows: + # (x Mg / 1 ha) * (1 ha / 10000 m^2) * (100 m^2) = (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_zero_rates(self): """Carbon: test with 0 discount and rate change.""" @@ -168,16 +228,16 @@ def test_carbon_zero_rates(self): # Add assertions for npv for future and REDD scenarios. # carbon change from cur to fut: - # -58 Mg/ha * .0001 ha/pixel * 43 $/Mg = -0.2494 $/pixel + # -58 Mg/ha * 43 $/Mg = -2494 $/ha assert_raster_equal_value( - os.path.join(args['workspace_dir'], 'npv_fut.tif'), -0.2494) + os.path.join(args['workspace_dir'], 'npv_fut.tif'), -2494) # carbon change from cur to redd: - # -78 Mg/ha * .0001 ha/pixel * 43 $/Mg = -0.3354 $/pixel + # -78 Mg/ha * 43 $/Mg = -3354 $/ha assert_raster_equal_value( - os.path.join(args['workspace_dir'], 'npv_redd.tif'), -0.3354) + os.path.join(args['workspace_dir'], 'npv_redd.tif'), -3354) - def test_carbon_future(self): - """Carbon: regression testing future scenario.""" + def test_carbon_without_redd(self): + """Carbon: regression testing for future scenario only (no REDD).""" from natcap.invest import carbon args = { 'workspace_dir': self.workspace_dir, @@ -201,10 +261,12 @@ def test_carbon_future(self): make_pools_csv(args['carbon_pools_path']) carbon.execute(args) - # Add assertions for npv for the future scenario. - # The npv was calculated based on _calculate_npv in carbon.py. + + # Ensure NPV calculations are correct. + # Valuation constant based on provided args is 59.00136. + # Future: 59.00136 * -58 Mg/ha = -3422.079 Mg/ha assert_raster_equal_value( - os.path.join(args['workspace_dir'], 'npv_fut.tif'), -0.3422078) + os.path.join(args['workspace_dir'], 'npv_fut.tif'), -3422.079) def test_carbon_missing_landcover_values(self): """Carbon: testing expected exception on missing LULC codes.""" @@ -261,12 +323,14 @@ def test_carbon_full_undefined_nodata(self): carbon.execute(args) - # Add assertions for npv for future and REDD scenarios. - # The npv was calculated based on _calculate_npv in carbon.py. + # Ensure NPV calculations are correct. + # Valuation constant based on provided args is 59.00136. + # Future: 59.00136 * -58 Mg/ha = -3422.079 Mg/ha assert_raster_equal_value( - os.path.join(args['workspace_dir'], 'npv_fut.tif'), -0.3422078) + os.path.join(args['workspace_dir'], 'npv_fut.tif'), -3422.079) + # REDD: 59.00136 * -78 Mg/ha = -4602.106 Mg/ha assert_raster_equal_value( - os.path.join(args['workspace_dir'], 'npv_redd.tif'), -0.4602106) + os.path.join(args['workspace_dir'], 'npv_redd.tif'), -4602.106) def test_generate_carbon_map(self): """Test `_generate_carbon_map`""" @@ -313,7 +377,7 @@ def _make_simple_lulc_raster(base_raster_path): band = actual_carbon_stock.GetRasterBand(1) actual_carbon_stock = band.ReadAsArray() - expected_carbon_stock = numpy.array([[0.5, 0.5], [0.006, 0.012]], + expected_carbon_stock = numpy.array([[5000, 5000], [60, 120]], dtype=numpy.float32) numpy.testing.assert_array_equal(actual_carbon_stock, diff --git a/tests/test_crop_production.py b/tests/test_crop_production.py index dc9b340c80..3372b70f1b 100644 --- a/tests/test_crop_production.py +++ b/tests/test_crop_production.py @@ -22,6 +22,20 @@ 'crop_production_model') +def _get_pixels_per_hectare(raster_path): + """Calculate number of pixels per hectare for a given raster. + + Args: + raster_path (str): full path to the raster. + + Returns: + A float representing the number of pixels per hectare. + """ + raster_info = pygeoprocessing.get_raster_info(raster_path) + pixel_area = abs(numpy.prod(raster_info['pixel_size'])) + return 10000 / pixel_area + + def make_aggregate_vector(path_to_shp): """ Generate shapefile with two overlapping polygons @@ -176,6 +190,41 @@ def test_crop_production_percentile(self): pandas.testing.assert_frame_equal( expected_result_table, result_table, check_dtype=False) + # Check raster outputs to make sure values are in Mg/ha. + # Raster sum is (Mg•px)/(ha•yr). + # Result table reports totals in Mg/yr. + # To convert from Mg/yr to (Mg•px)/(ha•yr), multiply by px/ha. + expected_raster_sums = {} + for (index, crop) in [(0, 'barley'), (1, 'soybean'), (2, 'wheat')]: + filename = crop + '_observed_production.tif' + pixels_per_hectare = _get_pixels_per_hectare( + os.path.join(args['workspace_dir'], filename)) + expected_raster_sums[filename] = ( + expected_result_table.loc[index]['production_observed'] + * pixels_per_hectare) + for percentile in ['25', '50', '75', '95']: + filename = ( + crop + '_yield_' + percentile + 'th_production.tif') + col_name = 'production_' + percentile + 'th' + pixels_per_hectare = _get_pixels_per_hectare( + os.path.join(args['workspace_dir'], filename)) + expected_raster_sums[filename] = ( + expected_result_table.loc[index][col_name] + * pixels_per_hectare) + + for filename in expected_raster_sums: + raster_path = os.path.join(args['workspace_dir'], filename) + raster_info = pygeoprocessing.get_raster_info(raster_path) + nodata = raster_info['nodata'][0] + raster_sum = 0.0 + for _, block in pygeoprocessing.iterblocks((raster_path, 1)): + raster_sum += numpy.sum( + block[~pygeoprocessing.array_equals_nodata( + block, nodata)], dtype=numpy.float32) + expected_sum = expected_raster_sums[filename] + numpy.testing.assert_allclose(raster_sum, expected_sum, + rtol=0, atol=0.1) + def test_crop_production_percentile_no_nodata(self): """Crop Production: test percentile model with undefined nodata raster. @@ -324,9 +373,6 @@ def test_crop_production_regression_bad_crop(self): 'model_data_path': MODEL_DATA_PATH, 'fertilization_rate_table_path': os.path.join( SAMPLE_DATA_PATH, 'crop_fertilization_rates.csv'), - 'nitrogen_fertilization_rate': 29.6, - 'phosphorus_fertilization_rate': 8.4, - 'potassium_fertilization_rate': 14.2, 'n_workers': '-1' } @@ -355,9 +401,6 @@ def test_crop_production_regression_missing_climate_bin(self): 'model_data_path': MODEL_DATA_PATH, 'fertilization_rate_table_path': os.path.join( SAMPLE_DATA_PATH, 'crop_fertilization_rates.csv'), - 'nitrogen_fertilization_rate': 29.6, - 'phosphorus_fertilization_rate': 8.4, - 'potassium_fertilization_rate': 14.2, 'n_workers': '-1' } @@ -409,15 +452,13 @@ def test_crop_production_regression(self): 'model_data_path': MODEL_DATA_PATH, 'fertilization_rate_table_path': os.path.join( SAMPLE_DATA_PATH, 'crop_fertilization_rates.csv'), - 'nitrogen_fertilization_rate': 29.6, - 'phosphorus_fertilization_rate': 8.4, - 'potassium_fertilization_rate': 14.2, } crop_production_regression.execute(args) expected_agg_result_table = pandas.read_csv( - os.path.join(TEST_DATA_PATH, 'expected_regression_aggregate_results.csv')) + os.path.join(TEST_DATA_PATH, + 'expected_regression_aggregate_results.csv')) agg_result_table = pandas.read_csv( os.path.join(args['workspace_dir'], 'aggregate_results.csv')) pandas.testing.assert_frame_equal( @@ -435,6 +476,38 @@ def test_crop_production_regression(self): pandas.testing.assert_frame_equal( expected_result_table, result_table, check_dtype=False) + # Check raster outputs to make sure values are in Mg/ha. + # Raster sum is (Mg•px)/(ha•yr). + # Result table reports totals in Mg/yr. + # To convert from Mg/yr to (Mg•px)/(ha•yr), multiply by px/ha. + expected_raster_sums = {} + for (index, crop) in [(0, 'barley'), (1, 'soybean'), (2, 'wheat')]: + filename = crop + '_observed_production.tif' + pixels_per_hectare = _get_pixels_per_hectare( + os.path.join(args['workspace_dir'], filename)) + expected_raster_sums[filename] = ( + expected_result_table.loc[index]['production_observed'] + * pixels_per_hectare) + filename = crop + '_regression_production.tif' + pixels_per_hectare = _get_pixels_per_hectare( + os.path.join(args['workspace_dir'], filename)) + expected_raster_sums[filename] = ( + expected_result_table.loc[index]['production_modeled'] + * pixels_per_hectare) + + for filename in expected_raster_sums: + raster_path = os.path.join(args['workspace_dir'], filename) + raster_info = pygeoprocessing.get_raster_info(raster_path) + nodata = raster_info['nodata'][0] + raster_sum = 0.0 + for _, block in pygeoprocessing.iterblocks((raster_path, 1)): + raster_sum += numpy.sum( + block[~pygeoprocessing.array_equals_nodata( + block, nodata)], dtype=numpy.float32) + expected_sum = expected_raster_sums[filename] + numpy.testing.assert_allclose(raster_sum, expected_sum, + rtol=0, atol=0.001) + def test_crop_production_regression_no_nodata(self): """Crop Production: test regression model with undefined nodata raster. @@ -453,9 +526,6 @@ def test_crop_production_regression_no_nodata(self): 'model_data_path': MODEL_DATA_PATH, 'fertilization_rate_table_path': os.path.join( SAMPLE_DATA_PATH, 'crop_fertilization_rates.csv'), - 'nitrogen_fertilization_rate': 29.6, - 'phosphorus_fertilization_rate': 8.4, - 'potassium_fertilization_rate': 14.2, } # Create a raster based on the test data geotransform, but smaller and @@ -504,12 +574,11 @@ def test_x_yield_op(self): lulc_array = numpy.array([[3, 3, 2], [3, -1, 3]]) fert_rate = 0.6 crop_lucode = 3 - pixel_area_ha = 10 actual_result = _x_yield_op(y_max, b_x, c_x, lulc_array, fert_rate, - crop_lucode, pixel_area_ha) - expected_result = numpy.array([[-1, -19.393047, -1], - [26.776089, -1, 15.1231]]) + crop_lucode) + expected_result = numpy.array([[-1, -1.9393047, -1], + [2.6776089, -1, 1.51231]]) numpy.testing.assert_allclose(actual_result, expected_result) @@ -544,13 +613,12 @@ def test_mask_observed_yield_op(self): landcover_nodata = -9999 crop_lucode = 3 - pixel_area_ha = 10 actual_result = _mask_observed_yield_op( lulc_array, observed_yield_array, observed_yield_nodata, - landcover_nodata, crop_lucode, pixel_area_ha) + landcover_nodata, crop_lucode) - expected_result = numpy.array([[-10, 0, -1], [80, -99990, 0]]) + expected_result = numpy.array([[-1, 0, -1], [8, -9999, 0]]) numpy.testing.assert_allclose(actual_result, expected_result) @@ -563,75 +631,75 @@ def _create_expected_results(): """Creates the expected results DataFrame.""" return pandas.DataFrame([ {'crop': 'corn', 'area (ha)': 20.0, - 'production_observed': 8.0, 'production_modeled': 4.0, - 'protein_modeled': 1562400.0, 'protein_observed': 3124800.0, - 'lipid_modeled': 297600.0, 'lipid_observed': 595200.0, - 'energy_modeled': 17707200.0, 'energy_observed': 35414400.0, - 'ca_modeled': 1004400.0, 'ca_observed': 2008800.0, - 'fe_modeled': 584040.0, 'fe_observed': 1168080.0, - 'mg_modeled': 10416000.0, 'mg_observed': 20832000.0, - 'ph_modeled': 26188800.0, 'ph_observed': 52377600.0, - 'k_modeled': 64244400.0, 'k_observed': 128488800.0, - 'na_modeled': 74400.0, 'na_observed': 148800.0, - 'zn_modeled': 182280.0, 'zn_observed': 364560.0, - 'cu_modeled': 70680.0, 'cu_observed': 141360.0, - 'fl_modeled': 297600.0, 'fl_observed': 595200.0, - 'mn_modeled': 107880.0, 'mn_observed': 215760.0, - 'se_modeled': 3720.0, 'se_observed': 7440.0, - 'vita_modeled': 111600.0, 'vita_observed': 223200.0, - 'betac_modeled': 595200.0, 'betac_observed': 1190400.0, - 'alphac_modeled': 85560.0, 'alphac_observed': 171120.0, - 'vite_modeled': 29760.0, 'vite_observed': 59520.0, - 'crypto_modeled': 59520.0, 'crypto_observed': 119040.0, - 'lycopene_modeled': 13392.0, 'lycopene_observed': 26784.0, - 'lutein_modeled': 2343600.0, 'lutein_observed': 4687200.0, - 'betat_modeled': 18600.0, 'betat_observed': 37200.0, - 'gammat_modeled': 78120.0, 'gammat_observed': 156240.0, - 'deltat_modeled': 70680.0, 'deltat_observed': 141360.0, - 'vitc_modeled': 252960.0, 'vitc_observed': 505920.0, - 'thiamin_modeled': 14880.0, 'thiamin_observed': 29760.0, - 'riboflavin_modeled': 66960.0, 'riboflavin_observed': 133920.0, - 'niacin_modeled': 305040.0, 'niacin_observed': 610080.0, - 'pantothenic_modeled': 33480.0, 'pantothenic_observed': 66960.0, - 'vitb6_modeled': 52080.0, 'vitb6_observed': 104160.0, - 'folate_modeled': 14322000.0, 'folate_observed': 28644000.0, - 'vitb12_modeled': 74400.0, 'vitb12_observed': 148800.0, - 'vitk_modeled': 1525200.0, 'vitk_observed': 3050400.0}, + 'production_observed': 80.0, 'production_modeled': 40.0, + 'protein_modeled': 15624000.0, 'protein_observed': 31248000.0, + 'lipid_modeled': 2976000.0, 'lipid_observed': 5952000.0, + 'energy_modeled': 177072000.0, 'energy_observed': 354144000.0, + 'ca_modeled': 10044000.0, 'ca_observed': 20088000.0, + 'fe_modeled': 5840400.0, 'fe_observed': 11680800.0, + 'mg_modeled': 104160000.0, 'mg_observed': 208320000.0, + 'ph_modeled': 261888000.0, 'ph_observed': 523776000.0, + 'k_modeled': 642444000.0, 'k_observed': 1284888000.0, + 'na_modeled': 744000.0, 'na_observed': 1488000.0, + 'zn_modeled': 1822800.0, 'zn_observed': 3645600.0, + 'cu_modeled': 706800.0, 'cu_observed': 1413600.0, + 'fl_modeled': 2976000.0, 'fl_observed': 5952000.0, + 'mn_modeled': 1078800.0, 'mn_observed': 2157600.0, + 'se_modeled': 37200.0, 'se_observed': 74400.0, + 'vita_modeled': 1116000.0, 'vita_observed': 2232000.0, + 'betac_modeled': 5952000.0, 'betac_observed': 11904000.0, + 'alphac_modeled': 855600.0, 'alphac_observed': 1711200.0, + 'vite_modeled': 297600.0, 'vite_observed': 595200.0, + 'crypto_modeled': 595200.0, 'crypto_observed': 1190400.0, + 'lycopene_modeled': 133920.0, 'lycopene_observed': 267840.0, + 'lutein_modeled': 23436000.0, 'lutein_observed': 46872000.0, + 'betat_modeled': 186000.0, 'betat_observed': 372000.0, + 'gammat_modeled': 781200.0, 'gammat_observed': 1562400.0, + 'deltat_modeled': 706800.0, 'deltat_observed': 1413600.0, + 'vitc_modeled': 2529600.0, 'vitc_observed': 5059200.0, + 'thiamin_modeled': 148800.0, 'thiamin_observed': 297600.0, + 'riboflavin_modeled': 669600.0, 'riboflavin_observed': 1339200.0, + 'niacin_modeled': 3050400.0, 'niacin_observed': 6100800.0, + 'pantothenic_modeled': 334800.0, 'pantothenic_observed': 669600.0, + 'vitb6_modeled': 520800.0, 'vitb6_observed': 1041600.0, + 'folate_modeled': 143220000.0, 'folate_observed': 286440000.0, + 'vitb12_modeled': 744000.0, 'vitb12_observed': 1488000.0, + 'vitk_modeled': 15252000.0, 'vitk_observed': 30504000.0}, {'crop': 'soybean', 'area (ha)': 40.0, - 'production_observed': 12.0, 'production_modeled': 7.0, - 'protein_modeled': 2102100.0, 'protein_observed': 3603600.0, - 'lipid_modeled': 127400.0, 'lipid_observed': 218400.0, - 'energy_modeled': 6306300.0, 'energy_observed': 10810800.0, - 'ca_modeled': 16370900.0, 'ca_observed': 28064400.0, - 'fe_modeled': 1000090.0, 'fe_observed': 1714440.0, - 'mg_modeled': 17836000.0, 'mg_observed': 30576000.0, - 'ph_modeled': 44844800.0, 'ph_observed': 76876800.0, - 'k_modeled': 12548900.0, 'k_observed': 21512400.0, - 'na_modeled': 127400.0, 'na_observed': 218400.0, - 'zn_modeled': 312130.0, 'zn_observed': 535080.0, - 'cu_modeled': 101920.0, 'cu_observed': 174720.0, - 'fl_modeled': 191100.0, 'fl_observed': 327600.0, - 'mn_modeled': 331240.0, 'mn_observed': 567840.0, - 'se_modeled': 19110.0, 'se_observed': 32760.0, - 'vita_modeled': 191100.0, 'vita_observed': 327600.0, - 'betac_modeled': 1019200.0, 'betac_observed': 1747200.0, - 'alphac_modeled': 63700.0, 'alphac_observed': 109200.0, - 'vite_modeled': 50960.0, 'vite_observed': 87360.0, - 'crypto_modeled': 38220.0, 'crypto_observed': 65520.0, - 'lycopene_modeled': 19110.0, 'lycopene_observed': 32760.0, - 'lutein_modeled': 3885700.0, 'lutein_observed': 6661200.0, - 'betat_modeled': 31850.0, 'betat_observed': 54600.0, - 'gammat_modeled': 146510.0, 'gammat_observed': 251160.0, - 'deltat_modeled': 76440.0, 'deltat_observed': 131040.0, - 'vitc_modeled': 191100.0, 'vitc_observed': 327600.0, - 'thiamin_modeled': 26754.0, 'thiamin_observed': 45864.0, - 'riboflavin_modeled': 52234.0, 'riboflavin_observed': 89544.0, - 'niacin_modeled': 777140.0, 'niacin_observed': 1332240.0, - 'pantothenic_modeled': 58604.0, 'pantothenic_observed': 100464.0, - 'vitb6_modeled': 343980.0, 'vitb6_observed': 589680.0, - 'folate_modeled': 19428500.0, 'folate_observed': 33306000.0, - 'vitb12_modeled': 191100.0, 'vitb12_observed': 327600.0, - 'vitk_modeled': 2675400.0, 'vitk_observed': 4586400.0}]) + 'production_observed': 120.0, 'production_modeled': 70.0, + 'protein_modeled': 21021000.0, 'protein_observed': 36036000.0, + 'lipid_modeled': 1274000.0, 'lipid_observed': 2184000.0, + 'energy_modeled': 63063000.0, 'energy_observed': 108108000.0, + 'ca_modeled': 163709000.0, 'ca_observed': 280644000.0, + 'fe_modeled': 10000900.0, 'fe_observed': 17144400.0, + 'mg_modeled': 178360000.0, 'mg_observed': 305760000.0, + 'ph_modeled': 448448000.0, 'ph_observed': 768768000.0, + 'k_modeled': 125489000.0, 'k_observed': 215124000.0, + 'na_modeled': 1274000.0, 'na_observed': 2184000.0, + 'zn_modeled': 3121300.0, 'zn_observed': 5350800.0, + 'cu_modeled': 1019200.0, 'cu_observed': 1747200.0, + 'fl_modeled': 1911000.0, 'fl_observed': 3276000.0, + 'mn_modeled': 3312400.0, 'mn_observed': 5678400.0, + 'se_modeled': 191100.0, 'se_observed': 327600.0, + 'vita_modeled': 1911000.0, 'vita_observed': 3276000.0, + 'betac_modeled': 10192000.0, 'betac_observed': 17472000.0, + 'alphac_modeled': 637000.0, 'alphac_observed': 1092000.0, + 'vite_modeled': 509600.0, 'vite_observed': 873600.0, + 'crypto_modeled': 382200.0, 'crypto_observed': 655200.0, + 'lycopene_modeled': 191100.0, 'lycopene_observed': 327600.0, + 'lutein_modeled': 38857000.0, 'lutein_observed': 66612000.0, + 'betat_modeled': 318500.0, 'betat_observed': 546000.0, + 'gammat_modeled': 1465100.0, 'gammat_observed': 2511600.0, + 'deltat_modeled': 764400.0, 'deltat_observed': 1310400.0, + 'vitc_modeled': 1911000.0, 'vitc_observed': 3276000.0, + 'thiamin_modeled': 267540.0, 'thiamin_observed': 458640.0, + 'riboflavin_modeled': 522340.0, 'riboflavin_observed': 895440.0, + 'niacin_modeled': 7771400.0, 'niacin_observed': 13322400.0, + 'pantothenic_modeled': 586040.0, 'pantothenic_observed': 1004640.0, + 'vitb6_modeled': 3439800.0, 'vitb6_observed': 5896800.0, + 'folate_modeled': 194285000.0, 'folate_observed': 333060000.0, + 'vitb12_modeled': 1911000.0, 'vitb12_observed': 3276000.0, + 'vitk_modeled': 26754000.0, 'vitk_observed': 45864000.0}]) nutrient_df = create_nutrient_df() @@ -675,76 +743,76 @@ def _create_expected_agg_table(): """Create expected output results""" # Define the new values manually return pandas.DataFrame([ - {"FID": 0, "corn_modeled": 1, "corn_observed": 4, - "soybean_modeled": 2, "soybean_observed": 5, - "protein_modeled": 991200, "protein_observed": 3063900, - "lipid_modeled": 110800, "lipid_observed": 388600, - "energy_modeled": 6228600, "energy_observed": 22211700, - "ca_modeled": 4928500, "ca_observed": 12697900, - "fe_modeled": 431750, "fe_observed": 1298390, - "mg_modeled": 7700000, "mg_observed": 23156000, - "ph_modeled": 19360000, "ph_observed": 58220800, - "k_modeled": 19646500, "k_observed": 73207900, - "na_modeled": 55000, "na_observed": 165400, - "zn_modeled": 134750, "zn_observed": 405230, - "cu_modeled": 46790, "cu_observed": 143480, - "fl_modeled": 129000, "fl_observed": 434100, - "mn_modeled": 121610, "mn_observed": 344480, - "se_modeled": 6390, "se_observed": 17370, - "vita_modeled": 82500, "vita_observed": 248100, - "betac_modeled": 440000, "betac_observed": 1323200, - "alphac_modeled": 39590, "alphac_observed": 131060, - "vite_modeled": 22000, "vite_observed": 66160, - "crypto_modeled": 25800, "crypto_observed": 86820, - "lycopene_modeled": 8808, "lycopene_observed": 27042, - "lutein_modeled": 1696100, "lutein_observed": 5119100, - "betat_modeled": 13750, "betat_observed": 41350, - "gammat_modeled": 61390, "gammat_observed": 182770, - "deltat_modeled": 39510, "deltat_observed": 125280, - "vitc_modeled": 117840, "vitc_observed": 389460, - "thiamin_modeled": 11364, "thiamin_observed": 33990, - "riboflavin_modeled": 31664, "riboflavin_observed": 104270, - "niacin_modeled": 298300, "niacin_observed": 860140, - "pantothenic_modeled": 25114, "pantothenic_observed": 75340, - "vitb6_modeled": 111300, "vitb6_observed": 297780, - "folate_modeled": 9131500, "folate_observed": 28199500, - "vitb12_modeled": 73200, "vitb12_observed": 210900, - "vitk_modeled": 1145700, "vitk_observed": 3436200}, - {"FID": 1, "corn_modeled": 4, "corn_observed": 8, - "soybean_modeled": 7, "soybean_observed": 12, - "protein_modeled": 3664500, "protein_observed": 6728400, - "lipid_modeled": 425000, "lipid_observed": 813600, - "energy_modeled": 24013500, "energy_observed": 46225200, - "ca_modeled": 17375300, "ca_observed": 30073200, - "fe_modeled": 1584130, "fe_observed": 2882520, - "mg_modeled": 28252000, "mg_observed": 51408000, - "ph_modeled": 71033600, "ph_observed": 129254400, - "k_modeled": 76793300, "k_observed": 150001200, - "na_modeled": 201800, "na_observed": 367200, - "zn_modeled": 494410, "zn_observed": 899640, - "cu_modeled": 172600, "cu_observed": 316080, - "fl_modeled": 488700, "fl_observed": 922800, - "mn_modeled": 439120, "mn_observed": 783600, - "se_modeled": 22830, "se_observed": 40200, - "vita_modeled": 302700, "vita_observed": 550800, - "betac_modeled": 1614400, "betac_observed": 2937600, - "alphac_modeled": 149260, "alphac_observed": 280320, - "vite_modeled": 80720, "vite_observed": 146880, - "crypto_modeled": 97740, "crypto_observed": 184560, - "lycopene_modeled": 32502, "lycopene_observed": 59544, - "lutein_modeled": 6229300, "lutein_observed": 11348400, - "betat_modeled": 50450, "betat_observed": 91800, - "gammat_modeled": 224630, "gammat_observed": 407400, - "deltat_modeled": 147120, "deltat_observed": 272400, - "vitc_modeled": 444060, "vitc_observed": 833520, - "thiamin_modeled": 41634, "thiamin_observed": 75624, - "riboflavin_modeled": 119194, "riboflavin_observed": 223464, - "niacin_modeled": 1082180, "niacin_observed": 1942320, - "pantothenic_modeled": 92084, "pantothenic_observed": 167424, - "vitb6_modeled": 396060, "vitb6_observed": 693840, - "folate_modeled": 33750500, "folate_observed": 61950000, - "vitb12_modeled": 265500, "vitb12_observed": 476400, - "vitk_modeled": 4200600, "vitk_observed": 7636800} + {"FID": 0, "corn_modeled": 10, "corn_observed": 40, + "soybean_modeled": 20, "soybean_observed": 50, + "protein_modeled": 9912000, "protein_observed": 30639000, + "lipid_modeled": 1108000, "lipid_observed": 3886000, + "energy_modeled": 62286000, "energy_observed": 222117000, + "ca_modeled": 49285000, "ca_observed": 126979000, + "fe_modeled": 4317500, "fe_observed": 12983900, + "mg_modeled": 77000000, "mg_observed": 231560000, + "ph_modeled": 193600000, "ph_observed": 582208000, + "k_modeled": 196465000, "k_observed": 732079000, + "na_modeled": 550000, "na_observed": 1654000, + "zn_modeled": 1347500, "zn_observed": 4052300, + "cu_modeled": 467900, "cu_observed": 1434800, + "fl_modeled": 1290000, "fl_observed": 4341000, + "mn_modeled": 1216100, "mn_observed": 3444800, + "se_modeled": 63900, "se_observed": 173700, + "vita_modeled": 825000, "vita_observed": 2481000, + "betac_modeled": 4400000, "betac_observed": 13232000, + "alphac_modeled": 395900, "alphac_observed": 1310600, + "vite_modeled": 220000, "vite_observed": 661600, + "crypto_modeled": 258000, "crypto_observed": 868200, + "lycopene_modeled": 88080, "lycopene_observed": 270420, + "lutein_modeled": 16961000, "lutein_observed": 51191000, + "betat_modeled": 137500, "betat_observed": 413500, + "gammat_modeled": 613900, "gammat_observed": 1827700, + "deltat_modeled": 395100, "deltat_observed": 1252800, + "vitc_modeled": 1178400, "vitc_observed": 3894600, + "thiamin_modeled": 113640, "thiamin_observed": 339900, + "riboflavin_modeled": 316640, "riboflavin_observed": 1042700, + "niacin_modeled": 2983000, "niacin_observed": 8601400, + "pantothenic_modeled": 251140, "pantothenic_observed": 753400, + "vitb6_modeled": 1113000, "vitb6_observed": 2977800, + "folate_modeled": 91315000, "folate_observed": 281995000, + "vitb12_modeled": 732000, "vitb12_observed": 2109000, + "vitk_modeled": 11457000, "vitk_observed": 34362000}, + {"FID": 1, "corn_modeled": 40, "corn_observed": 80, + "soybean_modeled": 70, "soybean_observed": 120, + "protein_modeled": 36645000, "protein_observed": 67284000, + "lipid_modeled": 4250000, "lipid_observed": 8136000, + "energy_modeled": 240135000, "energy_observed": 462252000, + "ca_modeled": 173753000, "ca_observed": 300732000, + "fe_modeled": 15841300, "fe_observed": 28825200, + "mg_modeled": 282520000, "mg_observed": 514080000, + "ph_modeled": 710336000, "ph_observed": 1292544000, + "k_modeled": 767933000, "k_observed": 1500012000, + "na_modeled": 2018000, "na_observed": 3672000, + "zn_modeled": 4944100, "zn_observed": 8996400, + "cu_modeled": 1726000, "cu_observed": 3160800, + "fl_modeled": 4887000, "fl_observed": 9228000, + "mn_modeled": 4391200, "mn_observed": 7836000, + "se_modeled": 228300, "se_observed": 402000, + "vita_modeled": 3027000, "vita_observed": 5508000, + "betac_modeled": 16144000, "betac_observed": 29376000, + "alphac_modeled": 1492600, "alphac_observed": 2803200, + "vite_modeled": 807200, "vite_observed": 1468800, + "crypto_modeled": 977400, "crypto_observed": 1845600, + "lycopene_modeled": 325020, "lycopene_observed": 595440, + "lutein_modeled": 62293000, "lutein_observed": 113484000, + "betat_modeled": 504500, "betat_observed": 918000, + "gammat_modeled": 2246300, "gammat_observed": 4074000, + "deltat_modeled": 1471200, "deltat_observed": 2724000, + "vitc_modeled": 4440600, "vitc_observed": 8335200, + "thiamin_modeled": 416340, "thiamin_observed": 756240, + "riboflavin_modeled": 1191940, "riboflavin_observed": 2234640, + "niacin_modeled": 10821800, "niacin_observed": 19423200, + "pantothenic_modeled": 920840, "pantothenic_observed": 1674240, + "vitb6_modeled": 3960600, "vitb6_observed": 6938400, + "folate_modeled": 337505000, "folate_observed": 619500000, + "vitb12_modeled": 2655000, "vitb12_observed": 4764000, + "vitk_modeled": 42006000, "vitk_observed": 76368000} ], dtype=float) workspace = self.workspace_dir @@ -765,25 +833,24 @@ def _create_expected_agg_table(): output_dir = os.path.join(workspace, "OUTPUT") os.makedirs(output_dir, exist_ok=True) file_suffix = 'test' - target_aggregate_table_path = '' # unused - - _create_crop_rasters(output_dir, crop_names, file_suffix) - - aggregate_regression_results_to_polygons( - base_aggregate_vector_path, target_aggregate_vector_path, - landcover_raster_projection, crop_names, - nutrient_df, output_dir, file_suffix, - target_aggregate_table_path) _AGGREGATE_TABLE_FILE_PATTERN = os.path.join( - '.','aggregate_results%s.csv') + '.', 'aggregate_results%s.csv') aggregate_table_path = os.path.join( output_dir, _AGGREGATE_TABLE_FILE_PATTERN % file_suffix) + pixel_area_ha = 10 + + _create_crop_rasters(output_dir, crop_names, file_suffix) + + aggregate_regression_results_to_polygons( + base_aggregate_vector_path, target_aggregate_vector_path, + aggregate_table_path, landcover_raster_projection, crop_names, + nutrient_df, pixel_area_ha, output_dir, file_suffix) + actual_aggregate_table = pandas.read_csv(aggregate_table_path, dtype=float) - print(actual_aggregate_table) expected_aggregate_table = _create_expected_agg_table() @@ -791,7 +858,6 @@ def _create_expected_agg_table(): actual_aggregate_table, expected_aggregate_table) - class CropValidationTests(unittest.TestCase): """Tests for the Crop Productions' MODEL_SPEC and validation.""" diff --git a/tests/test_forest_carbon_edge.py b/tests/test_forest_carbon_edge.py index cacbd13ae6..dd8f12f162 100644 --- a/tests/test_forest_carbon_edge.py +++ b/tests/test_forest_carbon_edge.py @@ -7,6 +7,8 @@ from osgeo import gdal import numpy +import pygeoprocessing + gdal.UseExceptions() REGRESSION_DATA = os.path.join( os.path.dirname(__file__), '..', 'data', 'invest-test-data', @@ -53,21 +55,31 @@ def test_carbon_full(self): args['workspace_dir']) self._assert_vector_results_close( - args['workspace_dir'], 'id', ['c_sum', 'c_ha_mean'], os.path.join( + 'id', + ['c_sum', 'c_ha_mean'], + os.path.join( args['workspace_dir'], 'aggregated_carbon_stocks.shp'), os.path.join(REGRESSION_DATA, 'agg_results_base.shp')) - expected_carbon_raster = gdal.OpenEx(os.path.join(REGRESSION_DATA, - 'carbon_map.tif')) - expected_carbon_band = expected_carbon_raster.GetRasterBand(1) - expected_carbon_array = expected_carbon_band.ReadAsArray() - actual_carbon_raster = gdal.OpenEx(os.path.join(REGRESSION_DATA, - 'carbon_map.tif')) - actual_carbon_band = actual_carbon_raster.GetRasterBand(1) - actual_carbon_array = actual_carbon_band.ReadAsArray() - self.assertTrue(numpy.allclose(expected_carbon_array, - actual_carbon_array)) - + # Check raster output to make sure values are in Mg/ha. + raster_path = os.path.join(args['workspace_dir'], 'carbon_map.tif') + raster_info = pygeoprocessing.get_raster_info(raster_path) + nodata = raster_info['nodata'][0] + raster_sum = 0.0 + for _, block in pygeoprocessing.iterblocks((raster_path, 1)): + raster_sum += numpy.sum( + block[~pygeoprocessing.array_equals_nodata( + block, nodata)], dtype=numpy.float64) + + # expected_sum_per_pixel_values is in Mg, calculated from raster + # generated by the model when each pixel value was in Mg/px. + # Since pixel values are now Mg/ha, raster sum is (Mg•px)/ha. + # To convert expected_sum_per_pixel_values from Mg, multiply by px/ha. + expected_sum_per_pixel_values = 21414391.997192383 + pixel_area = abs(numpy.prod(raster_info['pixel_size'])) + pixels_per_hectare = 10000 / pixel_area + expected_sum = expected_sum_per_pixel_values * pixels_per_hectare + numpy.testing.assert_allclose(raster_sum, expected_sum) def test_carbon_dup_output(self): """Forest Carbon Edge: test for existing output overlap.""" @@ -128,7 +140,8 @@ def test_carbon_no_forest_edge(self): REGRESSION_DATA, 'file_list_no_edge_effect.txt'), args['workspace_dir']) self._assert_vector_results_close( - args['workspace_dir'], 'id', ['c_sum', 'c_ha_mean'], + 'id', + ['c_sum', 'c_ha_mean'], os.path.join( args['workspace_dir'], 'aggregated_carbon_stocks_small_no_edge_effect.shp'), @@ -162,7 +175,7 @@ def test_carbon_bad_pool_value(self): expected_message = 'Could not interpret carbon pool value' actual_message = str(cm.exception) self.assertTrue(expected_message in actual_message, actual_message) - + def test_missing_lulc_value(self): """Forest Carbon Edge: test with missing LULC value.""" from natcap.invest import forest_carbon_edge_effect @@ -261,12 +274,11 @@ def _test_same_files(base_list_path, directory_path): '\n'.join(missing_files)) def _assert_vector_results_close( - self, workspace_dir, id_fieldname, field_list, result_vector_path, + self, id_fieldname, field_list, result_vector_path, expected_vector_path): """Test workspace state against expected aggregate results. Args: - workspace_dir (string): path to the completed model workspace id_fieldname (string): fieldname of the unique ID. field_list (list of string): list of fields to check near-equality. diff --git a/tests/test_ndr.py b/tests/test_ndr.py index f0ce7d780c..f70bdb9cad 100644 --- a/tests/test_ndr.py +++ b/tests/test_ndr.py @@ -1,5 +1,4 @@ """InVEST NDR model tests.""" -import collections import os import shutil import tempfile @@ -195,11 +194,12 @@ def test_no_nutrient_selected(self): self.assertEqual(len(validation_messages), 1) def test_base_regression(self): - """NDR base regression test on sample data. + """NDR base regression test on test data. - Execute NDR with sample data and checks that the output files are - generated and that the aggregate shapefile fields are the same as the - regression case. + Executes NDR with test data. Checks for accuracy of aggregate + values in summary vector, presence of drainage raster in + intermediate outputs, and accuracy of raster outputs (as + measured by the sum of their non-nodata pixel values). """ from natcap.invest.ndr import ndr @@ -222,14 +222,18 @@ def test_base_regression(self): mismatch_list = [] # these values were generated by manual inspection of regression # results - for field, expected_value in [ - ('p_surface_load', 41.826904), - ('p_surface_export', 5.870544), - ('n_surface_load', 2977.551270), - ('n_surface_export', 274.020844), - ('n_subsurface_load', 28.558048), - ('n_subsurface_export', 15.578484), - ('n_total_export', 289.599314)]: + expected_watershed_totals = { + 'p_surface_load': 41.826904, + 'p_surface_export': 5.870544, + 'n_surface_load': 2977.551270, + 'n_surface_export': 274.020844, + 'n_subsurface_load': 28.558048, + 'n_subsurface_export': 15.578484, + 'n_total_export': 289.599314 + } + + for field in expected_watershed_totals: + expected_value = expected_watershed_totals[field] val = result_feature.GetField(field) if not numpy.isclose(val, expected_value): mismatch_list.append( @@ -237,7 +241,7 @@ def test_base_regression(self): 'actual: %f' % val)) result_feature = None if mismatch_list: - raise RuntimeError("results not expected: %s" % mismatch_list) + raise AssertionError("results not expected: %s" % mismatch_list) # We only need to test that the drainage mask exists. Functionality # for that raster is tested in SDR. @@ -247,6 +251,28 @@ def test_base_regression(self): args['workspace_dir'], 'intermediate_outputs', 'what_drains_to_stream.tif'))) + # Check raster outputs to make sure values are in kg/ha/yr. + raster_info = pygeoprocessing.get_raster_info(args['dem_path']) + pixel_area = abs(numpy.prod(raster_info['pixel_size'])) + pixels_per_hectare = 10000 / pixel_area + for attr_name in ['p_surface_export', + 'n_surface_export', + 'n_subsurface_export', + 'n_total_export']: + # Since pixel values are kg/(ha•yr), raster sum is (kg•px)/(ha•yr), + # equal to the watershed total (kg/yr) * (pixels_per_hectare px/ha). + expected_sum = (expected_watershed_totals[attr_name] + * pixels_per_hectare) + raster_name = attr_name + '.tif' + raster_path = os.path.join(args['workspace_dir'], raster_name) + nodata = pygeoprocessing.get_raster_info(raster_path)['nodata'][0] + raster_sum = 0.0 + for _, block in pygeoprocessing.iterblocks((raster_path, 1)): + raster_sum += numpy.sum( + block[~pygeoprocessing.array_equals_nodata( + block, nodata)], dtype=numpy.float64) + numpy.testing.assert_allclose(raster_sum, expected_sum, rtol=1e-6) + def test_regression_undefined_nodata(self): """NDR test when DEM, LULC and runoff proxy have undefined nodata.""" from natcap.invest.ndr import ndr diff --git a/tests/test_sdr.py b/tests/test_sdr.py index 341076a085..0b517b3953 100644 --- a/tests/test_sdr.py +++ b/tests/test_sdr.py @@ -128,20 +128,21 @@ def test_sdr_validation_key_no_value(self): 'expected a validation error but didn\'t get one') def test_base_regression(self): - """SDR base regression test on sample data. + """SDR base regression test on test data. - Execute SDR with sample data and checks that the output files are - generated and that the aggregate shapefile fields are the same as the - regression case. + Executes SDR with test data. Checks for accuracy of aggregate + values in summary vector, presence of drainage raster in + intermediate outputs, absence of negative (non-nodata) values + in sed_deposition raster, and accuracy of raster outputs (as + measured by the sum of their non-nodata pixel values). """ from natcap.invest.sdr import sdr # use predefined directory so test can clean up files during teardown args = SDRTests.generate_base_args(self.workspace_dir) - # make args explicit that this is a base run of SWY sdr.execute(args) - expected_results = { + expected_watershed_totals = { 'usle_tot': 2.62457418442, 'sed_export': 0.09748090804, 'sed_dep': 1.71672844887, @@ -151,7 +152,8 @@ def test_base_regression(self): vector_path = os.path.join( args['workspace_dir'], 'watershed_results_sdr.shp') - assert_expected_results_in_vector(expected_results, vector_path) + assert_expected_results_in_vector(expected_watershed_totals, + vector_path) # We only need to test that the drainage mask exists. Functionality # for that raster is tested elsewhere @@ -175,6 +177,29 @@ def test_base_regression(self): self.assertEqual( numpy.count_nonzero(sed_dep_array[negative_non_nodata_mask]), 0) + # Check raster outputs to make sure values are in Mg/ha/yr. + raster_info = pygeoprocessing.get_raster_info(args['dem_path']) + pixel_area = abs(numpy.prod(raster_info['pixel_size'])) + pixels_per_hectare = 10000 / pixel_area + for (raster_name, + attr_name) in [('usle.tif', 'usle_tot'), + ('sed_export.tif', 'sed_export'), + ('sed_deposition.tif', 'sed_dep'), + ('avoided_export.tif', 'avoid_exp'), + ('avoided_erosion.tif', 'avoid_eros')]: + # Since pixel values are Mg/(ha•yr), raster sum is (Mg•px)/(ha•yr), + # equal to the watershed total (Mg/yr) * (pixels_per_hectare px/ha). + expected_sum = (expected_watershed_totals[attr_name] + * pixels_per_hectare) + raster_path = os.path.join(args['workspace_dir'], raster_name) + nodata = pygeoprocessing.get_raster_info(raster_path)['nodata'][0] + raster_sum = 0.0 + for _, block in pygeoprocessing.iterblocks((raster_path, 1)): + raster_sum += numpy.sum( + block[~pygeoprocessing.array_equals_nodata( + block, nodata)], dtype=numpy.float64) + numpy.testing.assert_allclose(raster_sum, expected_sum) + def test_regression_with_undefined_nodata(self): """SDR base regression test with undefined nodata values. @@ -288,8 +313,8 @@ def test_base_usle_c_too_large(self): with self.assertRaises(ValueError) as context: sdr.execute(args) self.assertIn( - f'A value in the biophysical table is not a number ' - f'within range 0..1.', str(context.exception)) + 'A value in the biophysical table is not a number ' + 'within range 0..1.', str(context.exception)) def test_base_usle_p_nan(self): """SDR test expected exception for USLE_P not a number.""" @@ -304,7 +329,7 @@ def test_base_usle_p_nan(self): with self.assertRaises(ValueError) as context: sdr.execute(args) self.assertIn( - f'could not be interpreted as ratios', str(context.exception)) + 'could not be interpreted as ratios', str(context.exception)) def test_lucode_not_a_number(self): """SDR test expected exception for invalid data in lucode column."""