From 5791a659991d23ff2b474d625d71360a7e80a2cf Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Tue, 2 Jun 2026 14:56:13 -0700 Subject: [PATCH 1/2] Validate MFD fraction values in public consumers (#2873) Add _validate_mfd_fractions to utils and call it after the shape check in the seven public MFD functions: flow_accumulation_mfd, flow_length_mfd, stream_order_mfd, stream_link_mfd, flow_path_mfd, hand_mfd, watershed_mfd. The check rejects negative fractions, per-cell band sums that are neither ~1.0 nor ~0.0, and partial-NaN band patterns, raising a clear ValueError before any hydrology math runs. numpy and cupy inputs are validated eagerly; dask inputs are skipped so laziness is preserved. --- xrspatial/hydro/flow_accumulation_mfd.py | 4 + xrspatial/hydro/flow_length_mfd.py | 4 + xrspatial/hydro/flow_path_mfd.py | 4 + xrspatial/hydro/hand_mfd.py | 4 + xrspatial/hydro/stream_link_mfd.py | 4 + xrspatial/hydro/stream_order_mfd.py | 4 + .../tests/test_validate_mfd_fractions.py | 188 ++++++++++++++++++ xrspatial/hydro/watershed_mfd.py | 4 + xrspatial/utils.py | 81 ++++++++ 9 files changed, 297 insertions(+) create mode 100644 xrspatial/hydro/tests/test_validate_mfd_fractions.py diff --git a/xrspatial/hydro/flow_accumulation_mfd.py b/xrspatial/hydro/flow_accumulation_mfd.py index 18d407177..48a11f0ea 100644 --- a/xrspatial/hydro/flow_accumulation_mfd.py +++ b/xrspatial/hydro/flow_accumulation_mfd.py @@ -30,6 +30,7 @@ class cupy: # type: ignore[no-redef] da = None from xrspatial.utils import ( + _validate_mfd_fractions, _validate_raster, cuda_args, has_cuda_and_cupy, @@ -842,6 +843,9 @@ def flow_accumulation_mfd(flow_dir_mfd: xr.DataArray, f"got shape {data.shape}" ) + _validate_mfd_fractions(data, func_name='flow_accumulation_mfd', + name='flow_dir_mfd') + if isinstance(data, np.ndarray): _check_memory(data.shape[1], data.shape[2]) out = _flow_accum_mfd_cpu( diff --git a/xrspatial/hydro/flow_length_mfd.py b/xrspatial/hydro/flow_length_mfd.py index d207b81e8..f39760061 100644 --- a/xrspatial/hydro/flow_length_mfd.py +++ b/xrspatial/hydro/flow_length_mfd.py @@ -25,6 +25,7 @@ from xrspatial.hydro._boundary_store import BoundaryStore from xrspatial.utils import ( + _validate_mfd_fractions, _validate_raster, get_dataarray_resolution, has_cuda_and_cupy, @@ -1026,6 +1027,9 @@ def flow_length_mfd(flow_dir_mfd: xr.DataArray, f"got shape {data.shape}" ) + _validate_mfd_fractions(data, func_name='flow_length_mfd', + name='flow_dir_mfd') + cellsize_x, cellsize_y = get_dataarray_resolution(flow_dir_mfd) if not (np.isfinite(cellsize_x) and cellsize_x != 0 and np.isfinite(cellsize_y) and cellsize_y != 0): diff --git a/xrspatial/hydro/flow_path_mfd.py b/xrspatial/hydro/flow_path_mfd.py index 6e992e9df..d390a6cbb 100644 --- a/xrspatial/hydro/flow_path_mfd.py +++ b/xrspatial/hydro/flow_path_mfd.py @@ -29,6 +29,7 @@ class cupy: # type: ignore[no-redef] da = None from xrspatial.utils import ( + _validate_mfd_fractions, _validate_raster, has_cuda_and_cupy, is_cupy_array, @@ -428,6 +429,9 @@ def flow_path_mfd(flow_dir_mfd: xr.DataArray, raise ValueError( f"flow_dir_mfd must have shape (8, H, W), got {data.shape}") + _validate_mfd_fractions(data, func_name='flow_path_mfd', + name='flow_dir_mfd') + _, H, W = data.shape if isinstance(data, np.ndarray): diff --git a/xrspatial/hydro/hand_mfd.py b/xrspatial/hydro/hand_mfd.py index 7d199fb48..a165b0f2b 100644 --- a/xrspatial/hydro/hand_mfd.py +++ b/xrspatial/hydro/hand_mfd.py @@ -27,6 +27,7 @@ _to_numpy_f64, ) from xrspatial.utils import ( + _validate_mfd_fractions, _validate_raster, has_cuda_and_cupy, is_cupy_array, @@ -708,6 +709,9 @@ def hand_mfd(flow_dir_mfd: xr.DataArray, raise ValueError( f"flow_dir_mfd must have shape (8, H, W), got {data.shape}") + _validate_mfd_fractions(data, func_name='hand_mfd', + name='flow_dir_mfd') + _, H, W = data.shape if isinstance(data, np.ndarray): diff --git a/xrspatial/hydro/stream_link_mfd.py b/xrspatial/hydro/stream_link_mfd.py index 85004e300..58ced5370 100644 --- a/xrspatial/hydro/stream_link_mfd.py +++ b/xrspatial/hydro/stream_link_mfd.py @@ -40,6 +40,7 @@ class cupy: # type: ignore[no-redef] da = None from xrspatial.utils import ( + _validate_mfd_fractions, _validate_raster, cuda_args, has_cuda_and_cupy, @@ -1025,6 +1026,9 @@ def stream_link_mfd(fractions: xr.DataArray, f"got shape {data.shape}" ) + _validate_mfd_fractions(data, func_name='stream_link_mfd', + name='fractions') + fa_data = flow_accum.data if isinstance(data, np.ndarray): diff --git a/xrspatial/hydro/stream_order_mfd.py b/xrspatial/hydro/stream_order_mfd.py index 12e13f9c3..6e7471491 100644 --- a/xrspatial/hydro/stream_order_mfd.py +++ b/xrspatial/hydro/stream_order_mfd.py @@ -37,6 +37,7 @@ class cupy: # type: ignore[no-redef] da = None from xrspatial.utils import ( + _validate_mfd_fractions, _validate_raster, cuda_args, has_cuda_and_cupy, @@ -1512,6 +1513,9 @@ def stream_order_mfd(fractions: xr.DataArray, f"got shape {frac_data.shape}" ) + _validate_mfd_fractions(frac_data, func_name='stream_order_mfd', + name='fractions') + if isinstance(frac_data, np.ndarray): _check_memory(frac_data.shape[1], frac_data.shape[2]) frac = frac_data.astype(np.float64) diff --git a/xrspatial/hydro/tests/test_validate_mfd_fractions.py b/xrspatial/hydro/tests/test_validate_mfd_fractions.py new file mode 100644 index 000000000..a0b1b1867 --- /dev/null +++ b/xrspatial/hydro/tests/test_validate_mfd_fractions.py @@ -0,0 +1,188 @@ +"""Tests for issue #2873: MFD public APIs validate fraction VALUES. + +The public MFD functions document a fraction-grid contract: each cell's +8 bands are in [0, 1] and sum to either 1.0 (flow) or 0.0 +(pit/flat/sink), with all-NaN bands at edges and nodata cells. Before +this change they only checked the (8, H, W) shape and ran hydrology math +on whatever values they were given. + +These tests pin the value validation: negative fractions, band sums that +are neither ~1.0 nor ~0.0, and partial-NaN band patterns all raise a +clear ValueError on the in-memory (numpy / cupy) backends. Valid grids +from ``flow_direction_mfd`` still pass, and dask inputs skip the eager +value check so laziness is preserved. +""" + +import numpy as np +import pytest +import xarray as xr + +from xrspatial.hydro.flow_direction_mfd import flow_direction_mfd +from xrspatial.hydro.flow_accumulation_mfd import flow_accumulation_mfd +from xrspatial.hydro.flow_length_mfd import flow_length_mfd +from xrspatial.hydro.stream_order_mfd import stream_order_mfd +from xrspatial.hydro.stream_link_mfd import stream_link_mfd +from xrspatial.hydro.flow_path_mfd import flow_path_mfd +from xrspatial.hydro.hand_mfd import hand_mfd +from xrspatial.hydro.watershed_mfd import watershed_mfd + + +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- + +def _bowl(n=7): + """A simple bowl whose center is the lowest cell.""" + y = np.arange(n, dtype=np.float64) - n // 2 + x = np.arange(n, dtype=np.float64) - n // 2 + yy, xx = np.meshgrid(y, x, indexing='ij') + return xr.DataArray(yy ** 2 + xx ** 2, dims=['y', 'x']) + + +def _valid_mfd(n=7): + """A valid MFD fraction grid produced by flow_direction_mfd.""" + return flow_direction_mfd(_bowl(n)) + + +def _flow_accum(mfd): + return flow_accumulation_mfd(mfd) + + +def _interior_flow_cell(mfd): + """Return (r, c) of an interior cell whose bands sum to ~1.0. + + Used to inject corruption into a cell the validator will inspect. + """ + vals = mfd.values + sums = np.nansum(vals, axis=0) + nan_count = np.isnan(vals).sum(axis=0) + rows, cols = np.where((nan_count == 0) & (np.abs(sums - 1.0) <= 1e-6)) + assert len(rows) > 0, "test grid has no normal flow cell" + return int(rows[0]), int(cols[0]) + + +def _corrupt(mfd, mutate): + """Copy *mfd* and apply *mutate(values)* in place, return DataArray.""" + bad = mfd.copy(deep=True) + mutate(bad.values) + return bad + + +# Each entry: (name, callable taking the (possibly corrupt) fraction grid) +def _callers(mfd): + fa = _flow_accum(_valid_mfd()) # valid accumulation as a secondary arg + sp = xr.DataArray(np.full(mfd.shape[1:], np.nan), dims=mfd.dims[1:]) + sp.values[0, 0] = 1.0 + pp = sp + elev = _bowl() + return [ + ("flow_accumulation_mfd", lambda g: flow_accumulation_mfd(g)), + ("flow_length_mfd", lambda g: flow_length_mfd(g)), + ("stream_order_mfd", lambda g: stream_order_mfd(g, fa, threshold=1)), + ("stream_link_mfd", lambda g: stream_link_mfd(g, fa, threshold=1)), + ("flow_path_mfd", lambda g: flow_path_mfd(g, sp)), + ("hand_mfd", lambda g: hand_mfd(g, fa, elev, threshold=1)), + ("watershed_mfd", lambda g: watershed_mfd(g, pp)), + ] + + +# --------------------------------------------------------------------------- +# Valid input still passes +# --------------------------------------------------------------------------- + +class TestValidInputPasses: + def test_all_consumers_accept_valid_grid(self): + mfd = _valid_mfd() + for fname, call in _callers(mfd): + # should not raise + call(mfd) + + +# --------------------------------------------------------------------------- +# Negative fractions +# --------------------------------------------------------------------------- + +class TestNegativeFractions: + def test_each_consumer_rejects_negative(self): + mfd = _valid_mfd() + r, c = _interior_flow_cell(mfd) + + def mutate(v): + v[0, r, c] = -0.5 + v[1, r, c] += 0.5 # keep the band sum at 1.0 + + bad = _corrupt(mfd, mutate) + for fname, call in _callers(mfd): + with pytest.raises(ValueError, match="negative"): + call(bad) + + +# --------------------------------------------------------------------------- +# Band sums outside {0, 1} +# --------------------------------------------------------------------------- + +class TestBandSums: + def test_each_consumer_rejects_sum_above_one(self): + mfd = _valid_mfd() + r, c = _interior_flow_cell(mfd) + bad = _corrupt(mfd, lambda v: v.__setitem__((slice(None), r, c), 0.5)) + # 8 bands * 0.5 = 4.0 + for fname, call in _callers(mfd): + with pytest.raises(ValueError, match="sum"): + call(bad) + + def test_sink_cell_sum_zero_is_accepted(self): + # The bowl center is a pit: all 8 bands are 0.0 (sum 0.0). This + # must remain valid. + mfd = _valid_mfd() + vals = mfd.values + nan_count = np.isnan(vals).sum(axis=0) + sums = np.nansum(vals, axis=0) + has_sink = np.any((nan_count == 0) & (np.abs(sums) <= 1e-6)) + assert has_sink, "test grid has no pit/sink cell" + flow_accumulation_mfd(mfd) # should not raise + + +# --------------------------------------------------------------------------- +# Partial-NaN band pattern +# --------------------------------------------------------------------------- + +class TestPartialNaN: + def test_each_consumer_rejects_partial_nan(self): + mfd = _valid_mfd() + r, c = _interior_flow_cell(mfd) + + def mutate(v): + # NaN one band, leave the rest finite -> partial NaN + v[0, r, c] = np.nan + + bad = _corrupt(mfd, mutate) + for fname, call in _callers(mfd): + with pytest.raises(ValueError, match="partial-NaN"): + call(bad) + + def test_all_nan_cell_is_accepted(self): + # Edge cells from flow_direction_mfd are all-NaN; that is valid. + mfd = _valid_mfd() + assert np.isnan(mfd.values[:, 0, 0]).all() + flow_accumulation_mfd(mfd) # should not raise + + +# --------------------------------------------------------------------------- +# Dask skips eager value validation (laziness preserved) +# --------------------------------------------------------------------------- + +class TestDaskSkipsValueCheck: + def test_dask_invalid_values_not_validated_eagerly(self): + dask = pytest.importorskip('dask.array') + mfd = _valid_mfd() + r, c = _interior_flow_cell(mfd) + bad = _corrupt(mfd, lambda v: v.__setitem__((0, r, c), -5.0)) + dask_bad = xr.DataArray( + dask.from_array(bad.values, chunks=(8, 5, 5)), + dims=mfd.dims, coords=mfd.coords, + ) + # Should not raise at validation time: dask value checks are + # deferred (laziness preserved), not run eagerly here. + out = flow_accumulation_mfd(dask_bad) + assert isinstance(out.data, dask.Array) diff --git a/xrspatial/hydro/watershed_mfd.py b/xrspatial/hydro/watershed_mfd.py index 6dfbd7838..f5a44136f 100644 --- a/xrspatial/hydro/watershed_mfd.py +++ b/xrspatial/hydro/watershed_mfd.py @@ -29,6 +29,7 @@ class cupy: # type: ignore[no-redef] from xrspatial.hydro._boundary_store import BoundaryStore from xrspatial.utils import ( + _validate_mfd_fractions, _validate_raster, has_cuda_and_cupy, is_cupy_array, @@ -678,6 +679,9 @@ def watershed_mfd(flow_dir_mfd: xr.DataArray, raise ValueError( f"flow_dir_mfd must have shape (8, H, W), got {data.shape}") + _validate_mfd_fractions(data, func_name='watershed_mfd', + name='flow_dir_mfd') + _, H, W = data.shape if isinstance(data, np.ndarray): diff --git a/xrspatial/utils.py b/xrspatial/utils.py index 3da402cf7..745faa359 100644 --- a/xrspatial/utils.py +++ b/xrspatial/utils.py @@ -183,6 +183,87 @@ def _validate_scalar( ) +def _validate_mfd_fractions(data, *, func_name: str, name: str = 'fractions', + atol: float = 1e-6): + """Validate the VALUES of a (8, H, W) MFD fraction grid. + + The public MFD functions document that each cell's 8 fraction + bands lie in ``[0, 1]`` and sum to either 1.0 (flow) or 0.0 + (pit/flat/sink), with all-NaN bands at edges and nodata cells. + This checks those value invariants and raises a clear error when + the input violates them, before any hydrology math runs. + + Three checks per cell: + + * No negative fractions. + * The band sum is 1.0 or 0.0 within *atol*. + * NaN bands are all-or-nothing: either all 8 directions are NaN + (edge/nodata) or none are. A partially-NaN cell is rejected. + + Only numpy and cupy (in-memory) arrays are validated. Dask arrays + are skipped so validation does not force computation; lazy + validation is handled separately. The shape is assumed to already + be ``(8, H, W)`` (callers check that first). + + Parameters + ---------- + data : numpy.ndarray, cupy.ndarray, or dask.array.Array + The fraction grid (``DataArray.data``). + func_name : str + Name of the calling function (for error messages). + name : str + Parameter name (for error messages). + atol : float + Absolute tolerance for the band-sum check. + + Raises + ------ + ValueError + If any cell has a negative fraction, a band sum that is + neither ~1.0 nor ~0.0, or a partial-NaN band pattern. + """ + if is_cupy_array(data): + xp = cupy + elif isinstance(data, np.ndarray): + xp = np + else: + # Dask (numpy- or cupy-backed) or other lazy types: skip value + # validation so we do not trigger computation. + return + + nan_mask = xp.isnan(data) + nan_count = nan_mask.sum(axis=0) + # Partial NaN: some but not all of the 8 bands are NaN. + if bool(((nan_count > 0) & (nan_count < 8)).any()): + raise ValueError( + f"{func_name}(): `{name}` has cells with a partial-NaN band " + f"pattern. Each cell must have all 8 direction bands NaN " + f"(edge/nodata) or none of them NaN." + ) + + finite = data[~nan_mask] + if finite.size: + if bool((finite < 0).any()): + raise ValueError( + f"{func_name}(): `{name}` contains negative flow " + f"fractions. Fractions must be in [0, 1]." + ) + + # Per-cell band sums, treating NaN bands as 0 so all-NaN cells sum + # to 0.0 and pass the sink check. + sums = xp.nansum(data, axis=0) + valid_cell = nan_count == 0 + bad_sum = valid_cell & ~( + (xp.abs(sums - 1.0) <= atol) | (xp.abs(sums) <= atol) + ) + if bool(bad_sum.any()): + raise ValueError( + f"{func_name}(): `{name}` has cells whose flow fractions do " + f"not sum to 1.0 (flow) or 0.0 (pit/flat/sink) within " + f"tolerance {atol}." + ) + + def _boundary_to_dask(boundary, is_cupy=False): """Convert a boundary mode string to the value expected by ``dask.array.map_overlap``'s *boundary* parameter.""" From 2b030e38be8059f827a34418830cc2cc3f6705ca Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Tue, 2 Jun 2026 14:59:12 -0700 Subject: [PATCH 2/2] Address review: drop copy in negative check, dedupe error prefix (#2873) --- xrspatial/utils.py | 29 ++++++++++++++--------------- 1 file changed, 14 insertions(+), 15 deletions(-) diff --git a/xrspatial/utils.py b/xrspatial/utils.py index 745faa359..9ceec07e7 100644 --- a/xrspatial/utils.py +++ b/xrspatial/utils.py @@ -231,23 +231,23 @@ def _validate_mfd_fractions(data, *, func_name: str, name: str = 'fractions', # validation so we do not trigger computation. return - nan_mask = xp.isnan(data) - nan_count = nan_mask.sum(axis=0) + prefix = f"{func_name}(): `{name}`" + + nan_count = xp.isnan(data).sum(axis=0) # Partial NaN: some but not all of the 8 bands are NaN. if bool(((nan_count > 0) & (nan_count < 8)).any()): raise ValueError( - f"{func_name}(): `{name}` has cells with a partial-NaN band " - f"pattern. Each cell must have all 8 direction bands NaN " - f"(edge/nodata) or none of them NaN." + f"{prefix} has cells with a partial-NaN band pattern. Each " + f"cell must have all 8 direction bands NaN (edge/nodata) or " + f"none of them NaN." ) - finite = data[~nan_mask] - if finite.size: - if bool((finite < 0).any()): - raise ValueError( - f"{func_name}(): `{name}` contains negative flow " - f"fractions. Fractions must be in [0, 1]." - ) + # NaN < 0 is False, so NaN cells never trip this (no copy needed). + if bool((data < 0).any()): + raise ValueError( + f"{prefix} contains negative flow fractions. Fractions must " + f"be in [0, 1]." + ) # Per-cell band sums, treating NaN bands as 0 so all-NaN cells sum # to 0.0 and pass the sink check. @@ -258,9 +258,8 @@ def _validate_mfd_fractions(data, *, func_name: str, name: str = 'fractions', ) if bool(bad_sum.any()): raise ValueError( - f"{func_name}(): `{name}` has cells whose flow fractions do " - f"not sum to 1.0 (flow) or 0.0 (pit/flat/sink) within " - f"tolerance {atol}." + f"{prefix} has cells whose flow fractions do not sum to 1.0 " + f"(flow) or 0.0 (pit/flat/sink) within tolerance {atol}." )