diff --git a/CLAUDE.md b/CLAUDE.md index 084680e..937d6bb 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -18,6 +18,8 @@ All commands run inside Docker. Use `./run.sh` as the primary interface. ./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 +./run.sh -g --no-viewer # Skip web viewer generation +./run.sh serve # Start web map server ./run.sh # Print help (no args) ``` @@ -37,7 +39,9 @@ docker run --rm --gpus all -v $(pwd)/input:/data/input:ro -v $(pwd)/output:/data - **`visualizations.py`** — 15 `generate_*` functions + 2 IGN overlay lambdas. All take `(dem_file, basename, vis_dir, resolution)` and return a TIF path or None. - **`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. -- **`rendering.py`** — `COLORMAPS` dict maps filename keywords to (cmap, title, legend, description). `tif_to_png()` converts TIF→WebP with legend/scale/north arrow. `generate_pdf_report()` creates A3 PDF. +- **`rendering.py`** — `COLORMAPS` dict maps filename keywords to (cmap, title, legend, description). `tif_to_png()` converts TIF→WebP with legend/scale/north arrow. `convert_to_cog()` converts TIF→Cloud Optimized GeoTIFF. `generate_cog_metadata()` creates metadata JSON for web viewer. `generate_pdf_report()` creates A3 PDF. +- **`viewer.py`** — Generates MapLibre GL JS HTML viewer with layer controls, opacity sliders, and IGN/OSM basemaps. +- **`server.py`** — TiTiler-based Starlette server for serving COG tiles and viewer HTML. Entry point via `python -m lidar_pipeline.server`. ### Adding a visualization @@ -68,6 +72,7 @@ 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 (not PNG). TIFF intermediates deleted unless `--keep-tif`. PDF reports use `PILImage.open().convert('RGB')`. +- **Output format**: Visualizations saved as WebP (not PNG). TIFF intermediates deleted unless `--keep-tif` or viewer enabled. COGs generated for web viewer by default. PDF reports use `PILImage.open().convert('RGB')`. +- **Web viewer**: MapLibre GL JS + TiTiler. COGs served as raster tiles. `./run.sh serve` starts server on port 8000. - **Flow accumulation**: Uses numba JIT for D8 accumulation loop. Falls back to pure Python if numba unavailable. - **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/Dockerfile b/Dockerfile index 91af8f8..19f598a 100644 --- a/Dockerfile +++ b/Dockerfile @@ -32,7 +32,12 @@ RUN pip3 install --no-cache-dir \ tqdm \ Pillow \ pytest \ - numba + numba \ + rio-cogeo \ + titiler.core \ + fastapi \ + uvicorn \ + piexif # Install CuPy for GPU acceleration (optional - will fallback to numpy if not available) RUN pip3 install --no-cache-dir cupy-cuda12x || echo "CuPy not available - GPU acceleration disabled" diff --git a/lidar_pipeline/cli.py b/lidar_pipeline/cli.py index 92b5aaa..5e24c16 100644 --- a/lidar_pipeline/cli.py +++ b/lidar_pipeline/cli.py @@ -116,6 +116,11 @@ 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"], @@ -170,7 +175,8 @@ Exemples: force=args.force, ground_method=args.ground_classification, force_classify=args.force_classification, - keep_tif=args.keep_tif + keep_tif=args.keep_tif, + no_viewer=args.no_viewer ) # If --file is specified, process only matching files diff --git a/lidar_pipeline/pipeline.py b/lidar_pipeline/pipeline.py index 6a35387..7e486d8 100644 --- a/lidar_pipeline/pipeline.py +++ b/lidar_pipeline/pipeline.py @@ -63,7 +63,8 @@ from .visualizations import ( ) from .gpu import gpu_cleanup from .ign import generate_ign_overlay -from .rendering import tif_to_png, generate_pdf_report +from .rendering import tif_to_png, generate_pdf_report, convert_to_cog, generate_cog_metadata +from .viewer import generate_viewer # Ordered list of visualization steps. @@ -105,7 +106,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): + 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): self.input_dir = Path(input_dir) self.output_dir = Path(output_dir) self.resolution = resolution @@ -114,6 +115,7 @@ 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(): @@ -138,6 +140,7 @@ 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.""" @@ -234,10 +237,17 @@ 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) 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(): - webp_file = tif_to_png(tif_file, file_vis_dir, self.resolution, keep_tif=self.keep_tif) + # 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) if webp_file: logger.info(f" ✓ {webp_file.name}") @@ -254,7 +264,7 @@ class LidarArchaeoPipeline: logger.info("=" * 60) # Step 1: Ground classification - logger.info("[1/4] Classification du sol...") + 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 @@ -264,7 +274,7 @@ class LidarArchaeoPipeline: logger.info(f" ✓ Classification terminée ({t_classif:.1f}s)") # Step 2: Generate DTM - logger.info("[2/4] Génération 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 @@ -274,17 +284,40 @@ class LidarArchaeoPipeline: logger.info(f" ✓ DTM terminé ({t_dtm:.1f}s)") # Step 3: Visualizations - logger.info("[3/4] Visualisations archéologiques...") + logger.info("[3/6] Visualisations archéologiques...") self.generate_all_visualizations(dtm_file, basename) # Step 4: PDF report file_vis_dir = self.vis_dir / basename - logger.info("[4/4] Rapport PDF A3...") + 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)") + else: + logger.info("[6/6] Viewer web: ignoré (--no-viewer)") + 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") @@ -316,7 +349,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): 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, self.no_viewer): laz_file for laz_file in files } done = 0 @@ -378,7 +411,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): +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): """Standalone function for multiprocessing — creates its own pipeline instance. Each worker gets its own temp directory to avoid file conflicts. @@ -394,7 +427,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) + 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) 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 41aae74..800a5c3 100644 --- a/lidar_pipeline/rendering.py +++ b/lidar_pipeline/rendering.py @@ -527,6 +527,146 @@ 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 new file mode 100644 index 0000000..4520697 --- /dev/null +++ b/lidar_pipeline/server.py @@ -0,0 +1,117 @@ +"""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( + '
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'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 new file mode 100644 index 0000000..100a3ca --- /dev/null +++ b/lidar_pipeline/viewer.py @@ -0,0 +1,416 @@ +"""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'