diff --git a/CLAUDE.md b/CLAUDE.md index 0a19c45..5c6ba44 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -16,9 +16,8 @@ All commands run inside Docker. Use `./run.sh` as the primary interface. ./run.sh -g -r 0.2 # High resolution (0.2m/px) ./run.sh --test # Run unit tests ./run.sh -g --file LHD_FXX_1000_6882_PTS_LAMB93_IGN69.copc # Single file -./run.sh --ground-classification pmf # Force PMF ground classification -./run.sh -g --keep-tif # Keep intermediate TIFF files (default: kept) -./run.sh -g --no-keep-tif # Delete intermediate TIFF files +./run.sh --ground-classification csf # Force CSF ground classification (complex terrain) +./run.sh -g --keep-tif # Keep TIFF files (allows WebP regeneration without recalculating DTM) ./run.sh # Print help (no args) ``` @@ -34,7 +33,7 @@ docker run --rm --gpus all -v $(pwd)/input:/data/input:ro -v $(pwd)/output:/data - **`cli.py`** — argparse + logging setup. Entry point via `python -m lidar_pipeline`. - **`pipeline.py`** — `LidarArchaeoPipeline` orchestrator. `VIZ_STEPS` registry maps names to generate functions. `FilePrefixFilter` for parallel logging. Creates `SharedDEM` once per file and passes it to all visualizations. -- **`dtm.py`** — PDAL ground classification (SMRF/PMF/CSF + auto-detection) and DTM generation via scipy `binned_statistic_2d`. +- **`dtm.py`** — PDAL ground classification (SMRF/CSF + auto-detection) and DTM generation via scipy `binned_statistic_2d`. - **`visualizations.py`** — 15 `generate_*` functions + 2 IGN overlay lambdas. All take `(dem_file, basename, vis_dir, resolution, shared=None)` and return a TIF path or None. `SharedDEM` class pre-computes gradient, NaN mask, LRM to avoid redundant I/O and computation. - **`gpu.py`** — CuPy/numpy abstraction: `HAS_GPU`, `to_gpu()`, `to_cpu()`, `xp_gaussian_filter()`, `xp_uniform_filter()`, `xp_minimum_filter()`, `gpu_cleanup()`. Falls back to CPU gracefully. - **`ign.py`** — IGN WMTS tile download + overlay generation for orthophoto and topographic maps. @@ -60,11 +59,11 @@ Three places must be updated: ### Ground classification Auto-detection in `dtm.py` `detect_ground_method()`: -- Single-return ratio > 0.6 → PMF (urban terrain) +- Single-return ratio > 0.6 → CSF (urban terrain, cloth simulation) - Height std > 30m → CSF (complex/mountainous terrain) - Default → SMRF (natural terrain) -Override with `--ground-classification {auto,smrf,pmf,csf}`. +Override with `--ground-classification {auto,smrf,csf}`. ### NaN handling @@ -83,6 +82,6 @@ Uses `ProcessPoolExecutor` with `'spawn'` start method (required for CUDA). Each - **Language**: UI messages and comments in French. Code identifiers in English. - **Logging**: Use `logger = logging.getLogger("lidar")`. Prefix per-file logs via `_file_filter.basename`. - **GPU pattern**: `arr_gpu = to_gpu(arr)` → compute → `result = to_cpu(arr_gpu)` → `gpu_cleanup()` between visualizations. -- **Output format**: Visualizations saved as WebP. DTM TIFF kept by default for reuse (use `--no-keep-tif` to delete). `--force` regenerates WebPs without re-classifying if DTM exists. No COGs or viewer — only WebP + PDF report remain. +- **Output format**: Visualizations saved as WebP. TIFF intermediates deleted by default. Use `--keep-tif` to keep DTM+TIF for WebP regeneration with `--force`. No COGs or viewer. - **Compression**: TIF intermediates use `deflate` compression (faster than LZW for float32 data). - **Tests**: Run only inside Docker via `./run.sh --test`. Synthetic DEM fixture in `tests/conftest.py`. \ No newline at end of file diff --git a/lidar_pipeline/cli.py b/lidar_pipeline/cli.py index 66800dd..8587892 100644 --- a/lidar_pipeline/cli.py +++ b/lidar_pipeline/cli.py @@ -118,15 +118,15 @@ Exemples: help="Reclassifier le sol même si le fichier .las existe déjà" ) parser.add_argument( - "--no-keep-tif", + "--keep-tif", action="store_true", - help="Supprimer les fichiers TIFF intermédiaires après conversion WebP (par défaut: conservés)" + help="Conserver les fichiers TIFF (DTM + visualisations) pour pouvoir régénérer les WebP sans recalculer" ) parser.add_argument( "--ground-classification", - choices=["auto", "smrf", "pmf", "csf"], + choices=["auto", "smrf", "csf"], default="auto", - help="Méthode de classification du sol : auto (détection), smrf, pmf, csf (défaut: auto)" + help="Méthode de classification du sol : auto (détection), smrf, csf (défaut: auto)" ) parser.add_argument( "--file", @@ -176,7 +176,7 @@ Exemples: force=args.force, ground_method=args.ground_classification, force_classify=args.force_classification, - keep_tif=not args.no_keep_tif, + keep_tif=args.keep_tif, ) # If --file is specified, process only matching files diff --git a/lidar_pipeline/dtm.py b/lidar_pipeline/dtm.py index 533f349..7dd0ddd 100644 --- a/lidar_pipeline/dtm.py +++ b/lidar_pipeline/dtm.py @@ -1,6 +1,6 @@ """DTM generation from classified LiDAR point clouds. -Handles ground classification via PDAL/SMRF and DTM rasterisation +Handles ground classification via PDAL (SMRF or CSF) and DTM rasterisation using scipy binned_statistic_2d. Zones without LiDAR data remain as NaN. """ @@ -27,13 +27,13 @@ def _create_ground_pipeline(input_laz, output_las, method): 1. Reset Classification to 0 2. ELM (Extended Local Minimum) — mark low outliers as noise (Classification=7) 3. Statistical outlier removal - 4. Ground classification (SMRF/PMF/CSF) + 4. Ground classification (SMRF or CSF) 5. Extract ground points (Classification=2) Args: input_laz: Path to input LAZ/LAS file. output_las: Path to output classified LAS file. - method: Ground classification method ('smrf', 'pmf', or 'csf'). + method: Ground classification method ('smrf' or 'csf'). Returns: JSON string of the PDAL pipeline. @@ -84,15 +84,6 @@ def _create_ground_pipeline(input_laz, output_las, method): "threshold": 0.5, "scalar": 1.25 } - elif method == 'pmf': - ground_step = { - "type": "filters.pmf", - "max_window": 33, - "slope": 0.15, - "max_distance": 2.5, - "initial_distance": 0.15, - "cell_size": 1.0 - } elif method == 'csf': ground_step = { "type": "filters.csf", @@ -128,11 +119,6 @@ def create_smrf_pipeline(input_laz, output_las): return _create_ground_pipeline(input_laz, output_las, 'smrf') -def create_pmf_pipeline(input_laz, output_las): - """Create a PDAL pipeline JSON for PMF ground classification.""" - return _create_ground_pipeline(input_laz, output_las, 'pmf') - - def create_csf_pipeline(input_laz, output_las): """Create a PDAL pipeline JSON for CSF ground classification.""" return _create_ground_pipeline(input_laz, output_las, 'csf') @@ -141,9 +127,9 @@ def create_csf_pipeline(input_laz, output_las): def detect_ground_method(laz_file): """Detect the best ground classification method based on point cloud statistics. - Auto-selects between SMRF (natural terrain) and PMF (urban) only. - CSF is available only via --ground-classification csf (slow but robust - on complex terrain). + Auto-selects between SMRF and CSF: + - SMRF: fast, robust for most natural terrain (PDAL recommended default) + - CSF: cloth simulation, better for complex/urban terrain Falls back to SMRF if the file cannot be read or attributes are missing. @@ -151,7 +137,7 @@ def detect_ground_method(laz_file): laz_file: Path to input LAZ/LAS file. Returns: - String: 'smrf', 'pmf', or 'csf' + String: 'smrf' or 'csf' """ import laspy @@ -182,13 +168,16 @@ def detect_ground_method(laz_file): f"ratio_retours_uniques={single_return_ratio:.2f}, " f"écart_Z={z_std:.1f}m, amplitude_Z={z_range:.1f}m") - # Decision logic (auto selects between SMRF and PMF only): - # - High single-return ratio (>0.6) → urban (buildings, roads) → PMF + # Decision logic: + # - High single-return ratio (>0.6) → urban (buildings, roads) → CSF (cloth simulation) + # - High elevation variance (>30m) → complex/mountainous terrain → CSF # - Default → SMRF (fast, robust for most natural terrain) - # Note: CSF is available only via --ground-classification csf (slow but robust on complex terrain) if single_return_ratio > 0.6: - method = 'pmf' + method = 'csf' reason = f"ratio retours uniques={single_return_ratio:.2f} > 0.6 → milieu urbain" + elif z_std > 30: + method = 'csf' + reason = f"écart_Z={z_std:.1f}m > 30m → terrain complexe" else: method = 'smrf' reason = f"terrain naturel standard" @@ -203,7 +192,7 @@ def classify_ground(laz_file, temp_dir, method='auto', force=False): Args: laz_file: Path to input LAZ/LAS file. temp_dir: Directory for temporary files (pipeline.json, ground.las). - method: Ground classification method ('auto', 'smrf', 'pmf', or 'csf'). + method: Ground classification method ('auto', 'smrf', or 'csf'). force: If True, reclassify even if output file already exists. Returns: diff --git a/lidar_pipeline/pipeline.py b/lidar_pipeline/pipeline.py index 82d462f..de6a72f 100644 --- a/lidar_pipeline/pipeline.py +++ b/lidar_pipeline/pipeline.py @@ -78,8 +78,8 @@ VIZ_STEPS = [ ('curvature', generate_curvature), ('svf', generate_svf), ('lrm', generate_lrm), - ('pos_open', lambda d, b, v, r: generate_openness(d, b, v, r, positive=True)), - ('neg_open', lambda d, b, v, r: generate_openness(d, b, v, r, positive=False)), + ('pos_open', lambda d, b, v, r, shared=None: generate_openness(d, b, v, r, positive=True, shared=shared)), + ('neg_open', lambda d, b, v, r, shared=None: generate_openness(d, b, v, r, positive=False, shared=shared)), ('mslrm', generate_mslrm), ('tpi', generate_tpi), ('sailore', generate_sailore), @@ -323,9 +323,6 @@ class LidarArchaeoPipeline: t_pdf = time.time() - t4 logger.info(f" ✓ Rapport PDF terminé ({t_pdf:.1f}s)") - # Step 5: Keep DTM TIF for reuse (regenerating WebPs skips classification) - # Use --no-keep-tif to delete DTM files after processing - t_total = time.time() - t_start logger.info(f"✓ {basename} terminé en {t_total:.1f}s") logger.debug(f" Détails: classification={t_classif:.1f}s, DTM={t_dtm:.1f}s, PDF={t_pdf:.1f}s") diff --git a/lidar_pipeline/visualizations.py b/lidar_pipeline/visualizations.py index a4121a0..b0dff01 100644 --- a/lidar_pipeline/visualizations.py +++ b/lidar_pipeline/visualizations.py @@ -100,8 +100,18 @@ def _filter_nanaware_from_filled(shared, filter_func, *args, **kwargs): return result -def _save_tif(output_path, data, transform, crs, dtype='float32', count=1, nodata=None): - """Helper to save a 2D or 3D array as GeoTIFF.""" +def _save_tif(output_path, data, transform, crs, dtype='float32', count=1, nodata=None, nan_mask=None): + """Helper to save a 2D or 3D array as GeoTIFF. + + Args: + nan_mask: Optional boolean mask (True=NaN) to apply before saving. + Restores NaN zones in gradient-derived products that were + computed on the filled DEM. + """ + if nan_mask is not None: + data = np.array(data, dtype=dtype, copy=True) + data[nan_mask] = np.nan + # Auto-detect nodata for float types with NaN if nodata is None and dtype.startswith('float') and np.any(np.isnan(data)): nodata = float('nan') @@ -238,7 +248,9 @@ def generate_hillshade(dem_file, basename, vis_dir, resolution, shared=None): combined_hillshade = xp.mean(xp.array(hillshades), axis=0) slope_shaded = cos_slope combined = 0.7 * combined_hillshade + 0.3 * slope_shaded - _save_tif(output, to_cpu(combined), transform, crs) + + nan_mask = shared.nan_mask if shared else np.isnan(to_cpu(dem_np)) + _save_tif(output, to_cpu(combined), transform, crs, nan_mask=nan_mask) logger.info(f" ✓ Hillshade terminé ({time.time()-t0:.1f}s){gpu_tag}") return output except Exception as e: @@ -258,6 +270,7 @@ def generate_slope(dem_file, basename, vis_dir, resolution, shared=None): transform = shared.transform crs = shared.crs slope = shared.slope_deg + nan_mask = shared.nan_mask if HAS_GPU: slope = to_gpu(slope) else: @@ -265,7 +278,8 @@ def generate_slope(dem_file, basename, vis_dir, resolution, shared=None): dem = to_gpu(dem_np) dy, dx = xp.gradient(dem) slope = xp.arctan(xp.sqrt(dx**2 + dy**2)) * 180 / xp.pi - _save_tif(output, to_cpu(slope) if HAS_GPU else slope, transform, crs) + nan_mask = np.isnan(dem_np) + _save_tif(output, to_cpu(slope) if HAS_GPU else slope, transform, crs, nan_mask=nan_mask) logger.info(f" ✓ Pente terminée ({time.time()-t0:.1f}s){gpu_tag}") return output except Exception as e: @@ -285,6 +299,7 @@ def generate_aspect(dem_file, basename, vis_dir, resolution, shared=None): transform = shared.transform crs = shared.crs aspect = shared.aspect + nan_mask = shared.nan_mask if HAS_GPU: aspect = to_gpu(aspect) else: @@ -293,7 +308,8 @@ def generate_aspect(dem_file, basename, vis_dir, resolution, shared=None): dy, dx = xp.gradient(dem) aspect = xp.arctan2(dy, dx) * 180 / xp.pi aspect = xp.mod(aspect, 360) - _save_tif(output, to_cpu(aspect) if HAS_GPU else aspect, transform, crs) + nan_mask = np.isnan(dem_np) + _save_tif(output, to_cpu(aspect) if HAS_GPU else aspect, transform, crs, nan_mask=nan_mask) logger.info(f" ✓ Aspect terminé ({time.time()-t0:.1f}s){gpu_tag}") return output except Exception as e: @@ -314,6 +330,7 @@ def generate_curvature(dem_file, basename, vis_dir, resolution, shared=None): crs = shared.crs dx = shared.dx dy = shared.dy + nan_mask = shared.nan_mask if HAS_GPU: dx = to_gpu(dx) dy = to_gpu(dy) @@ -321,10 +338,11 @@ def generate_curvature(dem_file, basename, vis_dir, resolution, shared=None): dem_np, transform, crs = _read_dem(dem_file) dem = to_gpu(dem_np) dy, dx = xp.gradient(dem) + nan_mask = np.isnan(dem_np) d2z_dx2 = xp.gradient(dx, axis=1) d2z_dy2 = xp.gradient(dy, axis=0) curvature = (d2z_dx2 + d2z_dy2) / 2 - _save_tif(output, to_cpu(curvature), transform, crs) + _save_tif(output, to_cpu(curvature), transform, crs, nan_mask=nan_mask) logger.info(f" ✓ Courbure terminée ({time.time()-t0:.1f}s){gpu_tag}") return output except Exception as e: diff --git a/run.sh b/run.sh index dbd166f..68d0bbe 100755 --- a/run.sh +++ b/run.sh @@ -9,10 +9,10 @@ # -v Mode verbeux (timestamps + niveaux) # --debug Mode debug (détails internes fichier:ligne) # -f / --force Régénérer tous les fichiers même si existants -# --no-keep-tif Supprimer les TIFF intermédiaires (par défaut: conservés) +# --keep-tif Conserver les fichiers TIFF pour régénérer les WebP # --force-classification # Reclassifier le sol même si le fichier .las existe déjà -# --ground-classification {auto,smrf,pmf,csf} +# --ground-classification {auto,smrf,csf} # Méthode de classification du sol (défaut: auto) # --file NOM... Traiter un ou plusieurs fichiers LAZ spécifiques # --test Exécuter les tests unitaires @@ -40,8 +40,8 @@ if [ $# -eq 0 ]; then echo " -f / --force Régénérer tous les fichiers même si les WebP existent" echo " --force-classification" echo " Reclassifier le sol même si le fichier .las existe" - echo " --no-keep-tif Supprimer les TIFF intermédiaires (défaut: conservés)" - echo " --ground-classification {auto,smrf,pmf,csf}" + echo " --keep-tif Conserver les TIFF pour régénérer les WebP" + echo " --ground-classification {auto,smrf,csf}" echo " Méthode de classification du sol (défaut: auto)" echo " --file NOM... Traiter un ou plusieurs fichiers LAZ (nom complet sans .laz/.las)" echo " --test Exécuter les tests unitaires" @@ -52,9 +52,9 @@ if [ $# -eq 0 ]; then echo " $0 -g -w 4 # GPU + 4 workers" echo " $0 -g -v # GPU + verbeux" echo " $0 -g -r 0.2 # Haute résolution" - echo " $0 -g --force # Régénérer les WebP" + echo " $0 -g --force # Régénérer WebP (DTM conservé si --keep-tif)" echo " $0 -g --force-classification # Reclassifier le sol seulement" - echo " $0 -g --ground-classification pmf # Forcer PMF" + echo " $0 -g --ground-classification csf # Forcer CSF (terrain complexe)" echo " $0 -g --file LHD_...IGN69.copc # Un fichier" exit 0 fi @@ -67,7 +67,7 @@ FORCE_FLAG="" FILE_ARGS="" GROUND_METHOD="" FORCE_CLASSIFY_FLAG="" -NO_KEEP_TIF_FLAG="" +KEEP_TIF_FLAG="" # Parse arguments manually (more robust than getopts for mixed short/long options) while [ $# -gt 0 ]; do @@ -80,7 +80,7 @@ while [ $# -gt 0 ]; do --debug) VERBOSE_FLAG="--debug"; shift ;; --force) FORCE_FLAG="--force"; shift ;; --force-classification) FORCE_CLASSIFY_FLAG="--force-classification"; shift ;; - --no-keep-tif) NO_KEEP_TIF_FLAG="--no-keep-tif"; shift ;; + --keep-tif) KEEP_TIF_FLAG="--keep-tif"; shift ;; --ground-classification) GROUND_METHOD="$2"; shift 2 ;; --ground-classification=*) GROUND_METHOD="${1#--ground-classification=}"; shift ;; --file) shift; while [ $# -gt 0 ] && [[ ! "$1" =~ ^- ]]; do FILE_ARGS="$FILE_ARGS $1"; shift; done ;; @@ -99,8 +99,8 @@ while [ $# -gt 0 ]; do echo " -f / --force Régénérer tous les fichiers même si les WebP existent" echo " --force-classification" echo " Reclassifier le sol même si le fichier .las existe" - echo " --no-keep-tif Supprimer les TIFF intermédiaires (défaut: conservés)" - echo " --ground-classification {auto,smrf,pmf,csf}" + echo " --keep-tif Conserver les TIFF pour régénérer les WebP" + echo " --ground-classification {auto,smrf,csf}" echo " Méthode de classification du sol (défaut: auto)" echo " --file NOM... Traiter un ou plusieurs fichiers LAZ (nom complet sans .laz/.las)" echo " --test Exécuter les tests unitaires" @@ -111,9 +111,9 @@ while [ $# -gt 0 ]; do echo " $0 -g -w 4 # GPU + 4 workers" echo " $0 -g -v # GPU + verbeux" echo " $0 -g -r 0.2 # Haute résolution" - echo " $0 -g --force # Régénérer les WebP" + echo " $0 -g --force # Régénérer WebP (DTM conservé si --keep-tif)" echo " $0 -g --force-classification # Reclassifier le sol seulement" - echo " $0 -g --ground-classification pmf # Forcer PMF" + echo " $0 -g --ground-classification csf # Forcer CSF (terrain complexe)" echo " $0 -g --file LHD_...IGN69.copc # Un fichier" exit 0 ;; @@ -159,14 +159,14 @@ echo " GPU : $([ -n "$GPU_FLAG" ] && echo 'OUI' || echo 'non')" echo " Verbeux : $([ -n "$VERBOSE_FLAG" ] && echo 'OUI' || echo 'non')" echo " Force : $([ -n "$FORCE_FLAG" ] && echo 'OUI' || echo 'non')" echo " Force classif.: $([ -n "$FORCE_CLASSIFY_FLAG" ] && echo 'OUI' || echo 'non')" -echo " Keep TIFF : $([ -n "$NO_KEEP_TIF_FLAG" ] && echo 'non' || echo 'OUI')" +echo " Keep TIFF : $([ -n "$KEEP_TIF_FLAG" ] && echo 'OUI' || echo 'non')" echo " Classification sol : $([ -n "$GROUND_METHOD" ] && echo "$GROUND_METHOD" || echo 'auto')" if [ -n "$FILE_ARGS" ]; then echo " Fichiers :${FILE_ARGS}" fi echo "============================================" -CMD_ARGS="-o /data/output -r $RESOLUTION -w $WORKERS $VERBOSE_FLAG $FORCE_FLAG $FORCE_CLASSIFY_FLAG $NO_KEEP_TIF_FLAG" +CMD_ARGS="-o /data/output -r $RESOLUTION -w $WORKERS $VERBOSE_FLAG $FORCE_FLAG $FORCE_CLASSIFY_FLAG $KEEP_TIF_FLAG" if [ -n "$GROUND_METHOD" ]; then CMD_ARGS="$CMD_ARGS --ground-classification $GROUND_METHOD" fi