Identifying Canopy Gaps Using Morphological Filters (scipy.ndimage Black Top-Hat)

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,
                )

import numpy as np
import rasterio

with rasterio.open("chm.tif") as chm, rasterio.open("gaps.tif") as gaps:
    # 1. Grid alignment — gaps must inherit the CHM transform and shape.
    assert chm.transform == gaps.transform, "Output is not co-registered with the CHM."
    assert chm.shape == gaps.shape, "Output grid differs from input grid."

    mask = gaps.read(1)
    # 2. Strictly binary.
    assert set(np.unique(mask)).issubset({0, 1}), "Mask is not binary."

    # 3. Sanity-check the gap fraction.
    gap_fraction = mask.mean()
    print(f"Gap fraction: {gap_fraction:.1%}")
    assert gap_fraction < 0.40, "Gap fraction implausibly high — re-check struct_radius_m."