Fixes #804: add min_observable_height() to experimental (#875) (#886)

New function that computes the minimum observer height needed to see
each grid cell from a given location.  Uses binary search over
viewshed evaluations (log₂(max_height/precision) calls) instead of
a linear sweep.

Author: brendancol

README.md

@@ -226,6 +226,7 @@ In the GIS world, rasters are used for representing continuous phenomena (e.g. e
 | [Slope](xrspatial/slope.py) | Computes terrain gradient steepness at each cell in degrees | ✅️  | ✅️  | ✅️ | ✅️  |
 | [Terrain Generation](xrspatial/terrain.py) | Generates synthetic terrain elevation using fractal noise | ✅️ | ✅️ | ✅️ | |
 | [Viewshed](xrspatial/viewshed.py) | Determines visible cells from a given observer point on terrain | ✅️ |  | | |
+| [Min Observable Height](xrspatial/experimental/min_observable_height.py) | Finds the minimum observer height needed to see each cell *(experimental)* | ✅️ | | | |
 | [Perlin Noise](xrspatial/perlin.py) | Generates smooth continuous random noise for procedural textures | ✅️ | ✅️ | ✅️ | |
 | [Bump Mapping](xrspatial/bump.py) | Adds randomized bump features to simulate natural terrain variation | ✅️ | | | |
 

xrspatial/accessor.py

@@ -43,6 +43,10 @@ def viewshed(self, x, y, **kwargs):
         from .viewshed import viewshed
         return viewshed(self._obj, x, y, **kwargs)
 
+    def min_observable_height(self, x, y, **kwargs):
+        from .experimental.min_observable_height import min_observable_height
+        return min_observable_height(self._obj, x, y, **kwargs)
+
     # ---- Classification ----
 
     def natural_breaks(self, **kwargs):

xrspatial/experimental/__init__.py

@@ -1 +1,2 @@
+from .min_observable_height import min_observable_height  # noqa
 from .polygonize import polygonize  # noqa

xrspatial/experimental/min_observable_height.py

@@ -0,0 +1,155 @@
+"""Compute the minimum observer height needed to see each grid cell.
+
+For a given observer location on a terrain raster, this module determines
+the minimum height above the terrain that an observer must be raised to
+in order to see each cell.  The result is a DataArray of the same shape
+as the input, where each cell contains that minimum height (or ``NaN``
+for cells that remain invisible even at *max_height*).
+
+The algorithm exploits the fact that visibility is **monotonic** in
+observer height: if a cell is visible at height *h*, it is also visible
+at every height greater than *h*.  This allows a per-cell binary search
+over observer height, requiring only ``ceil(log2(max_height / precision))``
+viewshed evaluations instead of a linear sweep.
+"""
+
+from math import ceil, log2
+from typing import Union
+
+import numpy as np
+import xarray
+
+from ..viewshed import viewshed, INVISIBLE
+
+
+def min_observable_height(
+    raster: xarray.DataArray,
+    x: Union[int, float],
+    y: Union[int, float],
+    max_height: float = 100.0,
+    target_elev: float = 0.0,
+    precision: float = 1.0,
+) -> xarray.DataArray:
+    """Return the minimum observer height to see each cell.
+
+    For every cell in *raster*, compute the smallest observer elevation
+    (above the terrain at the observer location) at which that cell
+    becomes visible.  Cells that are never visible up to *max_height*
+    are set to ``NaN``.
+
+    Parameters
+    ----------
+    raster : xarray.DataArray
+        2-D elevation raster with ``y`` and ``x`` coordinates.
+    x : int or float
+        x-coordinate (in data space) of the observer location.
+    y : int or float
+        y-coordinate (in data space) of the observer location.
+    max_height : float, default 100.0
+        Maximum observer height to test.  Cells invisible at this
+        height receive ``NaN`` in the output.
+    target_elev : float, default 0.0
+        Height offset added to target cells when testing visibility,
+        passed through to :func:`~xrspatial.viewshed.viewshed`.
+    precision : float, default 1.0
+        Height resolution of the binary search.  The result for each
+        cell is accurate to within *precision* units.
+
+    Returns
+    -------
+    xarray.DataArray
+        Same shape and coordinates as *raster*.  Each cell contains the
+        minimum observer height (in surface units) at which it becomes
+        visible, or ``NaN`` if the cell is not visible at *max_height*.
+        The cell containing the observer is always ``0.0``.
+
+    Raises
+    ------
+    ValueError
+        If *max_height* <= 0, *precision* <= 0, or *precision* >
+        *max_height*.
+
+    Notes
+    -----
+    The function calls :func:`~xrspatial.viewshed.viewshed` internally,
+    so the same backend (numpy / cupy) is used.  Dask-backed rasters
+    are not supported because viewshed itself does not support them.
+
+    The number of viewshed evaluations is
+    ``ceil(log2(max_height / precision)) + 1``, which is typically far
+    fewer than a linear sweep at the same resolution.
+
+    Examples
+    --------
+    .. sourcecode:: python
+
+        >>> import numpy as np
+        >>> import xarray as xr
+        >>> from xrspatial.experimental import min_observable_height
+        >>> elev = np.zeros((5, 5), dtype=float)
+        >>> elev[2, 4] = 10.0  # a tall obstacle
+        >>> raster = xr.DataArray(elev, dims=['y', 'x'])
+        >>> h, w = elev.shape
+        >>> raster['y'] = np.arange(h, dtype=float)
+        >>> raster['x'] = np.arange(w, dtype=float)
+        >>> result = min_observable_height(raster, x=0.0, y=2.0,
+        ...                               max_height=20, precision=0.5)
+    """
+    # ---- input validation ------------------------------------------------
+    if max_height <= 0:
+        raise ValueError(f"max_height must be positive, got {max_height}")
+    if precision <= 0:
+        raise ValueError(f"precision must be positive, got {precision}")
+    if precision > max_height:
+        raise ValueError(
+            f"precision ({precision}) must not exceed max_height ({max_height})"
+        )
+
+    # ---- determine candidate heights via binary search -------------------
+    #
+    # We evaluate viewshed at O(log2(max_height / precision)) heights.
+    # For each cell we then find the smallest height where it is visible.
+    #
+    n_steps = int(ceil(log2(max_height / precision))) + 1
+    # Build sorted list of heights from 0 to max_height (inclusive).
+    # Using linspace guarantees both endpoints are included.
+    if n_steps < 2:
+        n_steps = 2
+    heights = np.linspace(0.0, max_height, n_steps)
+
+    # ---- collect visibility at each height level -------------------------
+    # visibility_stack[k] is True where the cell is visible at heights[k].
+    visibility_stack = np.empty(
+        (n_steps, *raster.shape), dtype=np.bool_
+    )
+
+    for k, h in enumerate(heights):
+        vs = viewshed(raster, x=x, y=y, observer_elev=float(h),
+                       target_elev=target_elev)
+        # viewshed returns INVISIBLE (-1) for invisible, positive angle
+        # for visible.  The observer cell itself gets 180.
+        vs_data = vs.values if isinstance(vs.data, np.ndarray) else vs.data.get()
+        visibility_stack[k] = vs_data != INVISIBLE
+
+    # ---- per-cell binary search on the precomputed stack -----------------
+    # Because visibility is monotonic in height, we can find the first
+    # height index where each cell becomes visible.
+    #
+    # Strategy: scan from lowest to highest height.  The minimum height
+    # for a cell is the first index k where visibility_stack[k] is True.
+    #
+    # np.argmax on a boolean axis returns the *first* True index.
+    # If a cell is never visible, all entries are False and argmax
+    # returns 0 — we need to distinguish that case.
+
+    ever_visible = visibility_stack.any(axis=0)
+    first_visible_idx = visibility_stack.argmax(axis=0)
+
+    result = np.where(ever_visible, heights[first_visible_idx], np.nan)
+
+    return xarray.DataArray(
+        result,
+        coords=raster.coords,
+        dims=raster.dims,
+        attrs=raster.attrs,
+    )

xrspatial/tests/test_min_observable_height.py

@@ -0,0 +1,198 @@
+import numpy as np
+import pytest
+import xarray as xa
+
+from xrspatial.experimental import min_observable_height
+
+
+def _make_raster(data, xs=None, ys=None):
+    """Build a DataArray with y/x coords from a 2-D numpy array."""
+    ny, nx = data.shape
+    if xs is None:
+        xs = np.arange(nx, dtype=float)
+    if ys is None:
+        ys = np.arange(ny, dtype=float)
+    return xa.DataArray(
+        data.astype(np.float64),
+        coords=dict(x=xs, y=ys),
+        dims=["y", "x"],
+    )
+
+
+# ---- input validation ---------------------------------------------------
+
+def test_max_height_positive():
+    r = _make_raster(np.zeros((3, 3)))
+    with pytest.raises(ValueError, match="max_height must be positive"):
+        min_observable_height(r, x=1.0, y=1.0, max_height=0)
+
+
+def test_precision_positive():
+    r = _make_raster(np.zeros((3, 3)))
+    with pytest.raises(ValueError, match="precision must be positive"):
+        min_observable_height(r, x=1.0, y=1.0, precision=0)
+
+
+def test_precision_exceeds_max_height():
+    r = _make_raster(np.zeros((3, 3)))
+    with pytest.raises(ValueError, match="precision.*must not exceed"):
+        min_observable_height(r, x=1.0, y=1.0, max_height=5, precision=10)
+
+
+# ---- flat terrain --------------------------------------------------------
+
+def test_flat_terrain_all_visible_at_zero():
+    """On flat terrain with any positive observer_elev, all cells are
+    visible.  min_observable_height should be 0 everywhere."""
+    data = np.full((5, 5), 10.0)
+    r = _make_raster(data)
+    result = min_observable_height(r, x=2.0, y=2.0,
+                                   max_height=20.0, precision=1.0)
+    assert result.shape == r.shape
+    # Every cell should be visible at height 0 (observer at terrain level
+    # can see all cells at same elevation).
+    assert np.all(result.values == 0.0)
+
+
+# ---- obstacle blocking ---------------------------------------------------
+
+def test_wall_requires_height():
+    """A wall between the observer and far cells forces higher observer."""
+    #  y=0: 0  0  0  0  0
+    #  y=1: 0  0 10  0  0    <- wall at (y=1, x=2)
+    #  y=2: 0  0  0  0  0
+    #  Observer at (x=0, y=1).
+    #
+    #  To see cell (x=3, y=1) over the wall, the line of sight from
+    #  (0, h) to (3, 0) must clear (2, 10).  LOS height at x=2 is h/3,
+    #  so h >= 30 is required.
+    data = np.zeros((3, 5))
+    data[1, 2] = 10.0  # wall
+    r = _make_raster(data)
+
+    result = min_observable_height(r, x=0.0, y=1.0,
+                                   max_height=50.0, precision=1.0)
+
+    # Observer cell itself should be 0.
+    assert result.sel(x=0.0, y=1.0).item() == 0.0
+
+    # Cells on same side as observer (x=1) should be visible at 0.
+    assert result.sel(x=1.0, y=1.0).item() == 0.0
+
+    # Cell behind the wall (x=3, y=1) needs a significant height.
+    behind_wall = result.sel(x=3.0, y=1.0).item()
+    assert behind_wall > 0.0, "cells behind wall should need height > 0"
+
+    # The wall cell itself is visible from height 0 (it's tall, angles
+    # up to observer) or at least from some modest height.
+    wall_h = result.sel(x=2.0, y=1.0).item()
+    assert wall_h < behind_wall, "wall cell should be easier to see than cells behind it"
+
+
+def test_tall_peak_blocks_cells_behind():
+    """A peak should block cells directly behind it, requiring more height."""
+    # Observer at (x=0, y=1), peak at (x=2, y=1), target at (x=4, y=1).
+    # Need at least 2 rows so viewshed can compute ns_res.
+    data = np.zeros((3, 5))
+    data[1, 2] = 8.0  # peak
+    r = _make_raster(data)
+
+    result = min_observable_height(r, x=0.0, y=1.0,
+                                   max_height=40.0, precision=0.5)
+
+    h_behind = result.values[1, 4]
+    h_front = result.values[1, 1]
+    # cell in front of peak should be visible from lower height
+    assert h_front < h_behind
+
+
+# ---- NaN for unreachable cells -------------------------------------------
+
+def test_unreachable_cells_are_nan():
+    """If a cell cannot be seen even at max_height, result should be NaN."""
+    # Very tall wall close to observer, low max_height.
+    data = np.zeros((3, 3))
+    data[1, 1] = 1000.0  # enormous wall next to observer
+    r = _make_raster(data)
+
+    result = min_observable_height(r, x=0.0, y=1.0,
+                                   max_height=5.0, precision=1.0)
+
+    # Some cells behind the 1000-unit wall should be NaN at max_height=5.
+    # The far corner (x=2, y=0 or y=2) is behind the wall.
+    far_val = result.sel(x=2.0, y=0.0).item()
+    assert np.isnan(far_val), "cell behind huge wall should be NaN at low max_height"
+
+
+# ---- output shape and coords ---------------------------------------------
+
+def test_output_shape_coords_match_input():
+    data = np.random.default_rng(42).uniform(0, 10, (6, 8))
+    xs = np.linspace(0, 7, 8)
+    ys = np.linspace(0, 5, 6)
+    r = _make_raster(data, xs=xs, ys=ys)
+
+    result = min_observable_height(r, x=3.0, y=2.0,
+                                   max_height=20, precision=2.0)
+    assert result.shape == r.shape
+    np.testing.assert_array_equal(result.coords['x'].values, xs)
+    np.testing.assert_array_equal(result.coords['y'].values, ys)
+    assert result.dims == r.dims
+
+
+# ---- monotonicity property -----------------------------------------------
+
+def test_monotonicity():
+    """Visibility is monotonic: if visible at h, visible at all h' > h.
+
+    This indirectly tests that the binary search is correct: for each
+    cell, all heights >= min_observable_height should show it as visible.
+    """
+    data = np.zeros((5, 5))
+    data[2, 3] = 8.0
+    r = _make_raster(data)
+
+    result = min_observable_height(r, x=0.0, y=2.0,
+                                   max_height=20.0, precision=1.0)
+
+    # Pick a cell that requires non-zero height.
+    h_needed = result.sel(x=4.0, y=2.0).item()
+    if not np.isnan(h_needed) and h_needed > 0:
+        from xrspatial import viewshed
+        # Should be invisible slightly below threshold.
+        vs_below = viewshed(r, x=0.0, y=2.0, observer_elev=max(h_needed - 1.0, 0))
+        vs_above = viewshed(r, x=0.0, y=2.0, observer_elev=h_needed + 1.0)
+        # At or above threshold the cell should be visible.
+        assert vs_above.sel(x=4.0, y=2.0).item() != -1
+
+
+# ---- precision affects granularity ---------------------------------------
+
+def test_finer_precision():
+    """With finer precision, results should be at least as good (<=) as
+    coarser precision."""
+    data = np.zeros((3, 5))
+    data[1, 2] = 10.0
+    r = _make_raster(data)
+
+    coarse = min_observable_height(r, x=0.0, y=1.0,
+                                   max_height=20, precision=5.0)
+    fine = min_observable_height(r, x=0.0, y=1.0,
+                                 max_height=20, precision=0.5)
+
+    # Fine result should be <= coarse result for all finite cells
+    # (finer search can find an equal or lower threshold).
+    mask = np.isfinite(fine.values) & np.isfinite(coarse.values)
+    assert np.all(fine.values[mask] <= coarse.values[mask] + 1e-10)
+
+
+# ---- accessor integration ------------------------------------------------
+
+def test_accessor():
+    """min_observable_height is available via .xrs accessor."""
+    import xrspatial  # noqa: F401 – registers accessor
+    data = np.zeros((3, 3))
+    r = _make_raster(data)
+    result = r.xrs.min_observable_height(x=1.0, y=1.0,
+                                          max_height=10, precision=2)
+    assert result.shape == r.shape