Process Module

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=8, outlier_multiplier=3.0, cover_min_height=2.0, cover_k=0.5, pai_min_height=1.0, fhd_min_height=0.0, skip_existing=False, verbose=False, thin_radius=None, voxelgrid_cell=None, voxelgrid_mode='first', tile_indices=None, outliers_before_hag=False, rumple_min_height=None, rumple_resolution=None, rumple_method='two_triangle')

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.

Parameters:

Name	Type	Description	Default
ept_file	str	Path to the EPT file containing the point cloud data.	required
tile_size	tuple	Size of each tile as (tile_width, tile_height).	required
output_path	str	Directory where the output files will be saved.	required
metric	str	Metric to compute for each tile ("chm", "fhd", "pai", "cover", or "rumple").	required
voxel_size	tuple	Voxel resolution as (x_resolution, y_resolution, z_resolution).	required
voxel_height	float	Height of each voxel in meters. Required if metric is "fhd", "pai", or "cover".	1
buffer_size	float	Fractional buffer size relative to tile size (e.g., 0.1 for 10% buffer). Defaults to 0.1.	0.1
srs	str	Spatial Reference System for the output. If None, uses SRS from the EPT file.	None
hag	bool	If True, compute Height Above Ground using Delaunay triangulation. Defaults to False.	False
hag_dtm	bool	If True, compute Height Above Ground using a provided DTM raster. Defaults to False.	False
dtm	str	Path to the DTM raster file. Required if hag_dtm is True.	None
bounds	tuple	Spatial bounds to crop the data. Must be of the form ([xmin, xmax], [ymin, ymax], [zmin, zmax]) or ([xmin, xmax], [ymin, ymax]). If None, tiling is done over the entire dataset.	None
interpolation	str or None	Fill missing canopy heights for CHM or rumple using "linear", "cubic", "nearest", or None. For rumple, filling happens before the height mask and surface-area calculation. Other metrics do not use this option. Defaults to None.	None
remove_outliers	bool	Whether to remove statistical outliers before calculating metrics. Defaults to False.	False
outlier_mean_k	int	Number of nearest neighbors used by the statistical outlier filter. Used only when remove_outliers is True. Defaults to 8.	8
outlier_multiplier	float	Standard deviation multiplier used by the statistical outlier filter. Used only when remove_outliers is True. Defaults to 3.0.	3.0
cover_min_height	float	Height threshold (in meters) for canopy cover (used when metric == "cover"). Defaults to 2.0.	2.0
cover_k	float	Beer–Lambert extinction coefficient for canopy cover. Defaults to 0.5.	0.5
pai_min_height	float	Minimum height (m) to integrate PAI. Defaults to 1.0.	1.0
fhd_min_height	float	Minimum height (m) to include in FHD entropy. Defaults to 0.0.	0.0
skip_existing	bool	If True, skip tiles whose output file already exists. Defaults to False.	False
verbose	bool	If True, print warnings for empty/invalid tiles and buffer adjustments. Defaults to False.	False
thin_radius	float or None	If provided (> 0), apply Poisson radius-based thinning per tile before metrics. Units are in the same CRS as the data (e.g., meters). Defaults to None.	None
voxelgrid_cell	float or None	If provided (> 0), apply PDAL voxel-grid downsampling per tile before metrics using cell edge length of voxelgrid_cell. Defaults to None.	None
voxelgrid_mode	str	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".	'first'
tile_indices	set[tuple[int, int]] or None	If provided, process only these zero-based tile indices while preserving the full tiling grid and buffer behavior.	None
outliers_before_hag	bool	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.	False
rumple_min_height	float or None	Treat canopy heights strictly below this threshold as missing before calculating rumple. Defaults to None (use all finite heights). Missing heights follow rumple_method's rules and are not set to zero.	None
rumple_resolution	float or tuple	Required when metric is "rumple". Output cell size, either a scalar or (dx, dy), in coordinate units. Each size must be a positive integer multiple of the corresponding XY voxel_size. The voxel sizes set the fine CHM resolution; rumple_resolution sets the area over which surface and ground areas are summed. Tile sizes must be multiples of the output resolution. Bounds are expanded to whole fine cells, and the output grid is anchored at their northwest corner and padded east and south to whole output cells. Padding contributes no area. Rumple tile indices start at the northwest corner (j=0 is northmost). At least one fine neighboring cell is read around each tile, even when buffer_size=0. With interpolation enabled, use a buffer wide enough to include samples around gaps: filling uses only the buffered tile and can differ from a whole-cloud calculation.	None
rumple_method	str	Surface estimator and missing-data rule. Nonfinite and height-masked CHM values are treated as missing. "two_triangle" (default) splits each 2x2 CHM patch into two triangles and excludes patches with any missing height from both surface and ground area. "jenness" uses eight neighbors per cell: missing center cells are excluded, missing neighbors use the center height locally, and edge heights are repeated at the raster boundary. Cells without supported area produce NoData.	'two_triangle'

