From d6cb761699055f45540860610e77dfb0bf1e11ab Mon Sep 17 00:00:00 2001 From: James Douglass Date: Fri, 20 Oct 2023 20:18:12 -0700 Subject: [PATCH 1/8] Taking a first stab at implementing mutual masking. RE:#1433 --- src/natcap/invest/sdr/sdr.py | 118 ++++++++++++++++++++++++++++++----- 1 file changed, 103 insertions(+), 15 deletions(-) diff --git a/src/natcap/invest/sdr/sdr.py b/src/natcap/invest/sdr/sdr.py index 6eb954bd5..42756b68e 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', @@ -562,13 +618,45 @@ 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': lambda *x: 1, # any valid pixel gets a value of 1 + 'rasters': [f_reg['aligned_dem_path'], + f_reg['aligned_drainage_path'], + f_reg['aligned_erodibility_path'], + f_reg['aligned_erosivity_path'], + f_reg['aligned_lulc_path']], + '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 f_reg_key in ('dem_path', 'drainage_path', 'erodibility_path', + 'erosivity_path', 'lulc_path'): + aligned_key = f'aligned_{f_reg_key}' + masked_key = f'masked_{f_reg_key}' + mask_tasks[masked_key.replace('_path', '')] = task_graph.add_task( + func=pygeoprocessing.raster_map, + kwargs={ + 'op': lambda array, mask: array, + 'rasters': [f_reg[aligned_key], f_reg['mask_path']], + 'target_path': f_reg[masked_key], + }, + target_path_list=[f_reg[masked_key]], + dependent_task_list=[mutual_mask_task, align_task], + task_name=f'mask {f_reg_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( @@ -633,10 +721,10 @@ def execute(args): drainage_task = task_graph.add_task( func=_add_drainage( f_reg['stream_path'], - f_reg['aligned_drainage_path'], + f_reg['masked_drainage_path'], f_reg['stream_and_drainage_path']), 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) @@ -648,33 +736,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( @@ -705,8 +794,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 From 51eaed1a9020ba2aa190eebfa56dd430e57908bf Mon Sep 17 00:00:00 2001 From: James Douglass Date: Mon, 23 Oct 2023 09:45:00 -0700 Subject: [PATCH 2/8] Reworking mutual masking to accommodate drainage. RE:#1433 --- src/natcap/invest/sdr/sdr.py | 35 +++++++++++++++++------------------ 1 file changed, 17 insertions(+), 18 deletions(-) diff --git a/src/natcap/invest/sdr/sdr.py b/src/natcap/invest/sdr/sdr.py index 42756b68e..5e174e5ac 100644 --- a/src/natcap/invest/sdr/sdr.py +++ b/src/natcap/invest/sdr/sdr.py @@ -585,17 +585,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']) @@ -622,11 +627,7 @@ def execute(args): func=pygeoprocessing.raster_map, kwargs={ 'op': lambda *x: 1, # any valid pixel gets a value of 1 - 'rasters': [f_reg['aligned_dem_path'], - f_reg['aligned_drainage_path'], - f_reg['aligned_erodibility_path'], - f_reg['aligned_erosivity_path'], - f_reg['aligned_lulc_path']], + 'rasters': aligned_list, 'target_path': f_reg['mask_path'], 'target_nodata': 0, }, @@ -635,20 +636,18 @@ def execute(args): task_name='create mask') mask_tasks = {} # use a dict so we can put these in a loop - for f_reg_key in ('dem_path', 'drainage_path', 'erodibility_path', - 'erosivity_path', 'lulc_path'): - aligned_key = f'aligned_{f_reg_key}' - masked_key = f'masked_{f_reg_key}' - mask_tasks[masked_key.replace('_path', '')] = task_graph.add_task( + 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': lambda array, mask: array, - 'rasters': [f_reg[aligned_key], f_reg['mask_path']], - 'target_path': f_reg[masked_key], + 'rasters': [aligned_path, f_reg['mask_path']], + 'target_path': masked_path, }, - target_path_list=[f_reg[masked_key]], + target_path_list=[masked_path], dependent_task_list=[mutual_mask_task, align_task], - task_name=f'mask {f_reg_key}') + task_name=f'mask {key}') pit_fill_task = task_graph.add_task( func=pygeoprocessing.routing.fill_pits, From 1caecc9f51a2f0e6e580394e3d62c8c2e2048b87 Mon Sep 17 00:00:00 2001 From: James Douglass Date: Mon, 23 Oct 2023 10:04:25 -0700 Subject: [PATCH 3/8] Updating regression values. RE:#1433 --- tests/test_sdr.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/test_sdr.py b/tests/test_sdr.py index 436072fec..3eabb5390 100644 --- a/tests/test_sdr.py +++ b/tests/test_sdr.py @@ -238,10 +238,10 @@ def test_non_square_dem(self): sdr.execute(args) expected_results = { - 'sed_export': 0.08896198869, - 'usle_tot': 1.86480903625, - 'avoid_exp': 9203.955078125, - 'avoid_eros': 194212.28125, + 'sed_export': 0.08894859254, + 'usle_tot': 1.86440658569, + 'avoid_exp': 9202.60546875, + 'avoid_eros': 194171.578125, } vector_path = os.path.join( From b23a6eab5392d9047256dca8715171ce5ed0dfdb Mon Sep 17 00:00:00 2001 From: James Douglass Date: Mon, 23 Oct 2023 10:13:57 -0700 Subject: [PATCH 4/8] Noting change in HISTORY. RE#1433 --- HISTORY.rst | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/HISTORY.rst b/HISTORY.rst index 59523f494..fa2e45a23 100644 --- a/HISTORY.rst +++ b/HISTORY.rst @@ -50,6 +50,12 @@ Unreleased Changes condition when run with ``n_workers > 0``. https://github.com/natcap/invest/issues/1426 * 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 https://github.com/natcap/invest/issues/911`_) * RKLS, USLE, avoided erosion, and avoided export rasters will now have nodata in streams (`#1415 `_) * Wind Energy From 6743ff2a24f9e7f2499e1f46fe6da14575c0cc64 Mon Sep 17 00:00:00 2001 From: James Douglass Date: Mon, 23 Oct 2023 10:28:21 -0700 Subject: [PATCH 5/8] Correcting RST syntax. RE:#1433 --- HISTORY.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/HISTORY.rst b/HISTORY.rst index fa2e45a23..74ee40c4a 100644 --- a/HISTORY.rst +++ b/HISTORY.rst @@ -55,7 +55,7 @@ Unreleased Changes (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 https://github.com/natcap/invest/issues/911`_) + non-nodata. (`#911 `_) * RKLS, USLE, avoided erosion, and avoided export rasters will now have nodata in streams (`#1415 `_) * Wind Energy From 2af4a47cc39df75ac91cbe18d6c9ea836d588a82 Mon Sep 17 00:00:00 2001 From: James Douglass Date: Fri, 3 Nov 2023 17:19:03 -0700 Subject: [PATCH 6/8] Updating SDR regression values where affected by mutual-masking. RE:#1433 --- tests/test_sdr.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/test_sdr.py b/tests/test_sdr.py index 13c8ab477..e509837ec 100644 --- a/tests/test_sdr.py +++ b/tests/test_sdr.py @@ -238,10 +238,10 @@ def test_non_square_dem(self): sdr.execute(args) expected_results = { - 'sed_export': 0.08894859254, - 'usle_tot': 1.86440658569, - 'avoid_exp': 9202.60546875, - 'avoid_eros': 194171.578125, + 'sed_export': 0.08896198869, + 'usle_tot': 1.86480891705, + 'avoid_exp': 9203.955078125, + 'avoid_eros': 194212.28125, } vector_path = os.path.join( From f33cf709827dd1b426c4c7c08b8f59ca8898fa1a Mon Sep 17 00:00:00 2001 From: James Douglass Date: Mon, 6 Nov 2023 08:13:01 -0800 Subject: [PATCH 7/8] Replacing another masking lambda with a named function. RE:#1433 --- src/natcap/invest/sdr/sdr.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/natcap/invest/sdr/sdr.py b/src/natcap/invest/sdr/sdr.py index 938c02627..13802589d 100644 --- a/src/natcap/invest/sdr/sdr.py +++ b/src/natcap/invest/sdr/sdr.py @@ -641,7 +641,7 @@ def execute(args): mask_tasks[f"masked_{key}"] = task_graph.add_task( func=pygeoprocessing.raster_map, kwargs={ - 'op': lambda array, mask: array, + 'op': _mask_op, 'rasters': [aligned_path, f_reg['mask_path']], 'target_path': masked_path, }, From ee7ea1167cf01eb329cac379d083a07382d7e4bb Mon Sep 17 00:00:00 2001 From: James Douglass Date: Tue, 7 Nov 2023 10:55:23 -0800 Subject: [PATCH 8/8] Separating mutual mask from single-raster mask ops. RE:#1433 --- src/natcap/invest/sdr/sdr.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/natcap/invest/sdr/sdr.py b/src/natcap/invest/sdr/sdr.py index 13802589d..09a5786e4 100644 --- a/src/natcap/invest/sdr/sdr.py +++ b/src/natcap/invest/sdr/sdr.py @@ -626,7 +626,7 @@ def execute(args): mutual_mask_task = task_graph.add_task( func=pygeoprocessing.raster_map, kwargs={ - 'op': _mask_op, + 'op': _create_mutual_mask_op, 'rasters': aligned_list, 'target_path': f_reg['mask_path'], 'target_nodata': 0, @@ -641,7 +641,7 @@ def execute(args): mask_tasks[f"masked_{key}"] = task_graph.add_task( func=pygeoprocessing.raster_map, kwargs={ - 'op': _mask_op, + 'op': _mask_single_raster_op, 'rasters': [aligned_path, f_reg['mask_path']], 'target_path': masked_path, }, @@ -933,7 +933,11 @@ def execute(args): # raster_map op for building a mask where all pixels in the stack are valid. -def _mask_op(*arrays): return 1 +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):