diff --git a/HISTORY.rst b/HISTORY.rst index 4554c375cf..cf941dccf2 100644 --- a/HISTORY.rst +++ b/HISTORY.rst @@ -79,6 +79,12 @@ Unreleased Changes * Replaced custom kernel implementation with ``pygeoprocessing.kernels``. Convolution results may be slightly different (more accurate). * SDR + * Fixed an issue with SDR's sediment deposition where large regions would + become nodata in cases where the DEM has valid data but other inputs + (LULC, erosivity, erodibility) did not have valid pixels. Now, all + raster inputs are mutually masked so that only those pixel stacks + continue through to the model where all pixels in the stack are + non-nodata. (`#911 `_) * RKLS, USLE, avoided erosion, and avoided export rasters will now have nodata in streams (`#1415 `_) * Fixed an issue in SDR's sediment deposition where, on rasters with more diff --git a/src/natcap/invest/sdr/sdr.py b/src/natcap/invest/sdr/sdr.py index 9a8ad2d66b..f23526b4a1 100644 --- a/src/natcap/invest/sdr/sdr.py +++ b/src/natcap/invest/sdr/sdr.py @@ -390,12 +390,62 @@ }, "aligned_lulc.tif": { "about": gettext( - "Copy of the input drainage map, clipped to " + "Copy of the input Land Use Land Cover map, clipped to " "the extent of the other raster inputs and " "aligned to the DEM."), "bands": {1: {"type": "integer"}}, + }, + "mask.tif": { + "about": gettext( + "A raster aligned to the DEM and clipped to the " + "extent of the other raster inputs. Pixel values " + "indicate where a nodata value exists in the stack " + "of aligned rasters (pixel value of 0), or if all " + "values in the stack of rasters at this pixel " + "location are valid."), + "bands": {1: {"type": "integer"}}, + }, + "masked_dem.tif": { + "about": gettext( + "A copy of the aligned DEM, masked using the mask raster." + ), + "bands": {1: { + "type": "number", + "units": u.meter}}, + }, + "masked_drainage.tif": { + "about": gettext( + "A copy of the aligned drainage map, masked using the " + "mask raster." + ), + "bands": {1: {"type": "integer"}} + }, + "masked_erodibility.tif": { + "about": gettext( + "A copy of the aligned erodibility map, masked using " + "the mask raster." + ), + "bands": {1: { + "type": "number", + "units": u.metric_ton*u.hectare*u.hour/(u.hectare*u.megajoule*u.millimeter) + }}, + }, + "masked_erosivity.tif": { + "about": gettext( + "A copy of the aligned erosivity map, masked using " + "the mask raster."), + "bands": {1: { + "type": "number", + "units": u.megajoule*u.millimeter/(u.hectare*u.hour*u.year) + }}, + }, + "masked_lulc.tif": { + "about": gettext( + "A copy of the aligned Land Use Land Cover map, " + "masked using the mask raster."), + "bands": {1: {"type": "integer"}}, } - } + }, }, "taskgraph_cache": spec_utils.TASKGRAPH_DIR } @@ -421,6 +471,12 @@ 'aligned_erodibility_path': 'aligned_erodibility.tif', 'aligned_erosivity_path': 'aligned_erosivity.tif', 'aligned_lulc_path': 'aligned_lulc.tif', + 'mask_path': 'mask.tif', + 'masked_dem_path': 'masked_dem.tif', + 'masked_drainage_path': 'masked_drainage.tif', + 'masked_erodibility_path': 'masked_erodibility.tif', + 'masked_erosivity_path': 'masked_erosivity.tif', + 'masked_lulc_path': 'masked_lulc.tif', 'cp_factor_path': 'cp.tif', 'd_dn_path': 'd_dn.tif', 'd_up_path': 'd_up.tif', @@ -530,17 +586,22 @@ def execute(args): base_list = [] aligned_list = [] - for file_key in ['dem', 'lulc', 'erosivity', 'erodibility']: - base_list.append(args[file_key + "_path"]) - aligned_list.append(f_reg["aligned_" + file_key + "_path"]) - # all continuous rasters can use bilinaer, but lulc should be mode + masked_list = [] + input_raster_key_list = ['dem', 'lulc', 'erosivity', 'erodibility'] + for file_key in input_raster_key_list: + base_list.append(args[f"{file_key}_path"]) + aligned_list.append(f_reg[f"aligned_{file_key}_path"]) + masked_list.append(f_reg[f"masked_{file_key}_path"]) + # all continuous rasters can use bilinear, but lulc should be mode interpolation_list = ['bilinear', 'mode', 'bilinear', 'bilinear'] drainage_present = False if 'drainage_path' in args and args['drainage_path'] != '': drainage_present = True + input_raster_key_list.append('drainage') base_list.append(args['drainage_path']) aligned_list.append(f_reg['aligned_drainage_path']) + masked_list.append(f_reg['masked_drainage_path']) interpolation_list.append('near') dem_raster_info = pygeoprocessing.get_raster_info(args['dem_path']) @@ -563,13 +624,39 @@ def execute(args): target_path_list=aligned_list, task_name='align input rasters') + mutual_mask_task = task_graph.add_task( + func=pygeoprocessing.raster_map, + kwargs={ + 'op': _create_mutual_mask_op, + 'rasters': aligned_list, + 'target_path': f_reg['mask_path'], + 'target_nodata': 0, + }, + target_path_list=[f_reg['mask_path']], + dependent_task_list=[align_task], + task_name='create mask') + + mask_tasks = {} # use a dict so we can put these in a loop + for key, aligned_path, masked_path in zip(input_raster_key_list, + aligned_list, masked_list): + mask_tasks[f"masked_{key}"] = task_graph.add_task( + func=pygeoprocessing.raster_map, + kwargs={ + 'op': _mask_single_raster_op, + 'rasters': [aligned_path, f_reg['mask_path']], + 'target_path': masked_path, + }, + target_path_list=[masked_path], + dependent_task_list=[mutual_mask_task, align_task], + task_name=f'mask {key}') + pit_fill_task = task_graph.add_task( func=pygeoprocessing.routing.fill_pits, args=( - (f_reg['aligned_dem_path'], 1), + (f_reg['masked_dem_path'], 1), f_reg['pit_filled_dem_path']), target_path_list=[f_reg['pit_filled_dem_path']], - dependent_task_list=[align_task], + dependent_task_list=[mask_tasks['masked_dem']], task_name='fill pits') slope_task = task_graph.add_task( @@ -638,11 +725,11 @@ def execute(args): func=pygeoprocessing.raster_map, kwargs=dict( op=add_drainage_op, - rasters=[f_reg['stream_path'], f_reg['aligned_drainage_path']], + rasters=[f_reg['stream_path'], f_reg['masked_drainage_path']], target_path=f_reg['stream_and_drainage_path'], target_dtype=numpy.uint8), target_path_list=[f_reg['stream_and_drainage_path']], - dependent_task_list=[stream_task, align_task], + dependent_task_list=[stream_task, mask_tasks['masked_drainage']], task_name='add drainage') drainage_raster_path_task = ( f_reg['stream_and_drainage_path'], drainage_task) @@ -654,33 +741,34 @@ def execute(args): threshold_w_task = task_graph.add_task( func=_calculate_w, args=( - lulc_to_c, f_reg['aligned_lulc_path'], f_reg['w_path'], + lulc_to_c, f_reg['masked_lulc_path'], f_reg['w_path'], f_reg['thresholded_w_path']), target_path_list=[f_reg['w_path'], f_reg['thresholded_w_path']], - dependent_task_list=[align_task], + dependent_task_list=[mask_tasks['masked_lulc']], task_name='calculate W') lulc_to_cp = (biophysical_df['usle_c'] * biophysical_df['usle_p']).to_dict() cp_task = task_graph.add_task( func=_calculate_cp, args=( - lulc_to_cp, f_reg['aligned_lulc_path'], + lulc_to_cp, f_reg['masked_lulc_path'], f_reg['cp_factor_path']), target_path_list=[f_reg['cp_factor_path']], - dependent_task_list=[align_task], + dependent_task_list=[mask_tasks['masked_lulc']], task_name='calculate CP') rkls_task = task_graph.add_task( func=_calculate_rkls, args=( f_reg['ls_path'], - f_reg['aligned_erosivity_path'], - f_reg['aligned_erodibility_path'], + f_reg['masked_erosivity_path'], + f_reg['masked_erodibility_path'], drainage_raster_path_task[0], f_reg['rkls_path']), target_path_list=[f_reg['rkls_path']], dependent_task_list=[ - align_task, drainage_raster_path_task[1], ls_factor_task], + mask_tasks['masked_erosivity'], mask_tasks['masked_erodibility'], + drainage_raster_path_task[1], ls_factor_task], task_name='calculate RKLS') usle_task = task_graph.add_task( @@ -711,8 +799,7 @@ def execute(args): accumulation_path, out_bar_path), target_path_list=[accumulation_path, out_bar_path], dependent_task_list=[ - align_task, factor_task, flow_accumulation_task, - flow_dir_task], + factor_task, flow_accumulation_task, flow_dir_task], task_name='calculate %s' % bar_id) bar_task_map[bar_id] = bar_task @@ -846,6 +933,14 @@ def execute(args): task_graph.join() +# raster_map op for building a mask where all pixels in the stack are valid. +def _create_mutual_mask_op(*arrays): return 1 + + +# raster_map op for using a mask raster to mask out another raster. +def _mask_single_raster_op(source_array, mask_array): return source_array + + def _avoided_export_op(avoided_erosion, sdr, sed_deposition): """raster_map equation: calculate total retention. diff --git a/tests/test_sdr.py b/tests/test_sdr.py index 3c324934b0..e509837eca 100644 --- a/tests/test_sdr.py +++ b/tests/test_sdr.py @@ -239,7 +239,7 @@ def test_non_square_dem(self): expected_results = { 'sed_export': 0.08896198869, - 'usle_tot': 1.86480903625, + 'usle_tot': 1.86480891705, 'avoid_exp': 9203.955078125, 'avoid_eros': 194212.28125, }