diff --git a/CLAUDE.md b/CLAUDE.md
index 937d6bb..2a58a9f 100644
--- a/CLAUDE.md
+++ b/CLAUDE.md
@@ -18,8 +18,6 @@ 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)
```
@@ -34,19 +32,27 @@ docker run --rm --gpus all -v $(pwd)/input:/data/input:ro -v $(pwd)/output:/data
### Module responsibilities
- **`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.
+- **`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`.
-- **`visualizations.py`** — 15 `generate_*` functions + 2 IGN overlay lambdas. All take `(dem_file, basename, vis_dir, resolution)` and return a TIF path or None.
+- **`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.
-- **`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`.
+- **`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.
+
+### SharedDEM optimization
+
+`SharedDEM` pre-computes once per file:
+- DEM data (single I/O read)
+- NaN mask + filled DEM (single `_fill_nans` call, avoiding ~20 redundant calls)
+- Gradient components (dy, dx, slope, aspect) shared by hillshade, slope, aspect, curvature
+- LRM at 15m kernel (shared by lrm + anomalies)
+
+`_filter_nanaware_from_filled()` applies filters on the pre-filled DEM, skipping the expensive `_fill_nans` interpolation.
### Adding a visualization
Three places must be updated:
-1. `visualizations.py` — add `generate_X(dem_file, basename, vis_dir, resolution)` function
+1. `visualizations.py` — add `generate_X(dem_file, basename, vis_dir, resolution, shared=None)` function
2. `pipeline.py` `VIZ_STEPS` — add `('name', generate_X)` entry
3. `rendering.py` `COLORMAPS` — add entry keyed by the output filename keyword
@@ -61,7 +67,11 @@ Override with `--ground-classification {auto,smrf,pmf,csf}`.
### NaN handling
-DTM small gaps (< 1m from existing data) are filled using `rasterio.fill.fillnodata`. Large gaps remain as NaN. Visualization functions use `_fill_nans()` and `_filter_nanaware()` to avoid NaN propagation through filters.
+DTM small gaps (< 1m from existing data) are filled using `rasterio.fill.fillnodata`. Large gaps remain as NaN. `SharedDEM` fills NaN once; `_filter_nanaware_from_filled()` applies filters on the pre-filled array and restores the NaN mask.
+
+### Flow accumulation
+
+Uses priority-flood algorithm (Wang & Liu 2006) for sink filling, which is O(n log n) instead of iterative minimum_filter. D8 accumulation uses numba JIT; falls back to pure Python if numba unavailable.
### Parallel processing
@@ -72,7 +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 (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.
+- **Output format**: Visualizations saved as WebP. TIFF intermediates deleted after conversion unless `--keep-tif`. No COGs or viewer — only WebP + PDF report remain.
+- **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 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}
-
-