Detect Treetop

This module implements functions for detecting treetop locations from canopy height models (CHMs) derived from lidar data.

DetectionParams dataclass

Parameters for treetop detection algorithms.

Parameters:

Name	Type	Description	Default
detection_method	str	Method to use for treetop detection. Options: "prominence", "vws", "lmf".	'prominence'
pixel_size	float	Spatial resolution of the CHM in the same units as min_tree_distance.	0.25
min_height	float	Minimum height threshold for valid treetops (in same units as CHM).	3.0
min_tree_distance	float	Minimum distance between detected treetops (in same units as CHM).	1.5
prominence_height	float	Height above local background for prominence-based detection (in same units as CHM).	1.0
vws_detection_sigma	float	Sigma for Gaussian smoothing in VWS method. Set to 0 for no smoothing.	2.0
vws_distance_scale	float	Scaling factor for distance in VWS method (in same units as CHM).	0.12
vws_power	float	Power to which height is raised in VWS method for distance scaling.	1.0

Source code in src/phytospatial/lidar/detect_treetop.py

@dataclass
class DetectionParams:
    """
    Parameters for treetop detection algorithms.

    Args:
        detection_method (str): Method to use for treetop detection. Options: "prominence", "vws", "lmf".
        pixel_size (float): Spatial resolution of the CHM in the same units as min_tree_distance.
        min_height (float): Minimum height threshold for valid treetops (in same units as CHM).
        min_tree_distance (float): Minimum distance between detected treetops (in same units as CHM).
        prominence_height (float): Height above local background for prominence-based detection (in same units as CHM).
        vws_detection_sigma (float): Sigma for Gaussian smoothing in VWS method. Set to 0 for no smoothing.
        vws_distance_scale (float): Scaling factor for distance in VWS method (in same units as CHM).
        vws_power (float): Power to which height is raised in VWS method for distance scaling.
    """
    detection_method: str = "prominence"
    pixel_size: float = 0.25
    min_height: float = 3.0         
    min_tree_distance: float = 1.5         
    prominence_height: float = 1.0     
    vws_detection_sigma: float = 2.0   
    vws_distance_scale: float = 0.12   
    vws_power: float = 1.0
    lmf_window_size: int = 5

detect_treetops(chm_input, params=DetectionParams(), tile_mode='auto', tile_size=1024, overlap=64)

Detects treetop locations from a canopy height model (CHM) using specified detection parameters and memory-safe processing strategies.

Steps

1. Utilizes iter_core_halo to stream the CHM in partitioned spatial blocks. Each block includes a core processing area and a halo overlap to eliminate edge artifacts during algorithmic focal filtering.
2. Routes the active block array to the specified detection function (Prominence, VWS, or LMF) based on the DetectionParams configuration.
3. Receives localized row and column peak indices from the algorithmic subroutines.
4. Projects these localized pixel coordinates into global coordinate reference system (CRS) points via the block's Affine transform matrix.
5. Performs spatial filtering using core_box.contains() to discard any detected tops that occurred within the halo overlap region, preventing duplicate stem registrations across neighboring blocks.
6. Yields a standardized dictionary mapping the verified geometric point, associated elevation, and execution metadata back to the orchestrator.

Parameters:

Name	Type	Description	Default
chm_input	Union[str, Path, Raster]	Input CHM raster or path to CHM file.	required
params	DetectionParams	Parameters dictating the active treetop detection algorithm.	DetectionParams()
tile_mode	str	Tiling strategy for processing large rasters ('auto', 'in_memory', 'blocked', 'tiled').	'auto'
tile_size	int	Dimensions of the spatial partitions if streaming is necessary.	1024
overlap	int	Margin of pixels surrounding each block to ensure algorithm continuity.	64

Yields:

Type	Description
Dict[str, Any]	Dict[str, Any]: Structured payload containing 'height' (float), 'det_method' (str), and 'geometry' (shapely.Point).

Source code in src/phytospatial/lidar/detect_treetop.py

