Merge pull request #1907 from emlys/bugfix/1845

Bugfix: Don't use nodata pour points as real data in NDR effective retention
This commit is contained in:
James Douglass 2025-05-22 11:03:10 -07:00 committed by GitHub
commit 4f659d5236
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
3 changed files with 25 additions and 16 deletions

View File

@ -82,6 +82,13 @@ HRA
the input vectors, rather than using ``ogr.wkbUnknown``
(`#1881 <https://github.com/natcap/invest/issues/1881>`_).
NDR
===
* Fixed a bug in the effective retention calculation where nodata pour point
pixels were mistakenly used as real data. The effect of this change is most
pronounced along stream edges and should not affect the overall pattern of
results. (`#1845 <https://github.com/natcap/invest/issues/1845>`_)
3.15.1 (2025-05-06)
-------------------

View File

@ -177,6 +177,13 @@ void run_effective_retention(
neighbor.y < 0 or neighbor.y >= n_rows) {
continue;
}
neighbor_effective_retention = (
effective_retention_raster.get(
neighbor.x, neighbor.y));
if (is_close(neighbor_effective_retention, effective_retention_nodata)) {
continue;
}
if (neighbor.direction % 2 == 1) {
step_size = cell_size * 1.41421356237;
} else {
@ -189,10 +196,6 @@ void run_effective_retention(
current_step_factor = 0;
}
neighbor_effective_retention = (
effective_retention_raster.get(
neighbor.x, neighbor.y));
// Case 1: downslope neighbor is a stream pixel
if (neighbor_effective_retention == STREAM_EFFECTIVE_RETENTION) {
intermediate_retention = (
@ -222,7 +225,6 @@ void run_effective_retention(
}
}
// search upslope to see if we need to push a cell on the stack
// for i in range(8):
up_neighbors = UpslopeNeighbors<T>(Pixel<T>(flow_dir_raster, global_col, global_row));
for (auto neighbor: up_neighbors) {
neighbor_outflow_dir = INFLOW_OFFSETS[neighbor.direction];

View File

@ -152,10 +152,10 @@ class NDRTests(unittest.TestCase):
('p_surface_load', 41.826904),
('p_surface_export', 5.566120),
('n_surface_load', 2977.551270),
('n_surface_export', 274.020844),
('n_surface_export', 274.062129),
('n_subsurface_load', 28.558048),
('n_subsurface_export', 15.578484),
('n_total_export', 289.599314)]:
('n_total_export', 289.640609)]:
if not numpy.isclose(feature.GetField(field), value, atol=1e-2):
error_results[field] = (
'field', feature.GetField(field), value)
@ -226,12 +226,12 @@ class NDRTests(unittest.TestCase):
# results
expected_watershed_totals = {
'p_surface_load': 41.826904,
'p_surface_export': 5.870544,
'p_surface_export': 5.866880,
'n_surface_load': 2977.551270,
'n_surface_export': 274.020844,
'n_surface_export': 274.062129,
'n_subsurface_load': 28.558048,
'n_subsurface_export': 15.578484,
'n_total_export': 289.599314
'n_total_export': 289.640609
}
for field in expected_watershed_totals:
@ -306,12 +306,12 @@ class NDRTests(unittest.TestCase):
# results
for field, expected_value in [
('p_surface_load', 41.826904),
('p_surface_export', 4.915544),
('p_surface_export', 5.100640),
('n_surface_load', 2977.551914),
('n_surface_export', 320.082319),
('n_surface_export', 350.592891),
('n_subsurface_load', 28.558048),
('n_subsurface_export', 12.609187),
('n_total_export', 330.293407)]:
('n_total_export', 360.803969)]:
val = result_feature.GetField(field)
if not numpy.isclose(val, expected_value):
mismatch_list.append(
@ -361,12 +361,12 @@ class NDRTests(unittest.TestCase):
# results
for field, expected_value in [
('p_surface_load', 41.826904),
('p_surface_export', 5.870544),
('p_surface_export', 5.866880),
('n_surface_load', 2977.551270),
('n_surface_export', 274.020844),
('n_surface_export', 274.062129),
('n_subsurface_load', 28.558048),
('n_subsurface_export', 15.578484),
('n_total_export', 289.599314)]:
('n_total_export', 289.640609)]:
val = result_feature.GetField(field)
if not numpy.isclose(val, expected_value):
mismatch_list.append(