Calculate rumple on the voxel grid

docs/usage/forest-structure/rumple.md

@@ -17,39 +17,139 @@ Where:
 A flat canopy has a rumple value of 1.0. More structurally complex or
 corrugated canopies have values greater than 1.0.
 
-PyForestScan computes rumple from a Canopy Height Model (CHM) by treating the
-CHM as a triangulated surface over the raster grid and summing surface area
-over valid 2x2 CHM patches.
+PyForestScan returns one rumple value per XY voxel cell. It first takes the
+maximum height above ground in each column to form a canopy surface on the
+same grid as `assign_voxels`. It then uses the cell and its eight neighbors
+to estimate surface area following [Jenness (2004)](https://www.jennessent.com/downloads/WSB_32_3_Jenness.pdf).
+The eight triangles are clipped to the central cell, and their combined
+area is divided by that cell's ground area, `dx * dy`.
+
+The result has the same XY resolution and grid alignment as canopy cover,
+PAI, and FHD calculated with the same points and voxel resolution. The
+neighboring cells describe the surface across each cell; this is a local
+surface-area ratio, not a single average for the point cloud.
 
 ## Calculating Rumple
 
-To calculate rumple:
+Pass the point array and a `(dx, dy, dz)` voxel resolution. The function
+returns a 2D array and its spatial extent, ready to plot or save as a GeoTIFF:
 
 ```python
-from pyforestscan.handlers import read_lidar
-from pyforestscan.calculate import calculate_chm, calculate_rumple
+from pyforestscan.handlers import read_lidar, create_geotiff
+from pyforestscan.calculate import calculate_rumple
+from pyforestscan.visualize import plot_metric
 
 file_path = "../example_data/20191210_5QKB020880.laz"
 arrays = read_lidar(file_path, "EPSG:32605", hag=True)
 points = arrays[0]
 
-cell_resolution = (5.0, 5.0)
-chm, extent = calculate_chm(points, cell_resolution, interpolation="linear")
-rumple = calculate_rumple(chm, cell_resolution, min_height=2.0)
+voxel_resolution = (5.0, 5.0, 1.0)
+rumple, extent = calculate_rumple(points, voxel_resolution, min_height=2.0)
 
-print(f"Rumple: {rumple:.3f}")
+plot_metric("Rumple Index", rumple, extent, metric_name="Rumple", cmap="viridis")
+create_geotiff(rumple, "rumple.tif", "EPSG:32605", extent)
 ```
 
+## Filling Gaps
+
+Use `interpolation="linear"` to fill missing canopy heights before calculating
+rumple. Only empty cells are filled; measured heights are preserved.
+
+```python
+rumple, extent = calculate_rumple(
+    points, voxel_resolution, min_height=2.0, interpolation="linear",
+)
+```
+
+The default, `interpolation=None`, leaves gaps unfilled. The available methods
+match the CHM options:
+
+- `"linear"`: interpolate heights on triangles between observed cell centers.
+- `"cubic"`: use a smooth cubic surface between observed cell centers.
+- `"nearest"`: copy the closest observed height, using XY distances.
+
+Linear and cubic interpolation leave cells outside the observed cell centers'
+convex hull as NoData. They also leave gaps unfilled if there are too few
+noncollinear samples to form a surface. Nearest can fill outside that hull,
+but does not expand the raster extent. A complete 3x3 neighborhood is still
+needed for rumple, so the outermost row and column remain NoData.
+
+Interpolation happens before `min_height` is applied. Observed low or ground
+cells remain subject to the height mask; they are not replaced with higher
+canopy values. Interpolated heights below the threshold are masked too.
+
+This follows the same sequence as lidR: its
+[`p2r(na.fill = tin())`](https://search.r-project.org/CRAN/refmans/lidR/html/dsm_point2raster.html)
+option fills the canopy height model before `rumple_index()` is calculated.
+Here interpolation uses the gridded canopy maxima. For a numerical comparison,
+use the same filled height grid in both implementations. Filling gaps changes
+the estimated canopy surface and can affect roughness, so use consistent
+resolution and interpolation settings when comparing results.
+
+## Tiled GeoTIFF Output
+
+For large EPT point clouds, use `metric="rumple"`:
+
+```python
+from pyforestscan.process import process_with_tiles
+
+process_with_tiles(
+    ept_file="../example_data/ept/ept.json",
+    tile_size=(1000, 1000),
+    output_path="rumple_tiles",
+    metric="rumple",
+    voxel_size=(5.0, 5.0, 1.0),
+    srs="EPSG:32605",
+    hag=True,
+    rumple_min_height=2.0,
+    interpolation="linear",  # Optional: fill canopy gaps before calculating rumple.
+)
+```
+
+This writes `tile_<i>_<j>_rumple.tif` files. Tile dimensions must be multiples
+of the XY voxel sizes. Output bounds expand to whole voxel cells, and the
+last tile can be smaller than `tile_size` without changing pixel size.
+Each tile reads at least one extra cell on every available side before
+calculating rumple, then crops to its output grid. This also applies when
+`buffer_size=0`, so adjacent tiles retain the neighboring canopy data needed
+at their shared edges. `skip_existing` and `tile_indices` work as for the
+other tiled metrics.
+
+When interpolation is enabled, it uses the points in each buffered tile.
+Increase `buffer_size` enough to include observed canopy around gaps near tile
+edges. The one-cell minimum buffer supports the surface-area calculation, but
+may not be enough for interpolation across larger gaps. Filled values can
+differ from a whole-cloud calculation because fewer samples are available.
+The `interpolation` option applies to CHM and rumple only.
+
 ## Notes
 
-- `calculate_rumple` returns a single scalar value, not a raster.
-- `min_height` can be used to exclude low vegetation before calculating
-  canopy surface complexity.
-- Interpolating the CHM before calculating rumple may fill gaps, but it can
-  also smooth the canopy surface and reduce rumple slightly.
+- The array is shaped `(X, Y)`, with Y ordered north to south, matching
+  `assign_voxels` and `create_geotiff`.
+- `dx` and `dy` determine both the surface sampling and output pixel size.
+  Changing them changes the scale of canopy roughness being measured. Use
+  the same resolution when comparing sites.
+- `dz` is accepted as part of the shared voxel-resolution tuple. Actual
+  maximum heights are used without rounding them to vertical bins, so
+  changing `dz` alone does not change rumple.
+- XY coordinates and heights must use the same linear units, such as meters
+  in a projected CRS.
+- A complete 3x3 canopy neighborhood is required. The outermost row/column,
+  unfilled cells, and cells adjacent to missing or height-masked canopy
+  remain NaN, written as NoData in the GeoTIFF. Grids smaller than three
+  cells in either direction contain only NoData.
+- `min_height` masks canopy maxima strictly below the threshold; it does
+  not change the grid extent. Interpolation is optional and defaults to None.
+- This replaces the earlier scalar API. Update
+  `calculate_rumple(chm, cell_resolution)` calls to
+  `rumple, extent = calculate_rumple(points, voxel_resolution)`.
 
 ## References
 
+Jenness, Jeff S. 2004. "Calculating landscape surface area from digital
+elevation models." Wildlife Society Bulletin 32 (3): 829-839.
+<https://www.jennessent.com/downloads/WSB_32_3_Jenness.pdf>.
+
 McElhinny, Chris, Phillip Gibbons, Cris Brack, and Juergen Bauhus. 2005.
 "Forest and woodland stand structural complexity: Its definition and
 measurement." Forest Ecology and Management 218 (1-3): 1-24.

pyforestscan/calculate.py

@@ -1,6 +1,7 @@
 import numpy as np
 
 from scipy.interpolate import griddata
+from scipy.spatial import QhullError
 from scipy.stats import entropy
 from scipy import ndimage
 from typing import List, Tuple, Optional
@@ -584,74 +585,112 @@ def calculate_chm(arr, voxel_resolution, interpolation="linear",
     return chm, extent
 
 
-def calculate_rumple(chm: np.ndarray,
-                     cell_resolution: Tuple[float, float],
-                     min_height: float | None = None) -> float:
+def calculate_rumple(arr: np.ndarray,
+                     voxel_resolution: Tuple[float, float, float],
+                     min_height: float | None = None,
+                     interpolation: str | None = None) -> Tuple[np.ndarray, List]:
     """
-    Calculate the canopy rumple index from a Canopy Height Model (CHM).
+    Calculate a rumple raster on the same XY grid as ``assign_voxels``.
 
-    Rumple is defined here as the ratio of canopy surface area to planar
-    ground area. The CHM is treated as a triangulated surface over the raster
-    grid, and the surface area is summed over valid 2x2 CHM patches.
+    The highest HeightAboveGround value in each voxel column defines the
+    canopy surface. Following Jenness (2004), eight triangles connect each
+    cell center to its neighbors. The portions inside the central cell are
+    summed and divided by its planar area (dx * dy). Flat surfaces have a
+    rumple of 1; sloped or rough surfaces have values greater than 1.
 
     Args:
-        chm (np.ndarray): 2D array of canopy heights.
-        cell_resolution (tuple[float, float]): CHM cell size as (dx, dy).
-        min_height (float | None, optional): If provided, CHM cells below this
-            height threshold are masked before computing the rumple index.
-            Defaults to None.
+        arr (np.ndarray): Structured point array with X, Y, and
+            HeightAboveGround fields. Nonfinite points and points below
+            ground are ignored.
+        voxel_resolution (tuple[float, float, float]): Positive, finite
+            (dx, dy, dz) sizes in the same units as the point coordinates.
+            XY sets the raster resolution; dz does not quantize the canopy
+            heights or affect rumple.
+        min_height (float | None, optional): Mask canopy cells below this
+            height after interpolation and before calculating rumple.
+            Observed cells below this threshold are not filled. Defaults to None.
+        interpolation (str | None, optional): Fill missing canopy heights
+            using "linear", "cubic", or "nearest" interpolation before
+            calculating surface area. Observed heights are preserved.
+            Linear and cubic leave gaps outside the observed cell centers'
+            convex hull, or without enough noncollinear samples, as NaN.
+            Nearest can also fill outside that hull within the raster extent.
+            Defaults to None (no filling).
 
     Returns:
-        float: Rumple index (>= 1 for valid surfaces) or NaN if no valid 2x2
-            surface patches remain after masking.
+        tuple[np.ndarray, list]: Rumple array shaped (X, Y), with Y ordered
+            north to south, and extent [x_min, x_max, y_min, y_max]. A complete
+            3x3 canopy neighborhood is required. Outer edges and neighborhoods
+            with remaining missing or masked canopy are NaN, even when
+            interpolation is enabled.
 
     Raises:
-        ValueError: If the CHM is not 2D, if cell_resolution is invalid, or if
-            dx/dy are not positive.
+        ValueError: If the resolution, height threshold, or interpolation
+            method is invalid, or no finite points at or above ground remain.
+        KeyError: If a required point dimension is missing.
     """
-    chm = np.asarray(chm, dtype=float)
-    if chm.ndim != 2:
-        raise ValueError(f"chm must be a 2D array (got shape {chm.shape})")
+    resolution = np.asarray(voxel_resolution, dtype=float)
+    if resolution.shape != (3,) or not np.all(np.isfinite(resolution) & (resolution > 0)):
+        raise ValueError("voxel_resolution must contain three positive, finite sizes (dx, dy, dz)")
+    if min_height is not None and not np.isfinite(min_height):
+        raise ValueError("min_height must be finite or None")
+    if interpolation not in (None, "linear", "cubic", "nearest"):
+        raise ValueError("interpolation must be None, 'linear', 'cubic', or 'nearest'")
+
+    arr = np.asarray(arr)
+    required = ('X', 'Y', 'HeightAboveGround')
+    if arr.dtype.names is None or not all(name in arr.dtype.names for name in required):
+        raise KeyError("Input array must include X, Y, and HeightAboveGround fields")
+    if arr.ndim != 1:
+        raise ValueError("Input point array must be one-dimensional")
+    valid = np.isfinite(arr['X']) & np.isfinite(arr['Y']) & np.isfinite(arr['HeightAboveGround'])
+    points = arr[valid & (arr['HeightAboveGround'] >= 0)]
+    if points.size == 0:
+        raise ValueError("No finite points at or above ground are available")
+
+    chm, extent = calculate_voxel_stat(points, resolution, 'HeightAboveGround', 'max')
+    rumple = np.full(chm.shape, np.nan)
+    if min(chm.shape) < 3:
+        return rumple, extent
 
-    if len(cell_resolution) != 2:
-        raise ValueError("cell_resolution must be a (dx, dy) tuple")
-
-    dx, dy = map(float, cell_resolution)
-    if dx <= 0 or dy <= 0:
-        raise ValueError("cell_resolution components must be > 0")
+    if interpolation is not None:
+        valid_cells = np.isfinite(chm)
+        if not valid_cells.all():
+            # Use physical distances so rectangular pixels interpolate correctly.
+            # Fill heights before the threshold mask; masked ground is not a gap.
+            known = np.argwhere(valid_cells) * resolution[:2]
+            missing = np.argwhere(~valid_cells) * resolution[:2]
+            try:
+                chm[~valid_cells] = griddata(
+                    known, chm[valid_cells], missing, method=interpolation,
+                )
+            except QhullError:
+                # Linear/cubic need a 2D triangulation. Sparse or collinear
+                # samples cannot support one, so their gaps remain NoData.
+                pass
 
     if min_height is not None:
         chm = np.where(chm >= float(min_height), chm, np.nan)
 
-    z00 = chm[:-1, :-1]
-    z10 = chm[1:, :-1]
-    z01 = chm[:-1, 1:]
-    z11 = chm[1:, 1:]
-
-    valid = (
-        np.isfinite(z00) &
-        np.isfinite(z10) &
-        np.isfinite(z01) &
-        np.isfinite(z11)
-    )
-    if not np.any(valid):
-        return np.nan
-
-    # Approximate the CHM as a triangular mesh over each 2x2 raster patch.
-    tri1 = 0.5 * np.sqrt(
-        (dy * (z10 - z00)) ** 2 +
-        (dx * (z01 - z00)) ** 2 +
-        (dx * dy) ** 2
-    )
-    tri2 = 0.5 * np.sqrt(
-        (dy * (z01 - z11)) ** 2 +
-        (dx * (z11 - z10)) ** 2 +
-        (dx * dy) ** 2
-    )
-
-    surface_area = np.sum((tri1 + tri2)[valid], dtype=float)
-    planar_area = float(np.count_nonzero(valid)) * dx * dy
-    return surface_area / planar_area
+    dx, dy = resolution[:2]
+    center = chm[1:-1, 1:-1]
+    nx, ny = center.shape
+    offsets = [(1, 0), (1, 1), (0, 1), (-1, 1),
+               (-1, 0), (-1, -1), (0, -1), (1, -1)]
+    surface_ratio = np.zeros(center.shape)
+    for (ax, ay), (bx, by) in zip(offsets, offsets[1:] + offsets[:1]):
+        za = chm[1 + ax:1 + ax + nx, 1 + ay:1 + ay + ny] - center
+        zb = chm[1 + bx:1 + bx + nx, 1 + by:1 + by + ny] - center
+
+        # Cross-product area, normalized by dx*dy. Halving the two
+        # center-to-neighbor vectors clips each triangle to the cell:
+        # area = |cross| / 8. All eight projected areas sum to dx*dy.
+        slope_x = (ay * zb - by * za) / dx
+        slope_y = (bx * za - ax * zb) / dy
+        surface_ratio += np.hypot(np.hypot(slope_x, slope_y), 1.0) / 8.0
+
+    rumple[1:-1, 1:-1] = surface_ratio
+    return rumple, extent
 
 
 def _calc_valid_region_mask(arr):