From 9bb9853d7f61c9c8f49aa8a639c5c0c4e78446d4 Mon Sep 17 00:00:00 2001 From: iosefa Date: Mon, 18 May 2026 22:56:04 -1000 Subject: [PATCH 1/4] Add rumple metric and tiled processing controls --- docs/usage/forest-structure/intro.md | 1 + docs/usage/forest-structure/rumple.md | 62 ++++++++++++++++++++++++ pyforestscan/__init__.py | 2 + pyforestscan/calculate.py | 70 +++++++++++++++++++++++++++ pyforestscan/process.py | 44 +++++++++++++++-- tests/test_calculate.py | 69 +++++++++++++++++++++++++- 6 files changed, 243 insertions(+), 5 deletions(-) create mode 100644 docs/usage/forest-structure/rumple.md diff --git a/docs/usage/forest-structure/intro.md b/docs/usage/forest-structure/intro.md index 4c19bf7..6774074 100644 --- a/docs/usage/forest-structure/intro.md +++ b/docs/usage/forest-structure/intro.md @@ -10,6 +10,7 @@ Most calculations for the metrics follows the method outlined in Kamoske et al. ## Metrics * [Canopy Height Models (CHM)](chm.md) +* [Rumple Index](rumple.md) * [Plant Area Density (PAD)](pad.md) * [Plant Area Index (PAI)](pai.md) * [Foliage Height Diversity (FHD)](fhd.md) diff --git a/docs/usage/forest-structure/rumple.md b/docs/usage/forest-structure/rumple.md new file mode 100644 index 0000000..ee96078 --- /dev/null +++ b/docs/usage/forest-structure/rumple.md @@ -0,0 +1,62 @@ +# Rumple Index + +## Theory + +Rumple is a measure of canopy surface complexity. It is defined as the ratio +of canopy surface area to projected ground area: + +$$ +\text{Rumple} = \frac{A_{\text{surface}}}{A_{\text{planar}}} +$$ + +Where: + +- \( A_{\text{surface}} \) is the area of the canopy surface. +- \( A_{\text{planar}} \) is the projected ground area beneath that surface. + +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. + +## Calculating Rumple + +To calculate rumple: + +```python +from pyforestscan.handlers import read_lidar +from pyforestscan.calculate import calculate_chm, calculate_rumple + +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) + +print(f"Rumple: {rumple:.3f}") +``` + +## 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. + +## References + +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. +. + +Kane, Van R., Jonathan D. Bakker, Robert J. McGaughey, James A. Lutz, +Rolf F. Gersonde, and Jerry F. Franklin. 2010. "Examining conifer canopy +structural complexity across forest ages and elevations with LiDAR data." +Canadian Journal of Forest Research 40 (4): 774-787. +. diff --git a/pyforestscan/__init__.py b/pyforestscan/__init__.py index e168024..9635961 100644 --- a/pyforestscan/__init__.py +++ b/pyforestscan/__init__.py @@ -4,6 +4,7 @@ calculate_pai, calculate_fhd, calculate_chm, + calculate_rumple, calculate_point_density, calculate_voxel_stat, generate_dtm, @@ -16,6 +17,7 @@ "calculate_pai", "calculate_fhd", "calculate_chm", + "calculate_rumple", "calculate_point_density", "calculate_voxel_stat", "generate_dtm", diff --git a/pyforestscan/calculate.py b/pyforestscan/calculate.py index 635a0c5..bc05c4e 100644 --- a/pyforestscan/calculate.py +++ b/pyforestscan/calculate.py @@ -584,6 +584,76 @@ 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: + """ + Calculate the canopy rumple index from a Canopy Height Model (CHM). + + 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. + + 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. + + Returns: + float: Rumple index (>= 1 for valid surfaces) or NaN if no valid 2x2 + surface patches remain after masking. + + Raises: + ValueError: If the CHM is not 2D, if cell_resolution is invalid, or if + dx/dy are not positive. + """ + chm = np.asarray(chm, dtype=float) + if chm.ndim != 2: + raise ValueError(f"chm must be a 2D array (got shape {chm.shape})") + + 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 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 + + def _calc_valid_region_mask(arr): """Create valid region mask using morphological operations.""" kernel = ndimage.generate_binary_structure(2, 1) diff --git a/pyforestscan/process.py b/pyforestscan/process.py index 5ae2a66..d1f8917 100644 --- a/pyforestscan/process.py +++ b/pyforestscan/process.py @@ -8,7 +8,7 @@ from pyforestscan.calculate import calculate_fhd, calculate_pad, calculate_pai, assign_voxels, calculate_chm, calculate_canopy_cover from pyforestscan.filters import remove_outliers_and_clean, downsample_poisson, downsample_voxel from pyforestscan.handlers import create_geotiff -from pyforestscan.pipeline import _hag_raster, _hag_delaunay +from pyforestscan.pipeline import _filter_expression, _filter_statistical_outlier, _hag_raster, _hag_delaunay from pyforestscan.utils import get_bounds_from_ept, get_srs_from_ept @@ -51,13 +51,16 @@ def _crop_dtm(dtm_path, tile_min_x, tile_min_y, tile_max_x, tile_max_y): def process_with_tiles(ept_file, tile_size, output_path, metric, voxel_size, voxel_height=1, buffer_size=0.1, srs=None, hag=False, hag_dtm=False, dtm=None, bounds=None, interpolation=None, remove_outliers=False, + outlier_mean_k: int = 8, outlier_multiplier: float = 3.0, cover_min_height: float = 2.0, cover_k: float = 0.5, pai_min_height: float = 1.0, fhd_min_height: float = 0.0, skip_existing: bool = False, verbose: bool = False, thin_radius: float | None = None, voxelgrid_cell: float | None = None, - voxelgrid_mode: str = "first") -> None: + voxelgrid_mode: str = "first", + tile_indices: set[tuple[int, int]] | None = None, + outliers_before_hag: bool = False) -> None: """ Process a large EPT point cloud by tiling, compute CHM or other metrics for each tile, and write the results to the specified output directory. @@ -79,6 +82,10 @@ def process_with_tiles(ept_file, tile_size, output_path, metric, voxel_size, If None, tiling is done over the entire dataset. interpolation (str or None, optional): Interpolation method for CHM calculation ("linear", "cubic", "nearest", or None). remove_outliers (bool, optional): Whether to remove statistical outliers before calculating metrics. Defaults to False. + outlier_mean_k (int, optional): Number of nearest neighbors used by the statistical outlier filter. + Used only when remove_outliers is True. Defaults to 8. + outlier_multiplier (float, optional): Standard deviation multiplier used by the statistical outlier filter. + Used only when remove_outliers is True. Defaults to 3.0. cover_min_height (float, optional): Height threshold (in meters) for canopy cover (used when metric == "cover"). Defaults to 2.0. cover_k (float, optional): Beer–Lambert extinction coefficient for canopy cover. Defaults to 0.5. pai_min_height (float, optional): Minimum height (m) to integrate PAI. Defaults to 1.0. @@ -91,6 +98,11 @@ def process_with_tiles(ept_file, tile_size, output_path, metric, voxel_size, using cell edge length of ``voxelgrid_cell``. Defaults to None. voxelgrid_mode (str, optional): Representative selection for voxel-grid downsampling. Common values are "first" (keep first point unchanged) or "center" (snap kept point to voxel center). Defaults to "first". + tile_indices (set[tuple[int, int]] or None, optional): If provided, process only these + zero-based tile indices while preserving the full tiling grid and buffer behavior. + outliers_before_hag (bool, optional): If True with remove_outliers=True, apply PDAL's + statistical outlier filter and remove classification 7 points before HAG. This can + avoid Delaunay failures caused by bad points. Defaults to False. Returns: None @@ -124,6 +136,10 @@ def process_with_tiles(ept_file, tile_size, output_path, metric, voxel_size, with tqdm(total=total_tiles, desc="Processing tiles") as pbar: for i in range(num_tiles_x): for j in range(num_tiles_y): + if tile_indices is not None and (i, j) not in tile_indices: + pbar.update(1) + continue + # Apply buffer+crop for CHM and for PAI/COVER to avoid seam artifacts. if metric in ["chm", "pai", "cover"]: current_buffer_size = buffer_size @@ -163,6 +179,15 @@ def process_with_tiles(ept_file, tile_size, output_path, metric, voxel_size, tile_bounds = ([tile_min_x, tile_max_x], [tile_min_y, tile_max_y]) tile_pipeline_stages = [] + if remove_outliers and outliers_before_hag: + tile_pipeline_stages.extend([ + _filter_statistical_outlier( + mean_k=outlier_mean_k, + multiplier=outlier_multiplier, + ), + _filter_expression("!(Classification == 7)"), + ]) + if hag: tile_pipeline_stages.append(_hag_delaunay()) elif hag_dtm: @@ -194,8 +219,19 @@ def process_with_tiles(ept_file, tile_size, output_path, metric, voxel_size, pbar.update(1) continue - if remove_outliers: - tile_points = remove_outliers_and_clean(arrays)[0] + if remove_outliers and not outliers_before_hag: + cleaned_arrays = remove_outliers_and_clean( + arrays, + mean_k=outlier_mean_k, + multiplier=outlier_multiplier, + remove=True, + ) + tile_points = cleaned_arrays[0] + if tile_points.size == 0: + if verbose: + print(f"Warning: Tile ({i}, {j}) empty after outlier removal. Skipping.") + pbar.update(1) + continue else: tile_points = arrays[0] diff --git a/tests/test_calculate.py b/tests/test_calculate.py index 1dfab2e..a13901c 100644 --- a/tests/test_calculate.py +++ b/tests/test_calculate.py @@ -8,7 +8,8 @@ calculate_pai, calculate_fhd, calculate_chm, - calculate_canopy_cover + calculate_canopy_cover, + calculate_rumple ) import math @@ -628,6 +629,72 @@ def test_calculate_chm_large_heights(): assert np.all(chm >= 1000) +# ---------------------------- +# Tests for calculate_rumple +# ---------------------------- + +def test_calculate_rumple_flat_surface(): + chm = np.full((4, 4), 12.0) + rumple = calculate_rumple(chm, (2.0, 3.0)) + assert np.isclose(rumple, 1.0) + + +def test_calculate_rumple_planar_slope_matches_analytical_ratio(): + dx, dy = 2.0, 3.0 + a, b = 0.5, 0.25 + x = np.arange(5) * dx + y = np.arange(4) * dy + chm = a * x[:, None] + b * y[None, :] + + rumple = calculate_rumple(chm, (dx, dy)) + expected = math.sqrt(1.0 + a ** 2 + b ** 2) + assert np.isclose(rumple, expected) + + +def test_calculate_rumple_rough_surface_exceeds_flat_surface(): + flat = np.full((3, 3), 10.0) + rough = np.array([ + [10.0, 10.0, 10.0], + [10.0, 15.0, 10.0], + [10.0, 10.0, 10.0], + ]) + + flat_rumple = calculate_rumple(flat, (1.0, 1.0)) + rough_rumple = calculate_rumple(rough, (1.0, 1.0)) + assert np.isclose(flat_rumple, 1.0) + assert rough_rumple > flat_rumple + + +def test_calculate_rumple_min_height_can_remove_all_valid_patches(): + chm = np.array([ + [3.0, 3.0, 3.0], + [3.0, 1.0, 3.0], + [3.0, 3.0, 3.0], + ]) + rumple = calculate_rumple(chm, (1.0, 1.0), min_height=2.0) + assert np.isnan(rumple) + + +def test_calculate_rumple_returns_nan_when_no_valid_surface_exists(): + chm = np.array([ + [np.nan, np.nan], + [np.nan, np.nan], + ]) + rumple = calculate_rumple(chm, (1.0, 1.0)) + assert np.isnan(rumple) + + +def test_calculate_rumple_invalid_inputs(): + with pytest.raises(ValueError): + calculate_rumple(np.ones((2, 2, 2)), (1.0, 1.0)) + + with pytest.raises(ValueError): + calculate_rumple(np.ones((2, 2)), (0.0, 1.0)) + + with pytest.raises(ValueError): + calculate_rumple(np.ones((2, 2)), (1.0,)) + + # ---------------------------- # Tests for calculate_canopy_cover # ---------------------------- From f4b15cc738108d738a6b754219ffcf2130bd76b2 Mon Sep 17 00:00:00 2001 From: iosefa Date: Mon, 18 May 2026 22:59:55 -1000 Subject: [PATCH 2/4] Add rumple docs to navigation --- mkdocs.yml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/mkdocs.yml b/mkdocs.yml index 83a2b64..902970b 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -42,6 +42,7 @@ nav: - Forest Structure: - usage/forest-structure/intro.md - usage/forest-structure/chm.md + - usage/forest-structure/rumple.md - usage/forest-structure/pad.md - usage/forest-structure/pai.md - usage/forest-structure/fhd.md @@ -66,4 +67,4 @@ extra_css: - overrides/custom.css extra_javascript: - - https://cdnjs.cloudflare.com/ajax/libs/mathjax/3.2.0/es5/tex-mml-chtml.js \ No newline at end of file + - https://cdnjs.cloudflare.com/ajax/libs/mathjax/3.2.0/es5/tex-mml-chtml.js From 42aa95809a7c89be196b6f57bfdda78808cdde67 Mon Sep 17 00:00:00 2001 From: iosefa Date: Mon, 22 Jun 2026 13:36:15 -1000 Subject: [PATCH 3/4] Prepare v0.4.0 release --- README.md | 4 ++-- setup.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index 3c057e7..f58d91f 100644 --- a/README.md +++ b/README.md @@ -14,11 +14,11 @@ ## Overview -PyForestScan is a Python library designed for analyzing and visualizing forest structure using airborne 3D point cloud data. The library helps derive important forest metrics such as Canopy Height, Plant Area Index (PAI), Canopy Cover, Plant Area Density (PAD), and Foliage Height Diversity (FHD). +PyForestScan is a Python library designed for analyzing and visualizing forest structure using airborne 3D point cloud data. The library helps derive important forest metrics such as Canopy Height, Plant Area Index (PAI), Canopy Cover, Plant Area Density (PAD), Foliage Height Diversity (FHD), and Rumple Index. ## Features -- **Forest Metrics**: Calculate and visualize key metrics like Canopy Height, PAI, PAD, and FHD. +- **Forest Metrics**: Calculate and visualize key metrics like Canopy Height, PAI, PAD, FHD, Canopy Cover, and Rumple Index. - **Large Point Cloud Support**: Utilizes efficient data formats such as EPT for large point cloud processing. - **Visualization**: Create 2D and 3D visualizations of forest structure and structural metrics - **Extensibility**: Easily add custom filters and visualization techniques to suit your needs. diff --git a/setup.py b/setup.py index 39ce535..2871fad 100644 --- a/setup.py +++ b/setup.py @@ -5,7 +5,7 @@ setuptools.setup( name="pyforestscan", - version="0.3.8", + version="0.4.0", author="Joseph Emile Honour Percival", author_email="ipercival@gmail.com", description="Analyzing forest structure using aerial LiDAR data", From 82cc8334144ba43f4402d74ac0e64d8d2a1294cb Mon Sep 17 00:00:00 2001 From: iosefa Date: Sun, 28 Jun 2026 21:49:16 +1300 Subject: [PATCH 4/4] Bump version to 0.4.1 --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 2871fad..575a1ba 100644 --- a/setup.py +++ b/setup.py @@ -5,7 +5,7 @@ setuptools.setup( name="pyforestscan", - version="0.4.0", + version="0.4.1", author="Joseph Emile Honour Percival", author_email="ipercival@gmail.com", description="Analyzing forest structure using aerial LiDAR data",