Returns:

Type	Description
None	None

Raises:

Type	Description
ValueError	If an unsupported metric is requested, if buffer or voxel sizes are invalid, or required arguments are missing.
FileNotFoundError	If the EPT or DTM file does not exist, or a required file for processing is missing.

Source code in pyforestscan/process.py

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",
                       tile_indices: set[tuple[int, int]] | None = None,
                       outliers_before_hag: bool = False,
                       rumple_min_height: float | None = None,
                       rumple_resolution=None,
                       rumple_method: str = "two_triangle") -> 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.

    Args:
        ept_file (str): Path to the EPT file containing the point cloud data.
        tile_size (tuple): Size of each tile as (tile_width, tile_height).
        output_path (str): Directory where the output files will be saved.
        metric (str): Metric to compute for each tile ("chm", "fhd", "pai", "cover", or "rumple").
        voxel_size (tuple): Voxel resolution as (x_resolution, y_resolution, z_resolution).
        voxel_height (float, optional): Height of each voxel in meters. Required if metric is "fhd", "pai", or "cover".
        buffer_size (float, optional): Fractional buffer size relative to tile size (e.g., 0.1 for 10% buffer). Defaults to 0.1.
        srs (str, optional): Spatial Reference System for the output. If None, uses SRS from the EPT file.
        hag (bool, optional): If True, compute Height Above Ground using Delaunay triangulation. Defaults to False.
        hag_dtm (bool, optional): If True, compute Height Above Ground using a provided DTM raster. Defaults to False.
        dtm (str, optional): Path to the DTM raster file. Required if hag_dtm is True.
        bounds (tuple, optional): Spatial bounds to crop the data. Must be of the form
            ([xmin, xmax], [ymin, ymax], [zmin, zmax]) or ([xmin, xmax], [ymin, ymax]).
            If None, tiling is done over the entire dataset.
        interpolation (str or None, optional): Fill missing canopy heights for CHM
            or rumple using "linear", "cubic", "nearest", or None. For rumple,
            filling happens before the height mask and surface-area calculation.
            Other metrics do not use this option. Defaults to 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.
        fhd_min_height (float, optional): Minimum height (m) to include in FHD entropy. Defaults to 0.0.
        skip_existing (bool, optional): If True, skip tiles whose output file already exists. Defaults to False.
        verbose (bool, optional): If True, print warnings for empty/invalid tiles and buffer adjustments. Defaults to False.
        thin_radius (float or None, optional): If provided (> 0), apply Poisson radius-based thinning per tile before metrics.
            Units are in the same CRS as the data (e.g., meters). Defaults to None.
        voxelgrid_cell (float or None, optional): If provided (> 0), apply PDAL voxel-grid downsampling per tile before metrics
            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.
        rumple_min_height (float or None, optional): Treat canopy heights strictly
            below this threshold as missing before calculating rumple. Defaults to
            None (use all finite heights). Missing heights follow rumple_method's
            rules and are not set to zero.
        rumple_resolution (float or tuple, optional): Required when metric is "rumple".
            Output cell size, either a scalar or (dx, dy), in coordinate units. Each
            size must be a positive integer multiple of the corresponding XY
            voxel_size. The voxel sizes set the fine CHM resolution; rumple_resolution
            sets the area over which surface and ground areas are summed. Tile sizes
            must be multiples of the output resolution. Bounds are expanded to whole
            fine cells, and the output grid is anchored at their northwest corner and
            padded east and south to whole output cells. Padding contributes no area.
            Rumple tile indices start at the northwest corner (j=0 is northmost).
            At least one fine neighboring cell is read around each tile, even when
            buffer_size=0. With interpolation enabled, use a buffer wide enough to
            include samples around gaps: filling uses only the buffered tile and can
            differ from a whole-cloud calculation.
        rumple_method (str, optional): Surface estimator and missing-data rule.
            Nonfinite and height-masked CHM values are treated as missing.
            "two_triangle" (default) splits each 2x2 CHM patch into two triangles
            and excludes patches with any missing height from both surface and
            ground area. "jenness" uses eight neighbors per cell: missing center
            cells are excluded, missing neighbors use the center height locally,
            and edge heights are repeated at the raster boundary. Cells without
            supported area produce NoData.

    Returns:
        None

    Raises:
        ValueError: If an unsupported metric is requested, if buffer or voxel sizes are invalid, or required arguments are missing.
        FileNotFoundError: If the EPT or DTM file does not exist, or a required file for processing is missing.
    """
    if metric not in ["chm", "fhd", "pai", "cover", "rumple"]:
        raise ValueError(f"Unsupported metric: {metric}")

    if metric == "rumple":
        resolution = np.asarray(voxel_size, dtype=float)
        if resolution.shape != (3,) or not np.all(np.isfinite(resolution) & (resolution > 0)):
            raise ValueError("voxel_size must contain three positive, finite sizes (dx, dy, dz)")
        if rumple_resolution is None:
            raise ValueError("rumple_resolution is required for metric='rumple'")
        _, output_resolution, rumple_factors = _validate_rumple_grid_resolution(
            resolution[:2], rumple_resolution,
        )
        _validate_rumple_method(rumple_method)
        tile_dimensions = np.asarray(tile_size, dtype=float)
        if tile_dimensions.shape != (2,):
            raise ValueError("Rumple tile_size must contain a width and height")
        tile_cells = tile_dimensions / output_resolution
        if (not np.all(np.isfinite(tile_cells) & (tile_cells >= 1))
                or not np.allclose(tile_cells, np.round(tile_cells), rtol=0, atol=1e-8)):
            raise ValueError("Rumple tile_size must be positive multiples of rumple_resolution")
        if not np.isfinite(buffer_size) or buffer_size < 0:
            raise ValueError("buffer_size must be finite and nonnegative")
        if rumple_min_height is not None and not np.isfinite(rumple_min_height):
            raise ValueError("rumple_min_height must be finite or None")
        if interpolation not in (None, "linear", "cubic", "nearest"):
            raise ValueError("interpolation must be None, 'linear', 'cubic', or 'nearest'")

    (min_z, max_z) = (None, None)
    if bounds:
        if len(bounds) == 2:
            (min_x, max_x), (min_y, max_y) = bounds
        else:
            (min_x, max_x), (min_y, max_y), (min_z, max_z) = bounds
    else:
        min_x, max_x, min_y, max_y, min_z, max_z = get_bounds_from_ept(ept_file)

    if metric == "rumple":
        requested_bounds = np.asarray([min_x, max_x, min_y, max_y], dtype=float)
        if (not np.all(np.isfinite(requested_bounds))
                or max_x <= min_x or max_y <= min_y):
            raise ValueError("Rumple bounds must be finite and have positive width and height")
        data_min_x, data_max_x, data_min_y, data_max_y = map(float, requested_bounds)
        dx, dy = map(float, resolution[:2])
        grid_dx, grid_dy = map(float, output_resolution)
        fine_min_x = float(np.floor(min_x / dx) * dx)
        fine_max_x = float(np.ceil(max_x / dx) * dx)
        fine_min_y = float(np.floor(min_y / dy) * dy)
        fine_max_y = float(np.ceil(max_y / dy) * dy)
        fine_nx = int(round((fine_max_x - fine_min_x) / dx))
        fine_ny = int(round((fine_max_y - fine_min_y) / dy))
        # Anchor at the CHM's northwest corner. Extra output footprint at the
        # east and south edges has no supporting surface or ground area.
        min_x, max_y = fine_min_x, fine_max_y
        max_x = min_x + int(np.ceil(fine_nx / rumple_factors[0])) * grid_dx
        min_y = max_y - int(np.ceil(fine_ny / rumple_factors[1])) * grid_dy

    if not srs:
        srs = get_srs_from_ept(ept_file)

    num_tiles_x = int(np.ceil((max_x - min_x) / tile_size[0]))
    num_tiles_y = int(np.ceil((max_y - min_y) / tile_size[1]))
    total_tiles = num_tiles_x * num_tiles_y

    if not os.path.exists(output_path):
        os.makedirs(output_path)

    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

                # Read neighboring points before calculating buffered metrics.
                if metric in ["chm", "pai", "cover", "rumple"]:
                    current_buffer_size = buffer_size
                else:
                    current_buffer_size = 0.0

                buffer_x = current_buffer_size * tile_size[0]
                buffer_y = current_buffer_size * tile_size[1]
                if metric == "rumple":
                    # Read complete neighboring cells before computing area.
                    buffer_x = float(max(dx, np.ceil(buffer_x / dx) * dx))
                    buffer_y = float(max(dy, np.ceil(buffer_y / dy) * dy))
                if metric == "rumple":
                    core_min_x = min_x + i * tile_size[0]
                    core_max_x = min(max_x, core_min_x + tile_size[0])
                    core_max_y = max_y - j * tile_size[1]
                    core_min_y = max(min_y, core_max_y - tile_size[1])
                    core_extent = (core_min_x, core_max_x, core_min_y, core_max_y)
                    halo_min_x = max(fine_min_x, core_min_x - buffer_x)
                    halo_max_x = min(fine_max_x, core_max_x + buffer_x)
                    halo_min_y = max(fine_min_y, core_min_y - buffer_y)
                    halo_max_y = min(fine_max_y, core_max_y + buffer_y)
                    tile_min_x = max(data_min_x, halo_min_x)
                    tile_max_x = min(data_max_x, halo_max_x)
                    tile_min_y = max(data_min_y, halo_min_y)
                    tile_max_y = min(data_max_y, halo_max_y)
                else:
                    tile_min_x = max(min_x, min_x + i * tile_size[0] - buffer_x)
                    tile_max_x = min(max_x, min_x + (i + 1) * tile_size[0] + buffer_x)
                    tile_min_y = max(min_y, min_y + j * tile_size[1] - buffer_y)
                    tile_max_y = min(max_y, min_y + (j + 1) * tile_size[1] + buffer_y)

                if tile_max_x <= tile_min_x or tile_max_y <= tile_min_y:
                    if verbose:
                        print(f"Warning: Skipping tile ({i}, {j}) due to invalid spatial extent.")
                    pbar.update(1)
                    continue

                if metric == "chm":
                    result_file = os.path.join(output_path, f"tile_{i}_{j}_chm.tif")
                else:
                    result_file = os.path.join(output_path, f"tile_{i}_{j}_{metric}.tif")

                if skip_existing and os.path.isfile(result_file):
                    pbar.update(1)
                    continue

                if min_z is not None and max_z is not None:
                    tile_bounds = ([tile_min_x, tile_max_x], [tile_min_y, tile_max_y], [min_z, max_z])
                else:
                    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:
                    if not dtm or not os.path.isfile(dtm):
                        raise FileNotFoundError(f"DTM file is required for HAG calculation using DTM: {dtm}")
                    cropped_dtm_path = _crop_dtm(
                        dtm,
                        tile_min_x, tile_min_y,
                        tile_max_x, tile_max_y
                    )
                    tile_pipeline_stages.append(_hag_raster(cropped_dtm_path))
                base_pipeline = {
                    "type": "readers.ept",
                    "filename": ept_file,
                    "bounds": f"{tile_bounds}",
                }
                tile_pipeline_json = {
                    "pipeline": [base_pipeline] + tile_pipeline_stages
                }

                tile_pipeline = pdal.Pipeline(json.dumps(tile_pipeline_json))
                tile_pipeline.execute()

                # Extract points from pipeline output safely
                arrays = tile_pipeline.arrays if hasattr(tile_pipeline, "arrays") else []
                if not arrays or arrays[0].size == 0:
                    if verbose:
                        print(f"Warning: No data in tile ({i}, {j}). Skipping.")
                    pbar.update(1)
                    continue

                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]

                # Optional thinning before metrics
                if thin_radius is not None and thin_radius > 0:
                    thinned = downsample_poisson([tile_points], thin_radius=thin_radius)
                    tile_points = thinned[0] if thinned else tile_points

                    if tile_points.size == 0:
                        if verbose:
                            print(f"Warning: Tile ({i}, {j}) empty after thinning. Skipping.")
                        pbar.update(1)
                        continue

                if voxelgrid_cell is not None and voxelgrid_cell > 0:
                    try:
                        vthinned = downsample_voxel([tile_points], cell=float(voxelgrid_cell), mode=voxelgrid_mode)
                    except Exception as e:
                        # Fail soft per tile to keep processing going
                        if verbose:
                            print(f"Warning: Voxel-grid downsampling failed for tile ({i}, {j}): {e}. Proceeding without.")
                        vthinned = None
                    if vthinned:
                        tile_points = vthinned[0]
                    if tile_points.size == 0:
                        if verbose:
                            print(f"Warning: Tile ({i}, {j}) empty after voxel-grid thinning. Skipping.")
                        pbar.update(1)
                        continue

                buffer_pixels_x = int(np.ceil(buffer_x / voxel_size[0]))
                buffer_pixels_y = int(np.ceil(buffer_y / voxel_size[1]))

                if metric == "rumple":
                    valid = (
                        np.isfinite(tile_points['X']) & np.isfinite(tile_points['Y']) &
                        np.isfinite(tile_points['HeightAboveGround']) &
                        (tile_points['HeightAboveGround'] >= 0)
                    )
                    tile_points = tile_points[valid]
                    if tile_points.size == 0:
                        pbar.update(1)
                        continue

                    # Bin against one fixed global canvas rather than deriving
                    # each tile's CHM footprint from its available samples.
                    # Internal grid-line points belong to the east/north cell;
                    # the very last outer edge belongs to the last cell.
                    global_x = np.floor((tile_points['X'] - fine_min_x) / dx).astype(int)
                    global_y = np.floor((tile_points['Y'] - fine_min_y) / dy).astype(int)
                    global_x[tile_points['X'] == fine_max_x] = fine_nx - 1
                    global_y[tile_points['Y'] == fine_max_y] = fine_ny - 1
                    global_y = fine_ny - 1 - global_y
                    halo_x0 = int(round((halo_min_x - fine_min_x) / dx))
                    halo_y0 = int(round((fine_max_y - halo_max_y) / dy))
                    halo_shape = (int(round((halo_max_x - halo_min_x) / dx)),
                                  int(round((halo_max_y - halo_min_y) / dy)))
                    local_x, local_y = global_x - halo_x0, global_y - halo_y0
                    in_halo = ((local_x >= 0) & (local_x < halo_shape[0]) &
                               (local_y >= 0) & (local_y < halo_shape[1]))
                    if not np.any(in_halo):
                        pbar.update(1)
                        continue
                    chm = np.full(halo_shape, -np.inf)
                    np.maximum.at(chm, (local_x[in_halo], local_y[in_halo]),
                                  tile_points['HeightAboveGround'][in_halo])
                    chm[~np.isfinite(chm)] = np.nan
                    if interpolation is not None:
                        observed = np.isfinite(chm)
                        if np.any(~observed):
                            try:
                                chm[~observed] = griddata(
                                    np.argwhere(observed) * resolution[:2],
                                    chm[observed],
                                    np.argwhere(~observed) * resolution[:2],
                                    method=interpolation,
                                )
                            except QhullError:
                                # Linear/cubic need enough non-collinear samples.
                                # Leave unsupported heights missing.
                                pass

                    surface, ground = _rumple_cell_areas(
                        chm, (dx, dy), min_height=rumple_min_height,
                        method=rumple_method,
                    )
                    core_shape = (int(round((core_max_x - core_min_x) / dx)),
                                  int(round((core_max_y - core_min_y) / dy)))
                    core_x0 = int(round((core_min_x - halo_min_x) / dx))
                    core_y0 = int(round((halo_max_y - core_max_y) / dy))
                    overlap_x = min(core_shape[0], halo_shape[0] - core_x0)
                    overlap_y = min(core_shape[1], halo_shape[1] - core_y0)
                    totals = []
                    fx, fy = rumple_factors
                    output_shape = (core_shape[0] // fx, core_shape[1] // fy)
                    for areas in (surface, ground):
                        # Aggregate only supported fine cells. An output cell
                        # may be much larger than the entire CHM, so do not
                        # allocate its padded footprint at fine resolution.
                        supported = areas[
                            core_x0:core_x0 + overlap_x, core_y0:core_y0 + overlap_y,
                        ]
                        summed = np.add.reduceat(supported, np.arange(0, overlap_x, fx), axis=0)
                        summed = np.add.reduceat(summed, np.arange(0, overlap_y, fy), axis=1)
                        total = np.zeros(output_shape)
                        total[:summed.shape[0], :summed.shape[1]] = summed
                        totals.append(total)
                    result = np.full(totals[0].shape, np.nan)
                    np.divide(totals[0], totals[1], out=result, where=totals[1] > 0)
                    create_geotiff(result, result_file, srs, core_extent)
                elif metric == "chm":
                    chm, extent = calculate_chm(tile_points, voxel_size, interpolation=interpolation)

                    if buffer_pixels_x * 2 >= chm.shape[1] or buffer_pixels_y * 2 >= chm.shape[0]:
                        if verbose:
                            print(
                                f"Warning: Buffer size exceeds CHM dimensions for tile ({i}, {j}). Adjusting buffer size.")
                        buffer_pixels_x = max(0, chm.shape[1] // 2 - 1)
                        buffer_pixels_y = max(0, chm.shape[0] // 2 - 1)

                    # Safe crop: avoid Python's -0 slice behavior producing empty arrays.
                    if buffer_pixels_x < 0 or buffer_pixels_y < 0:
                        raise ValueError("Computed negative buffer pixels; check voxel and buffer sizes.")

                    start_x = buffer_pixels_x
                    end_x = chm.shape[1] - buffer_pixels_x if buffer_pixels_x > 0 else chm.shape[1]
                    start_y = buffer_pixels_y
                    end_y = chm.shape[0] - buffer_pixels_y if buffer_pixels_y > 0 else chm.shape[0]

                    # If cropping would yield an empty array due to small tiles, skip cropping
                    if end_x <= start_x or end_y <= start_y:
                        # Fall back to no buffer crop for this tile
                        core = chm
                    else:
                        core = chm[start_y:end_y, start_x:end_x]

                    chm = core

                    # Derive core extent from CHM grid extent and actual pixel crop to avoid stretching
                    dx, dy = voxel_size[0], voxel_size[1]
                    grid_extent = extent  # [x_min, x_max, y_min, y_max] as returned by calculate_chm
                    core_extent = (
                        grid_extent[0] + buffer_pixels_x * dx,
                        grid_extent[1] - buffer_pixels_x * dx,
                        grid_extent[2] + buffer_pixels_y * dy,
                        grid_extent[3] - buffer_pixels_y * dy,
                    )

                    create_geotiff(chm, result_file, srs, core_extent)
                elif metric in ["fhd", "pai", "cover"]:
                    voxels, spatial_extent = assign_voxels(tile_points, voxel_size)

                    if metric == "fhd":
                        if voxel_size[-1] <= 0:
                            raise ValueError("voxel_size z-resolution must be > 0 for FHD computation.")
                        result = calculate_fhd(
                            voxels,
                            voxel_height=voxel_size[-1],
                            min_height=fhd_min_height,
                        )
                    elif metric == "pai":
                        if not voxel_height:
                            raise ValueError(f"voxel_height is required for metric {metric}")

                        pad = calculate_pad(voxels, voxel_size[-1])

                        if np.all(pad == 0):
                            result = np.zeros((pad.shape[0], pad.shape[1]))
                        else:
                            # Guard against empty integration range when top height < default min_height
                            effective_max_height = pad.shape[2] * voxel_size[-1]
                            default_min_height = pai_min_height
                            if default_min_height >= effective_max_height:
                                result = np.zeros((pad.shape[0], pad.shape[1]))
                            else:
                                result = calculate_pai(pad, voxel_height, min_height=pai_min_height)
                    elif metric == "cover":
                        if not voxel_height:
                            raise ValueError(f"voxel_height is required for metric {metric}")

                        pad = calculate_pad(voxels, voxel_size[-1])
                        if np.all(pad == 0):
                            result = np.zeros((pad.shape[0], pad.shape[1]))
                        else:
                            result = calculate_canopy_cover(
                                pad,
                                voxel_height=voxel_height,
                                min_height=cover_min_height,
                                max_height=None,
                                k=cover_k,
                            )

                    if current_buffer_size > 0:
                        if buffer_pixels_x * 2 >= result.shape[1] or buffer_pixels_y * 2 >= result.shape[0]:
                            if verbose:
                                print(
                                    f"Warning: Buffer size exceeds {metric.upper()} dimensions for tile ({i}, {j}). "
                                    f"Adjusting buffer size."
                                )
                            buffer_pixels_x = max(0, result.shape[1] // 2 - 1)
                            buffer_pixels_y = max(0, result.shape[0] // 2 - 1)

                        # Safe crop (avoid -0 slicing)
                        start_x = buffer_pixels_x
                        end_x = result.shape[1] - buffer_pixels_x if buffer_pixels_x > 0 else result.shape[1]
                        start_y = buffer_pixels_y
                        end_y = result.shape[0] - buffer_pixels_y if buffer_pixels_y > 0 else result.shape[0]

                        if end_x > start_x and end_y > start_y:
                            result = result[start_y:end_y, start_x:end_x]

                    # Derive core extent from voxel grid extent and applied pixel crop to avoid stretching
                    dx, dy = voxel_size[0], voxel_size[1]
                    core_extent = (
                        spatial_extent[0] + buffer_pixels_x * dx,
                        spatial_extent[1] - buffer_pixels_x * dx,
                        spatial_extent[2] + buffer_pixels_y * dy,
                        spatial_extent[3] - buffer_pixels_y * dy,
                    )

                    if core_extent[1] <= core_extent[0] or core_extent[3] <= core_extent[2]:
                        if verbose:
                            print(f"Warning: Invalid core extent for tile ({i}, {j}): {core_extent}. Skipping.")
                        pbar.update(1)
                        continue

                    create_geotiff(result, result_file, srs, core_extent)
                else:
                    raise ValueError(f"Unsupported metric: {metric}")

                pbar.update(1)