Merge pull request #39 from iosefa/feat/improve-process-tiles

process: add `skip_existing` and `verbose` options to tile processing

Author: iosefa

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}")