Identifying Canopy Gaps Using Morphological Filters

Identifying canopy gaps using morphological filters provides a deterministic, computationally efficient alternative to heuristic height-thresholding or manual digitization. When processing airborne LiDAR-derived Canopy Height Models (CHMs), foresters and Python GIS developers routinely encounter fragmented gap boundaries, false positives from shrub layers, and memory bottlenecks during regional-scale analysis. Morphological image processing operations—specifically opening, closing, and white-top-hat transformations—treat the CHM as a continuous elevation surface, enabling precise isolation of depressions that meet strict structural and spatial criteria. This methodology integrates directly into established Canopy Height Modeling & Terrain Extraction pipelines, where normalized height rasters serve as the foundational input for structural forest metrics and disturbance mapping.

Morphological Theory Applied to Canopy Surfaces

A CHM represents the vertical distance between the first-return canopy surface and the interpolated ground. Ecologically meaningful gaps manifest as localized depressions where height values drop sharply relative to the surrounding crown envelope. Static thresholding (e.g., CHM < 2.0 m) fails to account for edge effects, small canopy protrusions, and terrain-induced slope artifacts. Morphological filters resolve this by applying structuring elements that probe the raster at defined spatial scales.

• Opening (erosion followed by dilation) suppresses isolated canopy peaks and high-frequency noise while preserving broader crown structures.
• Closing (dilation followed by erosion) bridges narrow canopy breaks and fills minor understory depressions.
• White-Top-Hat Transform (original - opened) explicitly extracts depressions that fall below the local canopy envelope. This operation is mathematically equivalent to identifying regions where the CHM is lower than the structurally smoothed surface.

By calibrating the structuring element radius and height differential threshold, developers can isolate gaps that match ecological definitions (e.g., >10 m² area, >3 m height drop). This approach is particularly effective for Forest Gap & Understory Analysis workflows where spatial continuity and scale-aware detection are required.

Production-Ready Python Implementation

The following workflow uses scipy.ndimage for morphological operations and rasterio for windowed I/O handling. It avoids loading entire rasters into memory by leveraging block-based processing with explicit overlap padding to prevent boundary artifacts. For detailed windowed I/O patterns, refer to the Rasterio windowed read/write documentation.

import numpy as np
import rasterio
from rasterio.windows import Window
from scipy.ndimage import grey_closing, label, sum as ndimage_sum
from skimage.morphology import disk

def detect_canopy_gaps(
    chm_path: str,
    output_path: str,
    gap_height_thresh: float = 2.0,
    min_gap_area_m2: float = 15.0,
    struct_radius_m: float = 5.0,
    overlap_px: int = 32,
) -> None:
    """
    Identify canopy gaps with a morphological black top-hat transform
    (grey_closing - chm), which highlights dark valleys — i.e. low spots
    surrounded by taller canopy. Processed in overlapping windows so edge
    pixels never see a truncated neighbourhood.
    """
    with rasterio.open(chm_path) as src:
        if src.count != 1:
            raise ValueError("Input CHM must be a single-band raster.")
        if src.crs is None or src.transform is None:
            raise ValueError("Input CHM must contain valid CRS and transform metadata.")

        res = src.res[0]  # Assume square pixels for simplicity
        struct_elem = disk(max(1, int(round(struct_radius_m / res))))
        area_thresh_px = min_gap_area_m2 / (res ** 2)

        meta = src.meta.copy()
        meta.update(dtype='uint8', count=1, nodata=0)

        with rasterio.open(output_path, 'w', **meta) as dst:
            for _ji, window in src.block_windows(1):
                # Expand window with overlap to prevent edge truncation.
                col_start_src = max(window.col_off - overlap_px, 0)
                row_start_src = max(window.row_off - overlap_px, 0)
                win = Window(
                    col_start_src,
                    row_start_src,
                    min(window.width  + 2 * overlap_px, src.width  - col_start_src),
                    min(window.height + 2 * overlap_px, src.height - row_start_src),
                )

                chm_block = src.read(1, window=win).astype(np.float32)

                # Treat non-positive samples as nodata for top-hat math.
                valid_mask = chm_block > 0
                if not np.any(valid_mask):
                    continue
                chm_block[~valid_mask] = 0

                # Black top-hat: closing(chm) - chm — large where chm is locally low.
                closed       = grey_closing(chm_block, structure=struct_elem)
                gap_surface  = closed - chm_block

                # Pixels are gaps if they sit in a deep depression AND fall below
                # the absolute height threshold (so tall canopy with a small dip
                # is not flagged as ground).
                raw_gap_mask = (
                    (gap_surface >= gap_height_thresh)
                    & (chm_block < gap_height_thresh)
                ).astype(np.uint8)

                # Filter connected components by minimum area.
                labeled, num_features = label(raw_gap_mask)
                if num_features > 0:
                    sizes = ndimage_sum(raw_gap_mask, labeled, range(1, num_features + 1))
                    keep_labels = np.where(sizes >= area_thresh_px)[0] + 1
                    final_mask = np.isin(labeled, keep_labels).astype(np.uint8)
                else:
                    final_mask = np.zeros_like(raw_gap_mask, dtype=np.uint8)

                # Trim the overlap padding before writing back to the un-padded window.
                # The expanded window starts `min(overlap_px, window.row_off)` rows
                # *above* the original — that prefix is the slice we discard.
                row_start = min(overlap_px, window.row_off)
                col_start = min(overlap_px, window.col_off)
                h = min(window.height, final_mask.shape[0] - row_start)
                w = min(window.width,  final_mask.shape[1] - col_start)

                dst.write(
                    final_mask[row_start:row_start + h, col_start:col_start + w],
                    1,
                    window=window,
                )

Spatial Constraints & Parameter Calibration

Morphological gap detection is highly sensitive to spatial resolution and structuring element geometry. Adhere to the following constraints to maintain ecological validity:

1. Projection Requirements: Input CHMs must be in a metric projected coordinate system (e.g., UTM, State Plane). Geographic coordinates (WGS84) will distort structuring element radii and invalidate area calculations.
2. Structuring Element Radius: The struct_radius_m should approximate the average crown radius of the dominant canopy species. Values between 3–8 m are typical for temperate broadleaf and coniferous stands. Oversized elements will merge adjacent gaps; undersized elements will retain canopy noise.
3. Height Threshold Selection: The gap_height_thresh defines the minimum vertical drop required to classify a depression as a gap. For closed-canopy forests, 2.0–3.5 m is standard. In open woodlands or savannas, lower thresholds (1.0–1.5 m) may be appropriate, but require stricter area filtering to exclude understory.
4. Area Filtering: Always convert min_gap_area_m2 to pixel units using the raster’s native resolution. Filtering post-thresholding eliminates isolated pixels and micro-depressions that lack ecological significance.

Targeted Troubleshooting

Integration Notes

Morphological gap masks serve as high-quality inputs for downstream ecological modeling. Binary gap layers can be vectorized for spatial statistics, intersected with species distribution models, or used to parameterize light availability indices. When scaling to landscape-level assessments, combine morphological detection with multi-temporal CHM differencing to track gap dynamics and recovery trajectories. For advanced structural metrics, integrate gap topology with Forest Gap & Understory Analysis frameworks to quantify edge-to-interior ratios and canopy closure gradients.