diff --git a/lidar_pipeline/cli.py b/lidar_pipeline/cli.py index 5e24c16..9b1578a 100644 --- a/lidar_pipeline/cli.py +++ b/lidar_pipeline/cli.py @@ -23,6 +23,12 @@ def setup_logging(verbose=False, debug=False): verbose: If True, include timestamps and level names. debug: If True, set level to DEBUG and add file:line info. """ + # Ensure UTF-8 output for French messages (é, è, ê, etc.) + if hasattr(sys.stdout, 'reconfigure'): + sys.stdout.reconfigure(encoding='utf-8', errors='replace') + if hasattr(sys.stderr, 'reconfigure'): + sys.stderr.reconfigure(encoding='utf-8', errors='replace') + if debug: level = logging.DEBUG fmt = "%(asctime)s.%(msecs)03d %(levelname)-5s [%(filename)s:%(lineno)d] %(message)s" @@ -116,11 +122,6 @@ Exemples: action="store_true", help="Conserver les fichiers TIFF intermédiaires (sinon supprimés après conversion WebP)" ) - parser.add_argument( - "--no-viewer", - action="store_true", - help="Ne pas générer le viewer web (COGs + HTML MapLibre)" - ) parser.add_argument( "--ground-classification", choices=["auto", "smrf", "pmf", "csf"], @@ -176,7 +177,6 @@ Exemples: ground_method=args.ground_classification, force_classify=args.force_classification, keep_tif=args.keep_tif, - no_viewer=args.no_viewer ) # If --file is specified, process only matching files diff --git a/lidar_pipeline/dtm.py b/lidar_pipeline/dtm.py index d76bf27..3038b0d 100644 --- a/lidar_pipeline/dtm.py +++ b/lidar_pipeline/dtm.py @@ -250,7 +250,7 @@ def classify_ground(laz_file, temp_dir, method='auto', force=False): return None -def create_dtm_fast(las_file, basename, dtm_dir, resolution): +def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False): """Create DTM using fast binning method with gap filling. Args: @@ -258,10 +258,17 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution): basename: Base name for output file. dtm_dir: Directory for output DTM GeoTIFF. resolution: Grid resolution in meters per pixel. + force: If True, regenerate even if DTM already exists. Returns: Path to output DTM GeoTIFF, or None on failure. """ + output_tif = dtm_dir / f"{basename}_dtm.tif" + + if output_tif.exists() and not force: + logger.info(f" DTM déjà existant — fichier réutilisé: {output_tif.name}") + return output_tif + import laspy logger.info(" → Génération DTM...") diff --git a/lidar_pipeline/pipeline.py b/lidar_pipeline/pipeline.py index 7e486d8..9373a1f 100644 --- a/lidar_pipeline/pipeline.py +++ b/lidar_pipeline/pipeline.py @@ -12,6 +12,7 @@ import multiprocessing import shutil import time from concurrent.futures import ProcessPoolExecutor, as_completed +from datetime import datetime from pathlib import Path import subprocess @@ -55,6 +56,7 @@ _file_filter = FilePrefixFilter() from .dtm import classify_ground, create_dtm_fast from .visualizations import ( + SharedDEM, generate_hillshade, generate_slope, generate_aspect, generate_curvature, generate_lrm, generate_svf, generate_openness, generate_mslrm, generate_tpi, generate_sailore, @@ -63,8 +65,7 @@ from .visualizations import ( ) from .gpu import gpu_cleanup from .ign import generate_ign_overlay -from .rendering import tif_to_png, generate_pdf_report, convert_to_cog, generate_cog_metadata -from .viewer import generate_viewer +from .rendering import tif_to_png, generate_pdf_report # Ordered list of visualization steps. @@ -106,7 +107,7 @@ VIZ_STEPS = [ class LidarArchaeoPipeline: """Orchestrates the LiDAR archaeological analysis pipeline.""" - def __init__(self, input_dir, output_dir, resolution=0.5, workers=1, force=False, ground_method='auto', force_classify=False, keep_tif=False, no_viewer=False): + def __init__(self, input_dir, output_dir, resolution=0.5, workers=1, force=False, ground_method='auto', force_classify=False, keep_tif=False): self.input_dir = Path(input_dir) self.output_dir = Path(output_dir) self.resolution = resolution @@ -115,7 +116,6 @@ class LidarArchaeoPipeline: self.ground_method = ground_method self.force_classify = force_classify self.keep_tif = keep_tif - self.no_viewer = no_viewer self.temp_dir = self.output_dir / "temp" if not self.input_dir.exists(): @@ -140,7 +140,6 @@ class LidarArchaeoPipeline: logger.info(f" Classification sol : {self.ground_method}") logger.info(f" Force classif.: {'OUI' if self.force_classify else 'non'}") logger.info(f" Keep TIFF : {'OUI' if self.keep_tif else 'non'}") - logger.info(f" Viewer web : {'non' if self.no_viewer else 'OUI'}") def find_laz_files(self): """Find all LAZ/LAS files in input directory.""" @@ -173,6 +172,12 @@ class LidarArchaeoPipeline: file_vis_dir = self.vis_dir / basename file_vis_dir.mkdir(exist_ok=True) + # Pre-compute shared DEM data (gradient, NaN mask, LRM) once for all visualizations + logger.info(" Pré-calcul données partagées (gradient, LRM)...") + t_shared = time.time() + shared = SharedDEM(dtm_file, self.resolution) + logger.info(f" ✓ Données partagées prêtes ({time.time()-t_shared:.1f}s)") + vis_results = {} total = len(VIZ_STEPS) elapsed_times = [] @@ -216,7 +221,11 @@ class LidarArchaeoPipeline: logger.info(f" [{idx}/{total}] {name}...") t0 = time.time() try: - result = func(dtm_file, basename, file_vis_dir, self.resolution) + # IGN overlays don't use SharedDEM (they download external data) + if name in ('ortho', 'topo'): + result = func(dtm_file, basename, file_vis_dir, self.resolution) + else: + result = func(dtm_file, basename, file_vis_dir, self.resolution, shared=shared) vis_results[name] = result elapsed = time.time() - t0 elapsed_times.append(elapsed) @@ -237,20 +246,23 @@ class LidarArchaeoPipeline: gpu_cleanup() # Convert to WebP (only newly generated TIFs, not skipped ones) - # Also generate COGs for web viewer if enabled logger.info(" Conversion images WebP:") - cog_dir = file_vis_dir / "cog" - if not self.no_viewer: - cog_dir.mkdir(exist_ok=True) + source_info = { + 'method': self.ground_method, + 'date': datetime.now().strftime('%Y-%m-%d'), + 'basename': basename, + } for name, tif_file in vis_results.items(): if tif_file and isinstance(tif_file, Path) and tif_file.suffix == '.tif' and tif_file.exists(): - # Generate COG before WebP conversion (which may delete the TIF) - if not self.no_viewer: - convert_to_cog(tif_file, cog_dir) - webp_file = tif_to_png(tif_file, file_vis_dir, self.resolution, keep_tif=self.keep_tif or not self.no_viewer) + webp_file = tif_to_png(tif_file, file_vis_dir, self.resolution, keep_tif=self.keep_tif, source_info=source_info) if webp_file: logger.info(f" ✓ {webp_file.name}") + # Clean up remaining TIF files unless --keep-tif + if not self.keep_tif: + for tif in file_vis_dir.glob("*.tif"): + tif.unlink(missing_ok=True) + return vis_results def process_file(self, laz_file): @@ -263,60 +275,55 @@ class LidarArchaeoPipeline: logger.info(f"FICHIER : {basename}") logger.info("=" * 60) - # Step 1: Ground classification - logger.info("[1/6] Classification du sol...") - t1 = time.time() - las_file = classify_ground(laz_file, self.temp_dir, method=self.ground_method, force=self.force_classify) - t_classif = time.time() - t1 - if not las_file: - logger.error(f" ✗ Échec classification ({t_classif:.1f}s)") - return False - logger.info(f" ✓ Classification terminée ({t_classif:.1f}s)") + # Skip ground classification + DTM if DTM already exists (unless --force) + dtm_path = self.dtm_dir / f"{basename}_dtm.tif" + if dtm_path.exists() and not self.force: + logger.info("[1/5] Classification du sol — sautée (DTM existant)") + logger.info("[2/5] Génération DTM — sautée (DTM existant)") + dtm_file = dtm_path + t_classif = 0 + t_dtm = 0 + else: + # Step 1: Ground classification + logger.info("[1/5] Classification du sol...") + t1 = time.time() + las_file = classify_ground(laz_file, self.temp_dir, method=self.ground_method, force=self.force_classify) + t_classif = time.time() - t1 + if not las_file: + logger.error(f" ✗ Échec classification ({t_classif:.1f}s)") + return False + logger.info(f" ✓ Classification terminée ({t_classif:.1f}s)") - # Step 2: Generate DTM - logger.info("[2/6] Génération DTM...") - t2 = time.time() - dtm_file = create_dtm_fast(las_file, basename, self.dtm_dir, self.resolution) - t_dtm = time.time() - t2 - if not dtm_file: - logger.error(f" ✗ Échec DTM ({t_dtm:.1f}s)") - return False - logger.info(f" ✓ DTM terminé ({t_dtm:.1f}s)") + # Step 2: Generate DTM + logger.info("[2/5] Génération DTM...") + t2 = time.time() + dtm_file = create_dtm_fast(las_file, basename, self.dtm_dir, self.resolution, force=self.force) + t_dtm = time.time() - t2 + if not dtm_file: + logger.error(f" ✗ Échec DTM ({t_dtm:.1f}s)") + return False + logger.info(f" ✓ DTM terminé ({t_dtm:.1f}s)") # Step 3: Visualizations - logger.info("[3/6] Visualisations archéologiques...") + logger.info("[3/5] Visualisations archéologiques...") self.generate_all_visualizations(dtm_file, basename) # Step 4: PDF report file_vis_dir = self.vis_dir / basename - logger.info("[4/6] Rapport PDF A3...") - t4 = time.time() - generate_pdf_report(basename, file_vis_dir, self.pdf_dir, self.resolution) - t_pdf = time.time() - t4 - logger.info(f" ✓ Rapport PDF terminé ({t_pdf:.1f}s)") - - # Step 5: COGs for web viewer - logger.info("[5/6] Génération métadonnées viewer web...") - t5 = time.time() - if not self.no_viewer: - # Convert DTM to COG as well - dtm_cog_dir = self.dtm_dir / "cog" - dtm_cog_dir.mkdir(exist_ok=True) - for dtm_file in sorted(self.dtm_dir.glob(f"{basename}_dtm.tif")): - convert_to_cog(dtm_file, dtm_cog_dir) - generate_cog_metadata(self.vis_dir, basename) - t_cog = time.time() - t5 - logger.info(f" ✓ Métadonnées viewer web terminées ({t_cog:.1f}s)") - - # Step 6: Web viewer - if not self.no_viewer: - logger.info("[6/6] Génération viewer web...") - t6 = time.time() - generate_viewer(basename, file_vis_dir, self.vis_dir) - t_viewer = time.time() - t6 - logger.info(f" ✓ Viewer web terminé ({t_viewer:.1f}s)") + pdf_file = self.pdf_dir / f"{basename}_rapport.pdf" + if pdf_file.exists() and not self.force: + logger.info(f"[4/5] Rapport PDF déjà existant — ignoré: {pdf_file.name}") else: - logger.info("[6/6] Viewer web: ignoré (--no-viewer)") + logger.info("[4/5] Rapport PDF A3...") + t4 = time.time() + generate_pdf_report(basename, file_vis_dir, self.pdf_dir, self.resolution) + t_pdf = time.time() - t4 + logger.info(f" ✓ Rapport PDF terminé ({t_pdf:.1f}s)") + + # Step 5: Clean up DTM TIF unless --keep-tif + if not self.keep_tif: + for dtm_file in sorted(self.dtm_dir.glob(f"{basename}_dtm.tif")): + dtm_file.unlink(missing_ok=True) t_total = time.time() - t_start logger.info(f"✓ {basename} terminé en {t_total:.1f}s") @@ -349,7 +356,7 @@ class LidarArchaeoPipeline: logger.info(f"Fichiers: {len(files)}") with ProcessPoolExecutor(max_workers=self.workers) as executor: future_to_file = { - executor.submit(_process_file_standalone, str(laz_file), str(self.input_dir), str(self.output_dir), self.resolution, self.force, self.ground_method, self.force_classify, self.keep_tif, self.no_viewer): laz_file + executor.submit(_process_file_standalone, str(laz_file), str(self.input_dir), str(self.output_dir), self.resolution, self.force, self.ground_method, self.force_classify, self.keep_tif): laz_file for laz_file in files } done = 0 @@ -411,7 +418,7 @@ class LidarArchaeoPipeline: logger.warning(f" Note: Impossible de supprimer les fichiers temporaires: {e}") -def _process_file_standalone(laz_file_str, input_dir, output_dir, resolution, force=False, ground_method='auto', force_classify=False, keep_tif=False, no_viewer=False): +def _process_file_standalone(laz_file_str, input_dir, output_dir, resolution, force=False, ground_method='auto', force_classify=False, keep_tif=False): """Standalone function for multiprocessing — creates its own pipeline instance. Each worker gets its own temp directory to avoid file conflicts. @@ -419,6 +426,11 @@ def _process_file_standalone(laz_file_str, input_dir, output_dir, resolution, fo # Configure logging in worker process (spawn doesn't inherit parent config) import logging import sys + # Ensure UTF-8 output — spawn workers may default to ASCII + if hasattr(sys.stdout, 'reconfigure'): + sys.stdout.reconfigure(encoding='utf-8', errors='replace') + if hasattr(sys.stderr, 'reconfigure'): + sys.stderr.reconfigure(encoding='utf-8', errors='replace') worker_logger = logging.getLogger("lidar") if not worker_logger.handlers: handler = logging.StreamHandler(sys.stdout) @@ -427,7 +439,7 @@ def _process_file_standalone(laz_file_str, input_dir, output_dir, resolution, fo worker_logger.addHandler(handler) worker_logger.addFilter(_file_filter) - pipeline = LidarArchaeoPipeline(input_dir, output_dir, resolution=resolution, workers=1, force=force, ground_method=ground_method, force_classify=force_classify, keep_tif=keep_tif, no_viewer=no_viewer) + pipeline = LidarArchaeoPipeline(input_dir, output_dir, resolution=resolution, workers=1, force=force, ground_method=ground_method, force_classify=force_classify, keep_tif=keep_tif) basename = _file_basename(laz_file_str) pipeline.temp_dir = pipeline.output_dir / "temp" / basename pipeline.temp_dir.mkdir(exist_ok=True) diff --git a/lidar_pipeline/rendering.py b/lidar_pipeline/rendering.py index 800a5c3..d4bed61 100644 --- a/lidar_pipeline/rendering.py +++ b/lidar_pipeline/rendering.py @@ -8,6 +8,7 @@ Contains: import logging import time +from datetime import datetime from pathlib import Path import numpy as np @@ -238,7 +239,31 @@ def _apply_colormap(data, tif_file): return data, 'terrain', title, 'Altitude normalisée', '', False -def tif_to_png(tif_file, vis_dir, resolution, keep_tif=False): +def _nice_scale(extent_m): + """Choose a nice round scale distance that fits well in the image. + + Returns (scale_m, label) where label is like '100 m' or '500 m' or '1 km'. + """ + nice_scales = [50, 100, 200, 500, 1000, 2000, 5000, 10000] + # Pick the largest scale that is <= 15% of the extent + for s in nice_scales: + if s <= extent_m * 0.15: + continue + # s is the first one > 15% — take the one below + break + else: + s = nice_scales[-1] + # Actually pick the largest scale <= 20% of extent + chosen = nice_scales[0] + for s in nice_scales: + if s <= extent_m * 0.20: + chosen = s + if chosen >= 1000: + return chosen, f"{chosen // 1000} km" + return chosen, f"{chosen} m" + + +def tif_to_png(tif_file, vis_dir, resolution, keep_tif=False, source_info=None): """Convert GeoTIFF to visualization WebP with GPS coordinates, legend, and scale bar. Args: @@ -440,20 +465,32 @@ def tif_to_png(tif_file, vis_dir, resolution, keep_tif=False): spine.set_color('black') spine.set_linewidth(0.8) - # North arrow - north_ax = inset_axes(ax, width="4%", height="7%", loc='upper right', - bbox_to_anchor=(-0.03, 0.12, 1, 1), bbox_transform=ax.transAxes) - north_ax.set_xlim(0, 1) - north_ax.set_ylim(0, 1) + # North arrow — compass rose style + north_ax = inset_axes(ax, width="5%", height="9%", loc='upper right', + bbox_to_anchor=(-0.03, 0.08, 1, 1), bbox_transform=ax.transAxes) + north_ax.set_xlim(-1.2, 1.2) + north_ax.set_ylim(-0.5, 1.5) north_ax.axis('off') - north_ax.plot([0.5, 0.5], [0.1, 0.65], color='black', linewidth=2.5, zorder=10) - north_ax.add_patch(MplPolygon([[0.5, 0.2], [0.35, 0.4], [0.5, 0.65], [0.65, 0.4]], - closed=True, facecolor='black', edgecolor='black', zorder=9)) - north_ax.text(0.5, 0.95, 'N', ha='center', va='top', - fontsize=13, fontweight='bold', color='black', zorder=11) + north_ax.set_aspect('equal') + # N arrow + north_ax.annotate('N', xy=(0, 1.3), fontsize=11, fontweight='bold', + ha='center', va='bottom', color='#b22222') + north_ax.plot([0, 0], [0.0, 1.0], color='#b22222', linewidth=2.0, zorder=10) + north_ax.add_patch(MplPolygon([[0, 0.3], [-0.2, 0.7], [0, 1.0], [0.2, 0.7]], + closed=True, facecolor='#b22222', edgecolor='#b22222', zorder=9)) + # Cardinal ticks + for angle, label in [(90, ''), (0, 'E'), (180, 'O'), (270, 'S')]: + rad = np.radians(angle) + r_text = 1.25 + north_ax.plot([0.85*np.cos(rad), 1.05*np.cos(rad)], + [0.85*np.sin(rad), 1.05*np.sin(rad)], + color='#555555', linewidth=0.8, zorder=5) + if label: + north_ax.text(r_text*np.cos(rad), r_text*np.sin(rad), label, + fontsize=7, ha='center', va='center', color='#555555') - # Bottom info bar — fixed position - info_ax = fig.add_axes([data_left, 0.02, data_width_frac + cbar_width + 0.02, 0.07]) + # Bottom info bar — enriched with source, method, date + info_ax = fig.add_axes([data_left, 0.015, data_width_frac + cbar_width + 0.02, 0.09]) info_ax.axis('off') extent_km_x = (max_x - min_x) / 1000 @@ -465,42 +502,73 @@ def tif_to_png(tif_file, vis_dir, resolution, keep_tif=False): alt_min = float(np.nanmin(valid_data)) if len(valid_data) > 0 else 0 alt_max = float(np.nanmax(valid_data)) if len(valid_data) > 0 else 0 + # Build info lines + line1_parts = [] if gps_coords: nw_lat, nw_lon = gps_coords['NW'] se_lat, se_lon = gps_coords['SE'] - info_text = ( - f"GPS: {nw_lat:.5f}N {nw_lon:.5f}E - {se_lat:.5f}N {se_lon:.5f}E | " - f"EPSG:2154 | Res: {resolution}m/px | " - f"Emprise: {extent_km_x:.1f}x{extent_km_y:.1f}km" - ) + line1_parts.append(f"GPS: {nw_lat:.5f}°N {nw_lon:.5f}°E — {se_lat:.5f}°N {se_lon:.5f}°E") else: - info_text = ( - f"EPSG:2154 | X: {min_x:.0f}-{max_x:.0f} Y: {min_y:.0f}-{max_y:.0f} | " - f"Res: {resolution}m/px | Emprise: {extent_km_x:.1f}x{extent_km_y:.1f}km" - ) + line1_parts.append(f"X: {min_x:.0f}–{max_x:.0f} Y: {min_y:.0f}–{max_y:.0f}") + line1_parts.append(f"EPSG:2154") + line1_parts.append(f"Res: {resolution}m/px") + line1_parts.append(f"Emprise: {extent_km_x:.1f}×{extent_km_y:.1f}km") + if not is_rgb: + line1_parts.append(f"Alt: {alt_min:.1f}–{alt_max:.1f}m") - info_ax.text(0.01, 0.5, info_text, - transform=info_ax.transAxes, fontsize=8.5, + line2_parts = [] + line2_parts.append("Source: LiDAR HD IGN") + if source_info: + if source_info.get('method'): + line2_parts.append(f"Classif.: {source_info['method'].upper()}") + if source_info.get('date'): + line2_parts.append(f"Date: {source_info['date']}") + else: + line2_parts.append(datetime.now().strftime("Date: %Y-%m-%d")) + + info_text_line1 = " | ".join(line1_parts) + info_text_line2 = " | ".join(line2_parts) + + info_ax.text(0.01, 0.7, info_text_line1, + transform=info_ax.transAxes, fontsize=8, verticalalignment='center', family='monospace', - bbox=dict(boxstyle='round,pad=0.3', facecolor='#f0f0f0', + bbox=dict(boxstyle='round,pad=0.2', facecolor='#f0f0f0', edgecolor='#aaaaaa', alpha=0.95)) + info_ax.text(0.01, 0.2, info_text_line2, + transform=info_ax.transAxes, fontsize=7.5, + verticalalignment='center', family='monospace', + color='#444444', + bbox=dict(boxstyle='round,pad=0.2', facecolor='#f8f8f8', + edgecolor='#cccccc', alpha=0.9)) - # Scale bar - scale_m = 100 + # Scale bar — adaptive with alternating black/white segments + extent_m_x = max_x - min_x + scale_m, scale_label = _nice_scale(extent_m_x) pixels_per_meter = 1.0 / pixel_size_x scale_px = int(scale_m * pixels_per_meter) - scale_bar_frac = scale_px / width - bar_left = 0.80 - bar_bottom = 0.15 - bar_width_frac = min(scale_bar_frac, 0.15) - bar_height = 0.35 + n_segments = 5 + segment_px = scale_px / n_segments + bar_bottom_y = 0.55 + bar_top_y = 0.85 + bar_height = bar_top_y - bar_bottom_y + scale_start_x = 0.80 - info_ax.add_patch(RectPatch((bar_left, bar_bottom), bar_width_frac, bar_height, - facecolor='black', edgecolor='black', linewidth=0.5, - transform=info_ax.transAxes, clip_on=False)) - info_ax.text(bar_left + bar_width_frac / 2, bar_bottom + bar_height + 0.12, - f"{scale_m} m", ha='center', va='bottom', fontsize=9, fontweight='bold', - transform=info_ax.transAxes) + for seg_i in range(n_segments): + color = 'black' if seg_i % 2 == 0 else 'white' + seg_left = scale_start_x + seg_i * segment_px / width + seg_width_frac = segment_px / width + info_ax.add_patch(RectPatch((seg_left, bar_bottom_y), seg_width_frac, bar_height, + facecolor=color, edgecolor='black', linewidth=0.5, + transform=info_ax.transAxes, clip_on=False)) + info_ax.text(scale_start_x + scale_px / (2 * width), bar_top_y + 0.12, + f"{scale_label}", ha='center', va='bottom', fontsize=8, fontweight='bold', + transform=info_ax.transAxes) + # Scale end ticks + info_ax.plot([scale_start_x, scale_start_x], [bar_bottom_y - 0.05, bar_top_y + 0.05], + color='black', linewidth=1, transform=info_ax.transAxes, clip_on=False) + info_ax.plot([scale_start_x + scale_px / width, scale_start_x + scale_px / width], + [bar_bottom_y - 0.05, bar_top_y + 0.05], + color='black', linewidth=1, transform=info_ax.transAxes, clip_on=False) fig.patch.set_facecolor('white') @@ -527,145 +595,6 @@ def tif_to_png(tif_file, vis_dir, resolution, keep_tif=False): return None -def convert_to_cog(tif_file, cog_dir): - """Convert a GeoTIFF to Cloud Optimized GeoTIFF for web map serving. - - Args: - tif_file: Path to input GeoTIFF. - cog_dir: Directory for output COG file. - - Returns: - Path to output COG file, or None on failure. - """ - cog_file = cog_dir / tif_file.name - - if cog_file.exists(): - logger.debug(f" COG déjà existant: {cog_file.name}") - return cog_file - - try: - from rio_cogeo.cogeo import cog_translate - from rio_cogeo.profiles import cog_profiles - - with rasterio.open(tif_file) as src: - is_rgb = src.count >= 3 - - # Use deflate compression profile - profile = dict(cog_profiles["deflate"]) # Make a mutable copy - - if not is_rgb: - # Single-band terrain data: keep float32 for continuous values - profile.update(dtype="float32") - - cog_translate(str(tif_file), str(cog_file), profile, quiet=True) - logger.debug(f" COG créé: {cog_file.name}") - return cog_file - - except ImportError: - logger.warning(" rio-cogeo non disponible — COG non généré") - return None - except Exception as e: - logger.warning(f" Erreur conversion COG: {e}") - return None - - -def generate_cog_metadata(vis_dir, basename): - """Generate metadata JSON for all visualization layers. - - Scans the COG directory for files and reads their geographic bounds. - Creates a metadata.json with WGS84 bounds and layer info for the web viewer. - - Args: - vis_dir: Directory containing COG files (vis_dir/basename/cog/). - basename: Base name for the LiDAR tile. - - Returns: - Path to metadata.json, or None on failure. - """ - import json - - cog_dir = vis_dir / basename / "cog" - if not cog_dir.exists(): - return None - - # Find the DTM to get geographic bounds - # The COG files inherit bounds from the DTM, so we can read any COG - cog_files = sorted(cog_dir.glob("*.tif")) - if not cog_files: - return None - - # Read bounds from first COG file - try: - ref_cog = cog_files[0] - with rasterio.open(ref_cog) as src: - bounds = src.bounds - crs = src.crs - if crs and HAS_WARP: - l93_xs = [bounds.left, bounds.right] - l93_ys = [bounds.top, bounds.bottom] - lons, lats = warp_transform(crs, 'EPSG:4326', l93_xs, l93_ys) - bounds_wgs84 = { - 'west': float(min(lons)), - 'south': float(min(lats)), - 'east': float(max(lons)), - 'north': float(max(lats)), - } - else: - # Fallback: use Lambert 93 bounds directly - bounds_wgs84 = { - 'west': float(bounds.left), - 'south': float(bounds.bottom), - 'east': float(bounds.right), - 'north': float(bounds.top), - } - - bounds_l93 = { - 'min_x': float(bounds.left), - 'min_y': float(bounds.bottom), - 'max_x': float(bounds.right), - 'max_y': float(bounds.top), - } - resolution = float(abs(src.transform.a)) - except Exception as e: - logger.warning(f" Erreur lecture bounds COG: {e}") - return None - - # Build layer list with titles - layers = [] - for cog_file in cog_files: - name = cog_file.stem.replace(basename + "_", "") - # Find matching title from COLORMAPS or RGB_LEGENDS - title = name.replace("_", " ").title() - for key in sorted(COLORMAPS.keys(), key=len, reverse=True): - if key in name: - title = COLORMAPS[key]['title'] - break - for key in RGB_LEGENDS: - if key in name: - title = RGB_LEGENDS[key]['title'] - break - - layers.append({ - 'name': name, - 'title': title, - 'file': cog_file.name, - }) - - metadata = { - 'basename': basename, - 'bounds_l93': bounds_l93, - 'bounds_wgs84': bounds_wgs84, - 'resolution': resolution, - 'layers': layers, - } - - metadata_file = vis_dir / basename / "metadata.json" - with open(metadata_file, 'w') as f: - json.dump(metadata, f, indent=2, ensure_ascii=False) - - logger.info(f" Métadonnées: {len(layers)} couches, bounds={bounds_wgs84}") - return metadata_file - def generate_pdf_report(basename, vis_dir, pdf_dir, resolution): """Generate A3 PDF report for a LiDAR file with all visualizations. diff --git a/lidar_pipeline/server.py b/lidar_pipeline/server.py deleted file mode 100644 index 4520697..0000000 --- a/lidar_pipeline/server.py +++ /dev/null @@ -1,117 +0,0 @@ -"""TiTiler-based web server for serving LiDAR COG visualizations. - -Provides: -- COG tile serving via TiTiler API -- Static file serving for viewer HTML -- CORS headers for local development - -Usage: - python -m lidar_pipeline.server /path/to/output [--port 8000] -""" - -import logging -import sys -from pathlib import Path - -logger = logging.getLogger("lidar") - - -def create_app(output_dir): - """Create the FastAPI application with TiTiler and static file serving. - - Args: - output_dir: Path to the output directory containing visualisations/ and viewer/. - - Returns: - FastAPI application instance. - """ - from fastapi import FastAPI - from fastapi.responses import HTMLResponse, FileResponse, JSONResponse - from fastapi.middleware.cors import CORSMiddleware - from titiler.core.factory import TilerFactory - - output_dir = Path(output_dir).resolve() - - # TiTiler COG endpoint - cog = TilerFactory(router_prefix="/cog") - - app = FastAPI(title="LiDAR Archéologique", docs_url="/docs") - - # CORS for local development - app.add_middleware(CORSMiddleware, allow_origins=["*"], allow_methods=["*"], allow_headers=["*"]) - - # Mount TiTiler routes - app.include_router(cog.router, prefix="/cog") - - @app.get("/") - async def index(): - """Landing page with links.""" - return HTMLResponse( - 'LiDAR Server' - '' - '' - '

Serveur LiDAR Archéologique

' - '

Visualisations

' - '

TiTiler API: /cog/

' - '' - ) - - @app.get("/viewer/{basename}") - @app.get("/viewer") - async def serve_viewer(basename: str = ""): - """Serve viewer HTML files.""" - if not basename: - # List available viewers - viewer_dir = output_dir / 'visualisations' / 'viewer' - if viewer_dir.exists(): - viewers = sorted(viewer_dir.glob('*.html')) - if viewers: - links = ''.join( - f'
  • {f.stem}
    • ' - for f in viewers - ) - return HTMLResponse( - f'LiDAR Viewer' - f'' - f'

    Visualisations LiDAR

    ' - ) - return HTMLResponse('

    Aucun viewer disponible

    ') - - viewer_file = output_dir / 'visualisations' / 'viewer' / f"{basename}.html" - if viewer_file.exists(): - return FileResponse(str(viewer_file), media_type='text/html') - return JSONResponse({'error': f'Viewer not found: {basename}'}, status_code=404) - - return app - - -def main(): - """Entry point for the LiDAR web server.""" - import argparse - import uvicorn - - parser = argparse.ArgumentParser(description='Serveur cartographique LiDAR') - parser.add_argument('output_dir', help='Répertoire de sortie contenant les visualisations') - parser.add_argument('--host', default='0.0.0.0', help='Hôte (défaut: 0.0.0.0)') - parser.add_argument('--port', type=int, default=8000, help='Port (défaut: 8000)') - args = parser.parse_args() - - output_dir = Path(args.output_dir) - if not output_dir.exists(): - logger.error(f"Répertoire introuvable: {output_dir}") - sys.exit(1) - - app = create_app(output_dir) - - print(f"Serveur LiDAR Archéologique") - print(f" Viewer: http://{args.host}:{args.port}/viewer") - print(f" TiTiler: http://{args.host}:{args.port}/cog/") - print(f" Répertoire: {output_dir}") - - uvicorn.run(app, host=args.host, port=args.port, log_level='info') - - -if __name__ == '__main__': - main() \ No newline at end of file diff --git a/lidar_pipeline/viewer.py b/lidar_pipeline/viewer.py deleted file mode 100644 index 100a3ca..0000000 --- a/lidar_pipeline/viewer.py +++ /dev/null @@ -1,416 +0,0 @@ -"""Web map viewer generator for LiDAR visualizations. - -Generates a self-contained HTML file with MapLibre GL JS that displays -all visualization layers with opacity controls and IGN/OSM basemaps. -Served by TiTiler for COG tile access. -""" - -import json -import logging -from pathlib import Path - -logger = logging.getLogger("lidar") - -# Layer ordering for the viewer panel (archaeological priority) -LAYER_ORDER = [ - 'hillshade_multi', - 'slope', - 'aspect', - 'curvature', - 'svf', - 'lrm', - 'positive_openness', - 'negative_openness', - 'mslrm', - 'tpi', - 'sailore', - 'roughness', - 'anomalies', - 'wavelet', - 'flow', - 'ortho', - 'topo', -] - -# Colormaps for TiTiler rendering of single-band COGs -LAYER_COLORMAPS = { - 'hillshade_multi': 'gray', - 'slope': 'inferno', - 'aspect': 'turbo', - 'curvature': 'RdYlBu_r', - 'svf': 'viridis', - 'lrm': 'RdBu_r', - 'positive_openness': 'YlOrBr', - 'negative_openness': 'GnBu_r', - 'mslrm': 'RdBu_r', - 'tpi': 'BrBG', - 'sailore': 'seismic', - 'roughness': 'magma', - 'anomalies': 'coolwarm', - 'wavelet': 'cividis', - 'flow': 'Blues', -} - - -def generate_viewer(basename, vis_dir, output_vis_dir, titiler_url='http://localhost:8000'): - """Generate a MapLibre GL JS viewer HTML file for the LiDAR tile. - - Args: - basename: Base name for the LiDAR tile. - vis_dir: Per-file visualization directory (vis_dir/basename/). - output_vis_dir: Parent visualization directory for viewer output. - titiler_url: Base URL of the TiTiler server. - - Returns: - Path to the generated HTML file, or None on failure. - """ - # Read metadata.json - metadata_file = vis_dir / "metadata.json" - if not metadata_file.exists(): - logger.error(f" Métadonnées manquantes: {metadata_file}") - return None - - with open(metadata_file) as f: - metadata = json.load(f) - - bounds_wgs84 = metadata['bounds_wgs84'] - resolution = metadata.get('resolution', 0.5) - - # Determine which files are RGB (ortho/topo) - layers = [] - for layer in metadata['layers']: - is_rgb = 'ortho' in layer['name'] or 'topo' in layer['name'] - layers.append({ - 'name': layer['name'], - 'title': layer['title'], - 'file': layer['file'], - 'is_rgb': is_rgb, - }) - - # Sort layers by archaeological priority - def layer_sort_key(l): - name = l['name'] - for i, key in enumerate(LAYER_ORDER): - if key in name: - return i - return len(LAYER_ORDER) - - layers.sort(key=layer_sort_key) - - # COG path prefix for TiTiler (absolute path inside Docker container) - cog_path_prefix = f'/data/output/visualisations/{basename}/cog' - - # Build layer controls HTML - layer_controls = [] - for layer in layers: - checked = 'checked' if layer['name'] == 'hillshade_multi' else '' - initial_opacity = 100 if layer['name'] == 'hillshade_multi' else 0 - layer_type = 'RGB' if layer['is_rgb'] else 'DEM' - layer_controls.append( - f'
    \n' - f' \n' - f' \n' - f'
    ' - ) - layer_controls_html = '\n'.join(layer_controls) - - # Serialize data for JS - bounds_json = json.dumps(bounds_wgs84) - layers_json = json.dumps(layers, ensure_ascii=False) - colormaps_json = json.dumps(LAYER_COLORMAPS) - center_lng = (bounds_wgs84['west'] + bounds_wgs84['east']) / 2 - center_lat = (bounds_wgs84['south'] + bounds_wgs84['north']) / 2 - - # Build the full HTML using string concatenation to avoid f-string issues with {z}/{x}/{y} - html_parts = [] - html_parts.append(f''' - - - - -LiDAR Archéologique — {basename} - - - - - - -
    -

    LiDAR Archéologique

    -
    - {basename}
    - Res: {resolution}m/px | EPSG:2154 -
    - -
    -

    Fond de carte

    - - - -
    - -

    Visualisations

    -
    -{layer_controls_html} -
    -
    - -
    -
    - - - -''') - - html = ''.join(html_parts) - - # Write viewer HTML - viewer_dir = output_vis_dir / 'viewer' - viewer_dir.mkdir(exist_ok=True) - viewer_file = viewer_dir / f"{basename}.html" - - with open(viewer_file, 'w', encoding='utf-8') as f: - f.write(html) - - logger.info(f" Viewer: {viewer_file}") - return viewer_file \ No newline at end of file diff --git a/lidar_pipeline/visualizations.py b/lidar_pipeline/visualizations.py index e6cbc25..bda618e 100644 --- a/lidar_pipeline/visualizations.py +++ b/lidar_pipeline/visualizations.py @@ -2,6 +2,10 @@ Each function takes (dem_file, basename, vis_dir, resolution) as explicit parameters and returns the path to the output GeoTIFF file, or None on error. + +When a SharedDEM object is provided via the `shared` parameter, pre-computed +data (gradient, NaN mask, LRM) is reused across visualizations to avoid +redundant I/O and computation. """ import logging @@ -26,6 +30,76 @@ else: xp = np +class SharedDEM: + """Pre-computed DEM data shared across all visualizations. + + Reads the DEM once and pre-computes: + - NaN mask and filled DEM (avoids 20+ calls to _fill_nans) + - Gradient components (shared by hillshade, slope, aspect, curvature) + - LRM at 15m kernel (shared by lrm + anomalies) + """ + + def __init__(self, dem_file, resolution): + dem_np, transform, crs = _read_dem(dem_file) + self.dem_file = dem_file + self.resolution = resolution + self.transform = transform + self.crs = crs + self.nan_mask = np.isnan(dem_np) + self.dem_np = dem_np.astype(np.float32) + + # Pre-fill NaNs once (saves ~20 calls to NearestNDInterpolator) + self.filled, _ = _fill_nans(self.dem_np) + + # Initialize GPU lazy caches before any filter calls + self._filled_gpu = None + self._dem_gpu = None + + # Pre-compute gradient (shared by hillshade, slope, aspect, curvature) + self.dy = np.gradient(self.filled, resolution, axis=0) + self.dx = np.gradient(self.filled, resolution, axis=1) + self.slope_rad = np.arctan(np.sqrt(self.dx**2 + self.dy**2)) + self.slope_deg = np.degrees(self.slope_rad) + self.aspect = np.mod(np.degrees(np.arctan2(self.dy, self.dx)), 360) + + # Pre-compute LRM at 15m (shared by lrm + anomalies) + sigma_15 = 15.0 / resolution + local_mean_15 = _filter_nanaware_from_filled(self, xp_gaussian_filter, sigma=sigma_15) + self.lrm_15 = self.dem_np - local_mean_15 + self.lrm_15[self.nan_mask] = np.nan + + @property + def filled_gpu(self): + """Lazy GPU copy of the filled DEM.""" + if self._filled_gpu is None and HAS_GPU: + self._filled_gpu = to_gpu(self.filled) + return self._filled_gpu + + @property + def dem_gpu(self): + """Lazy GPU copy of the DEM.""" + if self._dem_gpu is None and HAS_GPU: + self._dem_gpu = to_gpu(self.dem_np) + return self._dem_gpu + + +def _filter_nanaware_from_filled(shared, filter_func, *args, **kwargs): + """Apply filter on pre-filled DEM data (skips expensive _fill_nans). + + Uses the SharedDEM.filled array directly, then restores NaN mask. + If GPU is available, uses the lazy GPU copy to avoid CPU↔GPU transfers. + """ + if HAS_GPU and shared.filled_gpu is not None: + filled_gpu = to_gpu(shared.filled) + result_gpu = filter_func(filled_gpu, *args, **kwargs) + result = to_cpu(result_gpu) + gpu_cleanup() + else: + result = filter_func(shared.filled, *args, **kwargs) + result[shared.nan_mask] = np.nan + 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.""" # Auto-detect nodata for float types with NaN @@ -38,7 +112,7 @@ def _save_tif(output_path, data, transform, crs, dtype='float32', count=1, nodat output_path, 'w', driver='GTiff', height=height, width=width, count=count, dtype=dtype, crs=crs, transform=transform, - compress='lzw', nodata=nodata + compress='deflate', nodata=nodata ) as dst: dst.write(data.astype(dtype), 1) elif data.ndim == 3: @@ -47,7 +121,7 @@ def _save_tif(output_path, data, transform, crs, dtype='float32', count=1, nodat output_path, 'w', driver='GTiff', height=height, width=width, count=bands, dtype=dtype, crs=crs, transform=transform, - compress='lzw', nodata=nodata + compress='deflate', nodata=nodata ) as dst: for i in range(bands): dst.write(data[i].astype(dtype), i + 1) @@ -116,7 +190,7 @@ def _filter_nanaware(arr, filter_func, *args, use_gpu=True, **kwargs): # Core terrain visualizations # ============================================================ -def generate_hillshade(dem_file, basename, vis_dir, resolution): +def generate_hillshade(dem_file, basename, vis_dir, resolution, shared=None): """Generate multi-directional hillshade with slope shading — GPU if available. Combines 4-direction hillshade (NW, NE, SW, SE) with slope shading @@ -129,19 +203,29 @@ def generate_hillshade(dem_file, basename, vis_dir, resolution): output = vis_dir / f"{basename}_hillshade_multi.tif" try: - dem_np, transform, crs = _read_dem(dem_file) - dem = to_gpu(dem_np) - dy, dx = xp.gradient(dem) + if shared: + transform = shared.transform + crs = shared.crs + dem = to_gpu(shared.dem_np) + dy = to_gpu(shared.dy) if HAS_GPU else shared.dy + dx = to_gpu(shared.dx) if HAS_GPU else shared.dx + slope = to_gpu(shared.slope_rad) if HAS_GPU else shared.slope_rad + aspect = xp.arctan2(dy, dx) + sin_slope = xp.sin(slope) + cos_slope = xp.cos(slope) + else: + dem_np, transform, crs = _read_dem(dem_file) + dem = to_gpu(dem_np) + dy, dx = xp.gradient(dem) + slope = xp.arctan(xp.sqrt(dx**2 + dy**2)) + aspect = xp.arctan2(dy, dx) + sin_slope = xp.sin(slope) + cos_slope = xp.cos(slope) azimuts = [315, 45, 225, 135] altitude = 30 hillshades = [] - slope = xp.arctan(xp.sqrt(dx**2 + dy**2)) - aspect = xp.arctan2(dy, dx) - sin_slope = xp.sin(slope) - cos_slope = xp.cos(slope) - alt_rad = xp.radians(xp.array(altitude)) sin_alt = xp.sin(alt_rad) cos_alt = xp.cos(alt_rad) @@ -152,8 +236,7 @@ def generate_hillshade(dem_file, basename, vis_dir, resolution): hillshades.append(xp.clip(hs, 0, 1)) combined_hillshade = xp.mean(xp.array(hillshades), axis=0) - # Blend with slope shading for better micro-relief on flat terrain - slope_shaded = cos_slope # bright on flat, dark on steep + slope_shaded = cos_slope combined = 0.7 * combined_hillshade + 0.3 * slope_shaded _save_tif(output, to_cpu(combined), transform, crs) logger.info(f" ✓ Hillshade terminé ({time.time()-t0:.1f}s){gpu_tag}") @@ -163,7 +246,7 @@ def generate_hillshade(dem_file, basename, vis_dir, resolution): return None -def generate_slope(dem_file, basename, vis_dir, resolution): +def generate_slope(dem_file, basename, vis_dir, resolution, shared=None): """Generate slope map (degrees) — GPU if available.""" gpu_tag = " [GPU]" if HAS_GPU else "" logger.info(f" → Pente (Slope){gpu_tag}...") @@ -171,11 +254,18 @@ def generate_slope(dem_file, basename, vis_dir, resolution): output = vis_dir / f"{basename}_slope.tif" try: - dem_np, transform, crs = _read_dem(dem_file) - 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), transform, crs) + if shared: + transform = shared.transform + crs = shared.crs + slope = shared.slope_deg + if HAS_GPU: + slope = to_gpu(slope) + else: + dem_np, transform, crs = _read_dem(dem_file) + 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) logger.info(f" ✓ Pente terminée ({time.time()-t0:.1f}s){gpu_tag}") return output except Exception as e: @@ -183,7 +273,7 @@ def generate_slope(dem_file, basename, vis_dir, resolution): return None -def generate_aspect(dem_file, basename, vis_dir, resolution): +def generate_aspect(dem_file, basename, vis_dir, resolution, shared=None): """Generate aspect (slope orientation) map — GPU if available.""" gpu_tag = " [GPU]" if HAS_GPU else "" logger.info(f" → Aspect (Orientation){gpu_tag}...") @@ -191,12 +281,19 @@ def generate_aspect(dem_file, basename, vis_dir, resolution): output = vis_dir / f"{basename}_aspect.tif" try: - dem_np, transform, crs = _read_dem(dem_file) - dem = to_gpu(dem_np) - dy, dx = xp.gradient(dem) - aspect = xp.arctan2(dy, dx) * 180 / xp.pi - aspect = xp.mod(aspect, 360) - _save_tif(output, to_cpu(aspect), transform, crs) + if shared: + transform = shared.transform + crs = shared.crs + aspect = shared.aspect + if HAS_GPU: + aspect = to_gpu(aspect) + else: + dem_np, transform, crs = _read_dem(dem_file) + dem = to_gpu(dem_np) + 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) logger.info(f" ✓ Aspect terminé ({time.time()-t0:.1f}s){gpu_tag}") return output except Exception as e: @@ -204,7 +301,7 @@ def generate_aspect(dem_file, basename, vis_dir, resolution): return None -def generate_curvature(dem_file, basename, vis_dir, resolution): +def generate_curvature(dem_file, basename, vis_dir, resolution, shared=None): """Generate curvature (terrain concavity/convexity) map — GPU if available.""" gpu_tag = " [GPU]" if HAS_GPU else "" logger.info(f" → Courbure (Curvature){gpu_tag}...") @@ -212,12 +309,20 @@ def generate_curvature(dem_file, basename, vis_dir, resolution): output = vis_dir / f"{basename}_curvature.tif" try: - dem_np, transform, crs = _read_dem(dem_file) - dem = to_gpu(dem_np) - dz_dx = xp.gradient(dem, axis=1) - dz_dy = xp.gradient(dem, axis=0) - d2z_dx2 = xp.gradient(dz_dx, axis=1) - d2z_dy2 = xp.gradient(dz_dy, axis=0) + if shared: + transform = shared.transform + crs = shared.crs + dx = shared.dx + dy = shared.dy + if HAS_GPU: + dx = to_gpu(dx) + dy = to_gpu(dy) + else: + dem_np, transform, crs = _read_dem(dem_file) + dem = to_gpu(dem_np) + dy, dx = xp.gradient(dem) + 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) logger.info(f" ✓ Courbure terminée ({time.time()-t0:.1f}s){gpu_tag}") @@ -232,7 +337,7 @@ def generate_curvature(dem_file, basename, vis_dir, resolution): # GPU-accelerated visualizations # ============================================================ -def generate_lrm(dem_file, basename, vis_dir, resolution): +def generate_lrm(dem_file, basename, vis_dir, resolution, shared=None): """Local Relief Model - deviation from local mean (GPU if available).""" gpu_tag = " [GPU]" if HAS_GPU else "" logger.info(f" → Local Relief Model{gpu_tag}...") @@ -240,11 +345,16 @@ def generate_lrm(dem_file, basename, vis_dir, resolution): output = vis_dir / f"{basename}_lrm.tif" try: - dem_np, transform, crs = _read_dem(dem_file) - nan_mask = np.isnan(dem_np) - local_mean = _filter_nanaware(dem_np, xp_gaussian_filter, sigma=15/resolution) - lrm = dem_np - local_mean - lrm[nan_mask] = np.nan + if shared: + transform = shared.transform + crs = shared.crs + lrm = shared.lrm_15.copy() + else: + dem_np, transform, crs = _read_dem(dem_file) + nan_mask = np.isnan(dem_np) + local_mean = _filter_nanaware(dem_np, xp_gaussian_filter, sigma=15/resolution) + lrm = dem_np - local_mean + lrm[nan_mask] = np.nan _save_tif(output, lrm.astype(np.float32), transform, crs) logger.info(f" ✓ LRM terminé ({time.time()-t0:.1f}s){gpu_tag}") return output @@ -253,7 +363,7 @@ def generate_lrm(dem_file, basename, vis_dir, resolution): return None -def generate_svf(dem_file, basename, vis_dir, resolution): +def generate_svf(dem_file, basename, vis_dir, resolution, shared=None): """Sky-View Factor - ray-tracing on 16 azimuths (GPU if available). For each pixel, trace rays in N directions, find the max horizon @@ -266,23 +376,33 @@ def generate_svf(dem_file, basename, vis_dir, resolution): output = vis_dir / f"{basename}_svf.tif" try: - dem_np, transform, crs = _read_dem(dem_file) - rows, cols = dem_np.shape - res = resolution - - dem = to_gpu(dem_np) + if shared: + transform = shared.transform + crs = shared.crs + dem_np = shared.dem_np + rows, cols = dem_np.shape + res = resolution + dem = to_gpu(shared.filled) if HAS_GPU else shared.filled + nan_mask = shared.nan_mask + else: + dem_np, transform, crs = _read_dem(dem_file) + rows, cols = dem_np.shape + res = resolution + nan_mask = np.isnan(dem_np) + filled, _ = _fill_nans(dem_np) + dem = to_gpu(filled) if HAS_GPU else filled n_dirs = 16 angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False) - dx = np.cos(angles) - dy = np.sin(angles) + dx_dir = np.cos(angles) + dy_dir = np.sin(angles) max_dist = int(50 / res) padded = xp.pad(dem, max_dist, mode='constant', constant_values=xp.nan) svf = xp.zeros_like(dem) for d_idx in range(n_dirs): - ddx, ddy = dx[d_idx], dy[d_idx] + ddx, ddy = dx_dir[d_idx], dy_dir[d_idx] horizon = xp.zeros_like(dem) # Pre-compute all valid steps for this direction @@ -307,6 +427,7 @@ def generate_svf(dem_file, basename, vis_dir, resolution): svf /= n_dirs svf_np = to_cpu(svf).astype(np.float32) + svf_np[nan_mask] = np.nan _save_tif(output, svf_np, transform, crs) logger.info(f" ✓ SVF terminé ({time.time()-t0:.1f}s){gpu_tag}") return output @@ -315,7 +436,7 @@ def generate_svf(dem_file, basename, vis_dir, resolution): return None -def generate_openness(dem_file, basename, vis_dir, resolution, positive=True): +def generate_openness(dem_file, basename, vis_dir, resolution, positive=True, shared=None): """Positive/Negative Openness - true zenith/nadir angle computation (GPU if available). For each pixel, in 8 directions (N, NE, E, SE, S, SW, W, NW): @@ -330,23 +451,33 @@ def generate_openness(dem_file, basename, vis_dir, resolution, positive=True): output = vis_dir / f"{basename}_{name}.tif" try: - dem_np, transform, crs = _read_dem(dem_file) - rows, cols = dem_np.shape - res = resolution - - dem = to_gpu(dem_np) + if shared: + transform = shared.transform + crs = shared.crs + dem_np = shared.dem_np + rows, cols = dem_np.shape + res = resolution + dem = to_gpu(shared.filled) if HAS_GPU else shared.filled + nan_mask = shared.nan_mask + else: + dem_np, transform, crs = _read_dem(dem_file) + rows, cols = dem_np.shape + res = resolution + nan_mask = np.isnan(dem_np) + filled, _ = _fill_nans(dem_np) + dem = to_gpu(filled) if HAS_GPU else filled n_dirs = 8 angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False) - dx = np.cos(angles) - dy = np.sin(angles) + dx_dir = np.cos(angles) + dy_dir = np.sin(angles) max_dist = int(50 / res) padded = xp.pad(dem, max_dist, mode='constant', constant_values=xp.nan) openness_sum = xp.zeros_like(dem) for d_idx in range(n_dirs): - ddx, ddy = dx[d_idx], dy[d_idx] + ddx, ddy = dx_dir[d_idx], dy_dir[d_idx] max_angle = xp.zeros_like(dem) for step in range(1, max_dist + 1): @@ -370,6 +501,7 @@ def generate_openness(dem_file, basename, vis_dir, resolution, positive=True): openness_sum += max_angle openness_result = to_cpu(xp.degrees(openness_sum / n_dirs)).astype(np.float32) + openness_result[nan_mask] = np.nan _save_tif(output, openness_result, transform, crs) logger.info(f" ✓ {name} terminé ({time.time()-t0:.1f}s){gpu_tag}") return output @@ -378,7 +510,7 @@ def generate_openness(dem_file, basename, vis_dir, resolution, positive=True): return None -def generate_mslrm(dem_file, basename, vis_dir, resolution): +def generate_mslrm(dem_file, basename, vis_dir, resolution, shared=None): """Multi-Scale Relief Model (MSRM) - LRM at 5 scales combined (GPU if available).""" gpu_tag = " [GPU]" if HAS_GPU else "" logger.info(f" → Multi-Scale Relief Model (MSRM){gpu_tag}...") @@ -386,15 +518,24 @@ def generate_mslrm(dem_file, basename, vis_dir, resolution): output = vis_dir / f"{basename}_mslrm.tif" try: - dem_np, transform, crs = _read_dem(dem_file) - nan_mask = np.isnan(dem_np) + if shared: + transform = shared.transform + crs = shared.crs + dem_np = shared.dem_np + nan_mask = shared.nan_mask + else: + dem_np, transform, crs = _read_dem(dem_file) + nan_mask = np.isnan(dem_np) sigmas = [5, 10, 25, 50, 100] lrm_stack = [] for sigma in sigmas: sigma_px = sigma / resolution - local_mean = _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_px) + if shared: + local_mean = _filter_nanaware_from_filled(shared, xp_gaussian_filter, sigma=sigma_px) + else: + local_mean = _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_px) lrm = dem_np - local_mean lrm[nan_mask] = np.nan valid_lrm = lrm[~nan_mask] @@ -416,7 +557,7 @@ def generate_mslrm(dem_file, basename, vis_dir, resolution): return None -def generate_tpi(dem_file, basename, vis_dir, resolution): +def generate_tpi(dem_file, basename, vis_dir, resolution, shared=None): """Multi-Scale Topographic Position Index (GPU if available). TPI = elevation - mean(neighborhood). @@ -428,20 +569,32 @@ def generate_tpi(dem_file, basename, vis_dir, resolution): output = vis_dir / f"{basename}_tpi.tif" try: - dem_np, transform, crs = _read_dem(dem_file) - nan_mask = np.isnan(dem_np) + if shared: + transform = shared.transform + crs = shared.crs + dem_np = shared.dem_np + nan_mask = shared.nan_mask + else: + dem_np, transform, crs = _read_dem(dem_file) + nan_mask = np.isnan(dem_np) fine_size = int(5 / resolution) if fine_size % 2 == 0: fine_size += 1 - fine_mean = _filter_nanaware(dem_np, xp_uniform_filter, size=fine_size) + if shared: + fine_mean = _filter_nanaware_from_filled(shared, xp_uniform_filter, size=fine_size) + else: + fine_mean = _filter_nanaware(dem_np, xp_uniform_filter, size=fine_size) tpi_fine = dem_np - fine_mean tpi_fine[nan_mask] = np.nan broad_size = int(100 / resolution) if broad_size % 2 == 0: broad_size += 1 - broad_mean = _filter_nanaware(dem_np, xp_uniform_filter, size=broad_size) + if shared: + broad_mean = _filter_nanaware_from_filled(shared, xp_uniform_filter, size=broad_size) + else: + broad_mean = _filter_nanaware(dem_np, xp_uniform_filter, size=broad_size) tpi_broad = dem_np - broad_mean tpi_broad[nan_mask] = np.nan @@ -464,7 +617,7 @@ def generate_tpi(dem_file, basename, vis_dir, resolution): # SAILORE # ============================================================ -def generate_sailore(dem_file, basename, vis_dir, resolution): +def generate_sailore(dem_file, basename, vis_dir, resolution, shared=None): """SAILORE - Self-Adaptive Improved Local Relief Model (GPU if available). Kernel size adapts to local slope: flat areas get larger kernels, @@ -476,23 +629,40 @@ def generate_sailore(dem_file, basename, vis_dir, resolution): output = vis_dir / f"{basename}_sailore.tif" try: - dem_np, transform, crs = _read_dem(dem_file) - nan_mask = np.isnan(dem_np) - - gy, gx = np.gradient(dem_np, resolution) - slope = np.arctan(np.sqrt(gx**2 + gy**2)) - slope_deg = np.degrees(slope) - slope_deg[nan_mask] = np.nan + if shared: + transform = shared.transform + crs = shared.crs + dem_np = shared.dem_np + nan_mask = shared.nan_mask + slope_deg = shared.slope_deg + else: + dem_np, transform, crs = _read_dem(dem_file) + nan_mask = np.isnan(dem_np) + gy, gx = np.gradient(dem_np, resolution) + slope = np.arctan(np.sqrt(gx**2 + gy**2)) + slope_deg = np.degrees(slope) + slope_deg[nan_mask] = np.nan sigma_min = 2.0 / resolution sigma_max = 25.0 / resolution slope_norm = np.clip(slope_deg / 30.0, 0, 1) - lrm_fine = dem_np - _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_min) + if shared: + lrm_fine = dem_np - _filter_nanaware_from_filled(shared, xp_gaussian_filter, sigma=sigma_min) + else: + lrm_fine = dem_np - _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_min) lrm_fine[nan_mask] = np.nan - lrm_medium = dem_np - _filter_nanaware(dem_np, xp_gaussian_filter, sigma=(sigma_min + sigma_max) / 2) + + if shared: + lrm_medium = dem_np - _filter_nanaware_from_filled(shared, xp_gaussian_filter, sigma=(sigma_min + sigma_max) / 2) + else: + lrm_medium = dem_np - _filter_nanaware(dem_np, xp_gaussian_filter, sigma=(sigma_min + sigma_max) / 2) lrm_medium[nan_mask] = np.nan - lrm_coarse = dem_np - _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_max) + + if shared: + lrm_coarse = dem_np - _filter_nanaware_from_filled(shared, xp_gaussian_filter, sigma=sigma_max) + else: + lrm_coarse = dem_np - _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_max) lrm_coarse[nan_mask] = np.nan w_fine = slope_norm @@ -516,7 +686,7 @@ def generate_sailore(dem_file, basename, vis_dir, resolution): # Roughness # ============================================================ -def generate_roughness(dem_file, basename, vis_dir, resolution): +def generate_roughness(dem_file, basename, vis_dir, resolution, shared=None): """Surface roughness - standard deviation of elevation in a window (GPU-accelerated).""" gpu_tag = " [GPU]" if HAS_GPU else "" logger.info(f" → Rugosité de surface{gpu_tag}...") @@ -524,15 +694,28 @@ def generate_roughness(dem_file, basename, vis_dir, resolution): output = vis_dir / f"{basename}_roughness.tif" try: - dem_np, transform, crs = _read_dem(dem_file) - nan_mask = np.isnan(dem_np) + if shared: + transform = shared.transform + crs = shared.crs + dem_np = shared.dem_np + nan_mask = shared.nan_mask + else: + dem_np, transform, crs = _read_dem(dem_file) + nan_mask = np.isnan(dem_np) + window_size = int(5 / resolution) if window_size % 2 == 0: window_size += 1 # Vectorized std: sqrt(E[X²] - (E[X])²) via uniform_filter (NaN-aware) - local_mean = _filter_nanaware(dem_np.astype(np.float64), xp_uniform_filter, size=window_size) - local_mean_sq = _filter_nanaware(dem_np.astype(np.float64)**2, xp_uniform_filter, size=window_size) + if shared: + local_mean = _filter_nanaware_from_filled(shared, xp_uniform_filter, size=window_size) + # For local_mean_sq, we need to filter filled², not filled + local_mean_sq = _filter_nanaware(shared.filled.astype(np.float64)**2, xp_uniform_filter, size=window_size) + local_mean_sq[shared.nan_mask] = np.nan + else: + local_mean = _filter_nanaware(dem_np.astype(np.float64), xp_uniform_filter, size=window_size) + local_mean_sq = _filter_nanaware(dem_np.astype(np.float64)**2, xp_uniform_filter, size=window_size) roughness = np.sqrt(np.maximum(local_mean_sq - local_mean * local_mean, 0)) roughness[nan_mask] = np.nan @@ -549,7 +732,7 @@ def generate_roughness(dem_file, basename, vis_dir, resolution): # Anomalies # ============================================================ -def generate_anomalies(dem_file, basename, vis_dir, resolution): +def generate_anomalies(dem_file, basename, vis_dir, resolution, shared=None): """Statistical anomaly detection - z-score of local relief + Local Moran's I — GPU if available.""" gpu_tag = " [GPU]" if HAS_GPU else "" logger.info(f" → Détection anomalies statistiques{gpu_tag}...") @@ -557,12 +740,19 @@ def generate_anomalies(dem_file, basename, vis_dir, resolution): output = vis_dir / f"{basename}_anomalies.tif" try: - dem_np, transform, crs = _read_dem(dem_file) - nan_mask = np.isnan(dem_np) + if shared: + transform = shared.transform + crs = shared.crs + dem_np = shared.dem_np + nan_mask = shared.nan_mask + lrm = shared.lrm_15.copy() + else: + dem_np, transform, crs = _read_dem(dem_file) + nan_mask = np.isnan(dem_np) + lrm_mean_val = _filter_nanaware(dem_np, xp_gaussian_filter, sigma=15 / resolution) + lrm = dem_np - lrm_mean_val + lrm[nan_mask] = np.nan - lrm_mean_val = _filter_nanaware(dem_np, xp_gaussian_filter, sigma=15 / resolution) - lrm = dem_np - lrm_mean_val - lrm[nan_mask] = np.nan valid_lrm = lrm[~nan_mask] lrm_mean = np.nanmean(valid_lrm) if len(valid_lrm) > 0 else 0 lrm_std = max(np.nanstd(valid_lrm), 0.01) if len(valid_lrm) > 0 else 0.01 @@ -572,7 +762,10 @@ def generate_anomalies(dem_file, basename, vis_dir, resolution): if window % 2 == 0: window += 1 - local_mean = _filter_nanaware(z_score, xp_uniform_filter, size=window) + if shared: + local_mean = _filter_nanaware_from_filled(shared, xp_uniform_filter, size=window) + else: + local_mean = _filter_nanaware(z_score, xp_uniform_filter, size=window) z_mean = np.nanmean(valid_lrm) if len(valid_lrm) > 0 else 0 z_std = max(np.nanstd(z_score[~nan_mask]), 0.01) if np.any(~nan_mask) else 0.01 morans_i = z_score * (local_mean - z_mean) / z_std @@ -591,7 +784,7 @@ def generate_anomalies(dem_file, basename, vis_dir, resolution): # Wavelet # ============================================================ -def generate_wavelet(dem_file, basename, vis_dir, resolution): +def generate_wavelet(dem_file, basename, vis_dir, resolution, shared=None): """Mexican Hat wavelet multi-scale analysis (GPU if available). CWT 2D at multiple scales to detect circular features. @@ -602,15 +795,22 @@ def generate_wavelet(dem_file, basename, vis_dir, resolution): output = vis_dir / f"{basename}_wavelet.tif" try: - dem_np, transform, crs = _read_dem(dem_file) - nan_mask = np.isnan(dem_np) + if shared: + transform = shared.transform + crs = shared.crs + dem_np = shared.dem_np + nan_mask = shared.nan_mask + filled = shared.filled.astype(np.float64) + else: + dem_np, transform, crs = _read_dem(dem_file) + nan_mask = np.isnan(dem_np) + filled, _ = _fill_nans(dem_np.astype(np.float64)) scales = [2, 5, 10, 20, 50] wavelet_stack = [] for scale_m in scales: sigma_px = scale_m / resolution - filled, _ = _fill_nans(dem_np.astype(np.float64)) if HAS_GPU: from cupyx.scipy.ndimage import gaussian_laplace as gpu_gaussian_laplace response = -gpu_gaussian_laplace(to_gpu(filled), sigma=sigma_px) @@ -697,36 +897,68 @@ def _d8_accumulate_numba(flow_dir, nodata_mask, rows, cols): return None -def generate_flow(dem_file, basename, vis_dir, resolution): - """Flow accumulation using D8 algorithm — sink filling on GPU, accumulation via numba.""" +def _priority_flood(dem, nodata_mask): + """Priority-flood algorithm for sink filling (Wang & Liu 2006). + + O(n log n) compared to 50 iterations of minimum_filter. + Fills pits so water can flow downhill. + """ + import heapq + + rows, cols = dem.shape + filled = dem.copy() + closed = nodata_mask.copy() + open_queue = [] + + # Initialize border cells + for r in range(rows): + for c in [0, cols - 1]: + if not closed[r, c]: + heapq.heappush(open_queue, (filled[r, c], r, c)) + closed[r, c] = True + for c in range(1, cols - 1): + for r in [0, rows - 1]: + if not closed[r, c]: + heapq.heappush(open_queue, (filled[r, c], r, c)) + closed[r, c] = True + + dx8 = [1, 1, 0, -1, -1, -1, 0, 1] + dy8 = [0, 1, 1, 1, 0, -1, -1, -1] + + while open_queue: + elev, r, c = heapq.heappop(open_queue) + for d in range(8): + nr, nc = r + dy8[d], c + dx8[d] + if 0 <= nr < rows and 0 <= nc < cols and not closed[nr, nc]: + if filled[nr, nc] < elev: + filled[nr, nc] = elev # Fill the pit + closed[nr, nc] = True + heapq.heappush(open_queue, (filled[nr, nc], nr, nc)) + + return filled + + +def generate_flow(dem_file, basename, vis_dir, resolution, shared=None): + """Flow accumulation using D8 algorithm — priority-flood sink filling, accumulation via numba.""" gpu_tag = " [GPU]" if HAS_GPU else "" logger.info(f" → Accumulation de flux D8{gpu_tag}...") t0 = time.time() output = vis_dir / f"{basename}_flow.tif" try: - dem_np, transform, crs = _read_dem(dem_file) + if shared: + transform = shared.transform + crs = shared.crs + dem_np = shared.dem_np + nodata_mask = shared.nan_mask + else: + dem_np, transform, crs = _read_dem(dem_file) + nodata_mask = np.isnan(dem_np) + rows, cols = dem_np.shape - nodata_mask = np.isnan(dem_np) - # Sink filling — GPU-accelerated - dem_gpu = to_gpu(dem_np) - nodata_mask_gpu = xp.isnan(dem_gpu) - dem_filled = xp.copy(dem_gpu) - dem_filled[nodata_mask_gpu] = xp.nanmax(dem_gpu) + 1000 - - from scipy.ndimage import generate_binary_structure - struct = generate_binary_structure(2, 2) - - for _ in range(50): - neighbor_min = xp_minimum_filter(dem_filled, footprint=struct) - sinks = (dem_filled < neighbor_min) & ~nodata_mask_gpu - if not xp.any(sinks): - break - dem_filled = xp.where(sinks, neighbor_min, dem_filled) - - dem_filled[nodata_mask_gpu] = xp.nan - dem_filled_np = to_cpu(dem_filled) + # Sink filling — priority-flood (O(n log n), faster than 50× minimum_filter) + dem_filled_np = _priority_flood(dem_np, nodata_mask) # D8 slope — vectorized dx8 = np.array([1, 1, 0, -1, -1, -1, 0, 1], dtype=np.int32) diff --git a/run.sh b/run.sh index 10b071f..1609c2d 100755 --- a/run.sh +++ b/run.sh @@ -10,14 +10,12 @@ # --debug Mode debug (détails internes fichier:ligne) # -f / --force Régénérer tous les fichiers même si existants # --keep-tif Conserver les fichiers TIFF intermédiaires -# --no-viewer Ne pas générer le viewer web (COGs + HTML) # --force-classification # Reclassifier le sol même si le fichier .las existe déjà # --ground-classification {auto,smrf,pmf,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 -# serve Démarrer le serveur cartographique # -h Afficher l'aide complète set -e @@ -27,33 +25,11 @@ INPUT_DIR="${SCRIPT_DIR}/input" OUTPUT_DIR="${SCRIPT_DIR}/output" IMAGE_NAME="lidar-lidar" -# Serve command — start the web map server -if [ "$1" = "serve" ]; then - # Build l'image si elle n'existe pas - if ! docker image inspect "$IMAGE_NAME" >/dev/null 2>&1; then - echo "Build de l'image Docker..." - docker build -t "$IMAGE_NAME" "$SCRIPT_DIR" - fi - mkdir -p "$OUTPUT_DIR" - echo "============================================" - echo " Serveur cartographique LiDAR" - echo "============================================" - echo " Viewer: http://localhost:8000/viewer" - echo " TiTiler: http://localhost:8000/cog/" - echo "============================================" - docker run --rm -p 8000:8000 \ - -v "${OUTPUT_DIR}:/data/output" \ - "$IMAGE_NAME" \ - python3 -m lidar_pipeline.server /data/output - exit 0 -fi - # Afficher l'aide si aucun argument if [ $# -eq 0 ]; then echo "Pipeline LiDAR Archéologique" echo "" echo "Usage: $0 [options]" - echo " $0 serve # Démarrer le serveur cartographique" echo "" echo "Options:" echo " -r RESOLUTION Résolution en m/px (défaut: 0.5)" @@ -65,12 +41,10 @@ if [ $# -eq 0 ]; then echo " --force-classification" echo " Reclassifier le sol même si le fichier .las existe" echo " --keep-tif Conserver les fichiers TIFF intermédiaires" - echo " --no-viewer Ne pas générer le viewer web" echo " --ground-classification {auto,smrf,pmf,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" - echo " serve Démarrer le serveur cartographique" echo " -h Afficher cette aide" echo "" echo "Exemples:" @@ -82,7 +56,6 @@ if [ $# -eq 0 ]; then echo " $0 -g --force-classification # Reclassifier le sol seulement" echo " $0 -g --ground-classification pmf # Forcer PMF" echo " $0 -g --file LHD_...IGN69.copc # Un fichier" - echo " $0 serve # Démarrer le serveur web" exit 0 fi @@ -95,7 +68,6 @@ FILE_ARGS="" GROUND_METHOD="" FORCE_CLASSIFY_FLAG="" KEEP_TIF_FLAG="" -NO_VIEWER_FLAG="" # Parse arguments manually (more robust than getopts for mixed short/long options) while [ $# -gt 0 ]; do @@ -109,7 +81,6 @@ while [ $# -gt 0 ]; do --force) FORCE_FLAG="--force"; shift ;; --force-classification) FORCE_CLASSIFY_FLAG="--force-classification"; shift ;; --keep-tif) KEEP_TIF_FLAG="--keep-tif"; shift ;; - --no-viewer) NO_VIEWER_FLAG="--no-viewer"; 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 ;; @@ -118,7 +89,6 @@ while [ $# -gt 0 ]; do echo "Pipeline LiDAR Archéologique" echo "" echo "Usage: $0 [options]" - echo " $0 serve # Démarrer le serveur cartographique" echo "" echo "Options:" echo " -r RESOLUTION Résolution en m/px (défaut: 0.5)" @@ -130,12 +100,10 @@ while [ $# -gt 0 ]; do echo " --force-classification" echo " Reclassifier le sol même si le fichier .las existe" echo " --keep-tif Conserver les fichiers TIFF intermédiaires" - echo " --no-viewer Ne pas générer le viewer web" echo " --ground-classification {auto,smrf,pmf,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" - echo " serve Démarrer le serveur cartographique" echo " -h Afficher cette aide" echo "" echo "Exemples:" @@ -147,7 +115,6 @@ while [ $# -gt 0 ]; do echo " $0 -g --force-classification # Reclassifier le sol seulement" echo " $0 -g --ground-classification pmf # Forcer PMF" echo " $0 -g --file LHD_...IGN69.copc # Un fichier" - echo " $0 serve # Démarrer le serveur web" exit 0 ;; *) echo "Option invalide: $1" >&2; exit 1 ;; @@ -193,14 +160,13 @@ 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 "$KEEP_TIF_FLAG" ] && echo 'OUI' || echo 'non')" -echo " Viewer web : $([ -n "$NO_VIEWER_FLAG" ] && echo 'non' || echo 'OUI')" 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 $KEEP_TIF_FLAG $NO_VIEWER_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