def detect_treetops(
    chm_input: Union[str, Path, Raster], 
    params: DetectionParams = DetectionParams(),
    tile_mode: str = "auto",
    tile_size: int = 1024,
    overlap: int = 64
    ) -> Generator[Dict[str, Any], None, None]:
    """
    Detects treetop locations from a canopy height model (CHM) using specified 
    detection parameters and memory-safe processing strategies.

    Steps:
        1. Utilizes `iter_core_halo` to stream the CHM in partitioned spatial blocks. 
           Each block includes a core processing area and a halo overlap to eliminate 
           edge artifacts during algorithmic focal filtering.
        2. Routes the active block array to the specified detection function 
           (Prominence, VWS, or LMF) based on the `DetectionParams` configuration.
        3. Receives localized row and column peak indices from the algorithmic subroutines.
        4. Projects these localized pixel coordinates into global coordinate reference 
           system (CRS) points via the block's Affine transform matrix.
        5. Performs spatial filtering using `core_box.contains()` to discard any 
           detected tops that occurred within the halo overlap region, preventing 
           duplicate stem registrations across neighboring blocks.
        6. Yields a standardized dictionary mapping the verified geometric point, 
           associated elevation, and execution metadata back to the orchestrator.

    Args:
        chm_input (Union[str, Path, Raster]): Input CHM raster or path to CHM file.
        params (DetectionParams): Parameters dictating the active treetop detection algorithm.
        tile_mode (str): Tiling strategy for processing large rasters ('auto', 'in_memory', 'blocked', 'tiled').
        tile_size (int): Dimensions of the spatial partitions if streaming is necessary.
        overlap (int): Margin of pixels surrounding each block to ensure algorithm continuity.

    Yields:
        Dict[str, Any]: Structured payload containing 'height' (float), 'det_method' (str), and 'geometry' (shapely.Point).
    """
    # iter_core_halo is a generator that yields core tiles of the CHM along with their halo (overlap) and the affine transform for georeferencing.
    # We use it to process the CHM in manageable chunks, which allows us to handle larger rasters that may not fit into memory all at once.
    # Without the halo, we might miss treetops that are near the edges of the tiles, so the overlap ensures we can detect those properly.
    for chm, transform, core_box, _ in iter_core_halo(chm_input, tile_mode, tile_size, overlap):

        # We call the appropriate method to detect peaks in the current tile of the CHM based on the specified detection method in the parameters.
        if params.detection_method == "prominence":
            peaks_r, peaks_c = _detect_peaks_prominence(chm, params)
        elif params.detection_method == "vws":
            chm_detect = ndimage.gaussian_filter(chm, sigma=params.vws_detection_sigma) if params.vws_detection_sigma > 0 else chm
            peaks_r, peaks_c = _detect_peaks_vws(
                chm_detect, 
                params.min_tree_distance / params.pixel_size, 
                params.vws_distance_scale / params.pixel_size, 
                params.min_height,
                params.vws_power
            )
        elif params.detection_method == "lmf":
            peaks_r, peaks_c = _detect_peaks_lmf(chm, params)
        else: 
            raise ValueError(f"Unknown detection method: {params.detection_method}")

        # The detected peaks are returned as row and column indices relative to the current tile.
        # We need to convert these to spatial coordinates using the affine transform, and we also filter them based on the core box
        # to ensure we only return valid detections that are not in the halo area.    
        apex_ref = np.column_stack((peaks_r, peaks_c))
        if len(apex_ref) == 0:
            continue

        for r, c in apex_ref:
            top_x, top_y = transform * (c + 0.5, r + 0.5)
            geom = Point(top_x, top_y)

            if core_box is not None and not core_box.contains(geom):
                continue

            # Finally, we yield a metadata dictionary for each detected treetop containing its height, detection method used, and its geometry (Point).    
            yield {
                "height": float(chm[r, c]),
                "det_method": params.detection_method,
                "geometry": geom
            }