Skip to content
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
48 changes: 39 additions & 9 deletions pyforestscan/process.py
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,8 @@ def _crop_dtm(dtm_path, tile_min_x, tile_min_y, tile_max_x, tile_max_y):
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,
cover_min_height: float = 2.0, cover_k: float = 0.5) -> None:
cover_min_height: float = 2.0, cover_k: float = 0.5,
skip_existing: bool = False, verbose: bool = False) -> 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.
Expand All @@ -75,6 +76,8 @@ def process_with_tiles(ept_file, tile_size, output_path, metric, voxel_size,
remove_outliers (bool, optional): Whether to remove statistical outliers before calculating metrics. Defaults to False.
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.
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.

Returns:
None
Expand Down Expand Up @@ -126,7 +129,17 @@ def process_with_tiles(ept_file, tile_size, output_path, metric, voxel_size,
tile_max_y = min(max_y, tile_max_y)

if tile_max_x <= tile_min_x or tile_max_y <= tile_min_y:
print(f"Warning: Skipping tile ({i}, {j}) due to invalid spatial extent.")
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

Expand Down Expand Up @@ -164,7 +177,8 @@ def process_with_tiles(ept_file, tile_size, output_path, metric, voxel_size,
tile_points = tile_pipeline.arrays[0]

if tile_points.size == 0:
print(f"Warning: No data in tile ({i}, {j}). Skipping.")
if verbose:
print(f"Warning: No data in tile ({i}, {j}). Skipping.")
pbar.update(1)
continue

Expand All @@ -175,12 +189,29 @@ def process_with_tiles(ept_file, tile_size, output_path, metric, voxel_size,
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]:
print(
f"Warning: Buffer size exceeds CHM dimensions for tile ({i}, {j}). Adjusting buffer size.")
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)

chm = chm[buffer_pixels_y:-buffer_pixels_y, buffer_pixels_x:-buffer_pixels_x]
# 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

core_extent = (
tile_min_x + buffer_x,
Expand All @@ -189,7 +220,6 @@ def process_with_tiles(ept_file, tile_size, output_path, metric, voxel_size,
tile_max_y - buffer_y,
)

result_file = os.path.join(output_path, f"tile_{i}_{j}_chm.tif")
create_geotiff(chm, result_file, srs, core_extent)
elif metric in ["fhd", "pai", "cover"]:
voxels, spatial_extent = assign_voxels(tile_points, voxel_size)
Expand Down Expand Up @@ -243,11 +273,11 @@ def process_with_tiles(ept_file, tile_size, output_path, metric, voxel_size,
)

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

result_file = os.path.join(output_path, f"tile_{i}_{j}_{metric}.tif")
create_geotiff(result, result_file, srs, core_extent)
else:
raise ValueError(f"Unsupported metric: {metric}")
Expand Down
Loading