diff --git a/pyforestscan/process.py b/pyforestscan/process.py index b243c11..01c8144 100644 --- a/pyforestscan/process.py +++ b/pyforestscan/process.py @@ -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. @@ -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 @@ -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 @@ -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 @@ -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, @@ -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) @@ -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}")