Delineate Crown

This module implements functions for delineating tree crowns from canopy height models and treetop locations.

DelineationParams dataclass

Parameters for tree crown delineation methods.

Parameters:

Name	Type	Description	Default
delineation_method	str	Method to use for delineation ("watershed" or "region_growing").	'watershed'
pixel_size	float	Resolution of the input canopy height model in meters.	0.25
min_height	float	Minimum height threshold for including pixels in crowns.	3.0
watershed_sigma	float	Sigma value for Gaussian smoothing in watershed method.	0.5
max_crown_radius	float	Maximum expected crown radius in meters for region growing method.	10.0
apex_inclusion	float	Minimum relative height of pixels to be included in the crown compared to the apex height.	0.45
crown_threshold	float	Minimum relative height of pixels to be included in the crown compared to the current average height of the crown during region growing.	0.55

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

@dataclass
class DelineationParams:
    """
    Parameters for tree crown delineation methods.

    Args:
        delineation_method: Method to use for delineation ("watershed" or "region_growing").
        pixel_size: Resolution of the input canopy height model in meters.
        min_height: Minimum height threshold for including pixels in crowns.
        watershed_sigma: Sigma value for Gaussian smoothing in watershed method.
        max_crown_radius: Maximum expected crown radius in meters for region growing method.
        apex_inclusion: Minimum relative height of pixels to be included in the crown compared to the apex height.
        crown_threshold: Minimum relative height of pixels to be included in the crown compared to the 
                        current average height of the crown during region growing.
    """
    delineation_method: str = "watershed"
    pixel_size: float = 0.25
    min_height: float = 3.0
    watershed_sigma: float = 0.5
    max_crown_radius: float = 10.0
    apex_inclusion: float = 0.45
    crown_threshold: float = 0.55
    smoothing_sigma: float = 0.5
    max_canopy_spread: float = 10.0

delineate_crowns(chm_input, treetops, params=DelineationParams(), tile_mode='auto', tile_size=1024, overlap=64)

Delineates tree crowns from a canopy height model (CHM) and treetop locations using specified parameters and processing strategy.

Parameters:

Name	Type	Description	Default
chm_input	Union[str, Path, Raster]	Input canopy height model, either as a file path or a Raster object.	required
treetops	Vector	Vector layer containing the treetop locations as Point geometries.	required
params	DelineationParams	DelineationParams object containing parameters for the delineation methods.	DelineationParams()
tile_mode	str	Tiling strategy for processing large rasters ('auto', 'in_memory', 'blocked', 'tiled').	'auto'
tile_size	int	Size of tiles to use if tiling is necessary.	1024
overlap	int	Number of pixels to overlap between tiles to avoid edge effects.	64

Yields:

Type	Description
Dict[str, Any]	Dict[str, Any]: A dictionary containing the tree ID, height, delineation method used, and geometry of each detected tree crown as a Polygon.

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

def delineate_crowns(
    chm_input: Union[str, Path, Raster], 
    treetops: Vector, 
    params: DelineationParams = DelineationParams(),
    tile_mode: str = "auto",
    tile_size: int = 1024,
    overlap: int = 64
    ) -> Generator[Dict[str, Any], None, None]:
    """
    Delineates tree crowns from a canopy height model (CHM) and treetop locations using specified parameters and processing strategy.

    Args:
        chm_input: Input canopy height model, either as a file path or a Raster object.
        treetops: Vector layer containing the treetop locations as Point geometries.
        params: DelineationParams object containing parameters for the delineation methods.
        tile_mode: Tiling strategy for processing large rasters ('auto', 'in_memory', 'blocked', 'tiled').
        tile_size: Size of tiles to use if tiling is necessary.
        overlap: Number of pixels to overlap between tiles to avoid edge effects.

    Yields:
        Dict[str, Any]: A dictionary containing the tree ID, height, delineation method used, and geometry of each detected tree crown as a Polygon.
    """
    # We first check the coordinate reference system (CRS) of the input CHM and treetop vector layer. 
    # If they differ, we reproject the treetop geometries to match the CHM's CRS to ensure spatial alignment during processing.
    raster_crs = chm_input.crs if isinstance(chm_input, Raster) else None

    # We extract the treetop geometries (points) and their associated attributes into a GeoDataFrame
    tt_gdf = treetops.data
    if tt_gdf.empty:
        return

    # If the treetops and raster have different CRSs, we reproject the treetop geometries to match the raster's CRS.
    # This is much more efficient than reprojecting the raster.
    if raster_crs and treetops.crs != raster_crs:
        tt_gdf = to_crs(treetops, raster_crs, inplace=False).data

    sindex = tt_gdf.sindex

    # We use the iter_core_halo function to process the CHM in manageable tiles, 
    # which allows us to handle larger rasters that may not fit into memory all at once.
    for chm, transform, core_box, read_box in iter_core_halo(chm_input, tile_mode, tile_size, overlap):

        # We compare the spatial index of the treetop points with the current tile's bounding box 
        # to quickly identify which treetops fall within the tile (including the halo).
        if sindex is not None and read_box is not None:
            possible_matches_index = list(sindex.intersection(read_box.bounds))
            local_trees = tt_gdf.iloc[possible_matches_index]
            local_trees = local_trees[local_trees.intersects(read_box)]
        else:
            local_trees = tt_gdf

        if local_trees.empty:
            continue # pass if no treetops in this tile

        # We then convert the spatial coordinates of the local treetops to row and column indices relative to the current tile using the affine transform.
        x_coords = local_trees.geometry.x.values
        y_coords = local_trees.geometry.y.values
        cols, rows = (~transform) * (x_coords, y_coords)

        # We filter out any treetops that fall outside the bounds of the current tile (including the halo)
        valid_mask = (rows >= 0) & (rows < chm.shape[0]) & (cols >= 0) & (cols < chm.shape[1])
        if not np.any(valid_mask):
            continue

        valid_rows = rows[valid_mask].astype(int)
        valid_cols = cols[valid_mask].astype(int)
        apex_p = np.column_stack((valid_rows, valid_cols))

        # We also extract the tree IDs and geometries for the valid treetops to associate them with the delineated crowns later on.
        tree_ids = local_trees['tree_id'].values[valid_mask] if 'tree_id' in local_trees.columns else local_trees.index.values[valid_mask]
        apex_geoms = local_trees.geometry.values[valid_mask]

        # Finally, we call the appropriate delineation method (watershed or region growing)
        if params.delineation_method == "watershed":
            final_labels = _run_watershed(chm, apex_p, params)
        elif params.delineation_method == "region_growing":
            final_labels = _run_region_growing(chm, apex_p, params)
        else: 
            raise ValueError(f"Unknown delineation method: {params.delineation_method}")

        # We use the shapes function from rasterio to convert the labeled array of delineated crowns into vector geometries (Polygons),
        # and we associate each polygon with the corresponding tree ID and height information for output.
        for geometry, value in shapes(final_labels, mask=final_labels > 0, transform=transform):
            idx = int(value) - 1
            apex_geom = apex_geoms[idx]

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

            region_mask = (final_labels == value)

            # We yield a metadata dictionary for each delineated crown containing the tree ID, height, delineation method and  geometry of the crown (Polygon).
            # the height of the tree can be approximated as the maximum CHM value within the delineated crown region,
            # which gives us an estimate of the tree's height.
            yield {
                "tree_id": int(tree_ids[idx]), 
                "height": float(chm[region_mask].max()),
                "del_method": params.delineation_method,
                "geometry": shape(geometry)
            }