19 Commits
ok1 ... main

Author SHA1 Message Date
d334892880 Improve visualizations: adaptive scales, revert z-score to std normalization
- MSRM/TPI/roughness/anomalies: revert z-score (x-mean)/std to std normalization x/std
  to preserve contrast and visibility of linear features (paths, ditches, trenches)
- MSRM: adaptive scales based on resolution, archaeological weight combination
- TPI: extend from 2 to 4 scales (3m/15m/50m/200m) with weighted combination
- Hillshade: 8 directions instead of 4, altitude 35° instead of 30°
- LRM: adaptive sigma based on resolution
- Openness: doubled radius (100m instead of 50m)
- Roughness: multi-scale (3m fine + 15m broad) instead of single 5x5 window
- Anomalies: uses MSRM multi-scale relief instead of single LRM 15m
- Wavelet: 8 adaptive scales, std normalization, archaeological weights
- Remove svf (Sky-View Factor) and local_dominance visualizations
- Add AVIF format support (default), quality 98
- Add multi-resolution support (-r 0.5,0.2)
- Improve Ctrl+C handling for immediate process termination
- Update rendering.py descriptions for all modified visualizations

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
2026-05-14 23:12:08 +02:00
ac56ba8084 Add multi-resolution support and remove PDF generation
- --resolution now accepts comma-separated values (e.g. 0.5,0.2)
- Additional resolutions get suffixed output dirs: basename_r0p2/
- DTM files are named basename_dtm_r0p2.tif for extra resolutions
- Ground classification is done once and shared across resolutions
- PDF report generation removed per user request
- Fix --file argument to accept full filenames with extensions

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
2026-05-14 21:29:45 +02:00
53b6369a1b Fix --file argument to accept full filenames with extensions
Previously --file LHD_FXX_...copc.laz would fail because it appended
extensions. Now tries exact filename match first, then falls back to
adding extensions.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
2026-05-14 21:17:03 +02:00
34b79ac2c2 Add WebP quality control and selective visualization (--only / --skip)
- WebP output now uses quality=85 by default (down from lossless),
  reducing file size by ~75% (35MB → 5-8MB per visualization)
- Added --quality N (1-100) and --lossless flags in CLI and run.sh
- Added --only and --skip to select/exclude specific visualizations
  (e.g., --only hillshade,svf,lrm or --skip ortho,topo)
- VIZ_STEPS filtering is done in LidarArchaeoPipeline.__init__
- SharedDEM is skipped when all selected visualizations already exist
- Invalid visualization names are validated at startup with clear error

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
2026-05-14 21:15:21 +02:00
6ed4972afc Handle empty point clouds gracefully at every pipeline stage
Fixes "zero-size array to reduction operation" crash on corrupt/incomplete
LAZ files. Added checks at each step:

- validate_laz(): check point_count > 0 via laspy header, parse PDAL
  info JSON for point count when using PDAL fallback
- detect_ground_method(): return 'smrf' default if point cloud is empty
  after PDAL conversion instead of crashing on np.max(empty_array)
- _read_with_pdal(): log warning and return None if converted file has
  0 points
- create_dtm_fast(): fail gracefully if ground file has 0 points
- classify_ground(): check output file size after PDAL pipeline to
  catch empty ground classifications early

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
2026-05-14 20:56:52 +02:00
ab9d694dd4 List all processed LAZ files with status in pipeline summary
Show each file with ✓/✗ before the success/fail counts, so the user
can see at a glance which files succeeded and which failed.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
2026-05-14 20:49:30 +02:00
c6c804749e Skip SharedDEM computation when all visualizations already exist
Two optimizations to avoid ~2min wasted per file on re-runs:

1. pipeline.py: Check which visualizations need regeneration before
   computing SharedDEM. If all WebP outputs exist, skip SharedDEM
   entirely. If only IGN overlays need updating, also skip SharedDEM.

2. visualizations.py: Make SharedDEM attributes lazy (filled, gradient,
   lrm_15) so only the data actually needed is computed. For example,
   if only hillshade is regenerated, LRM at 15m is never calculated.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
2026-05-14 20:40:51 +02:00
cf3e680b02 Add download script for IGN LiDAR HD LAZ files
Supports tile download by coordinates, zone-based bulk download,
integrity checking, and auto-search across zones with resume support.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
2026-05-14 20:34:50 +02:00
02218b2cfc Upgrade PDAL to 2.10 via conda-forge, add COPC v1.1 support
- Dockerfile: install PDAL 2.10.1 from conda-forge (was 2.3 from apt)
  Ubuntu 22.04's PDAL 2.3 cannot read COPC v1.1 files from IGN LiDAR HD
- dtm.py: add _read_with_pdal() fallback for COPC files that laspy can't read
- dtm.py: validate_laz() now tries PDAL when laspy fails
- dtm.py: create_dtm_fast() and detect_ground_method() use PDAL fallback
- ign.py: auto-retry at lower zoom on 404 errors
- pipeline.py: check DTM resolution mismatch and regenerate if needed
- pipeline.py: propagate actual DTM resolution to visualizations
- pipeline.py: add --init to docker run for proper Ctrl+C signal handling
- Remove RRIM and Multi-Hillshade RGB visualizations

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
2026-05-14 19:01:05 +02:00
7ac08f75dc Add LAZ integrity check to skip corrupted files early
Validate file readability before PDAL classification. Corrupted/truncated
files are detected instantly via laspy header read and skipped with a clear
error message pointing to re-download, instead of wasting time on PDAL
and repair attempts that will fail anyway.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
2026-05-14 17:59:32 +02:00
f988ddb76d Auto-retry IGN tiles at lower zoom on 404
When zoom 20 tiles are unavailable (rural areas), fall back to zoom 19, 18,
etc. down to 15. Breaks out immediately on first-tile 404 to avoid wasting
requests at unsupported zoom levels.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
2026-05-14 02:21:47 +02:00
e2bd6b2536 Remove RRIM and Multi-Hillshade RGB, fix DTM resolution reuse bug, add --init to docker run
- Remove generate_rrim, generate_multi_hillshade, _compute_openness_both
- Remove corresponding VIZ_STEPS entries, COLORMAPS, RGB_LEGENDS, and tests
- Fix DTM resolution mismatch: existing DTM at different resolution is now
  regenerated instead of silently reused
- Propagate actual DTM resolution to visualizations and rendering
- Add --init to docker run commands for proper signal handling on Ctrl+C
- Add .playwright-mcp/ to .gitignore

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
2026-05-14 02:19:42 +02:00
bf17ca4662 Remove ETA display from visualization progress logs 2026-05-14 01:18:52 +02:00
08dbc73869 Fix lambda signatures for rrim and multi_hillshade to accept shared kwarg
Same issue as pos_open/neg_open — lambdas in VIZ_STEPS need shared=None
parameter since the pipeline passes shared=shared to all visualization functions.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
2026-05-14 01:07:48 +02:00
8c2065801b Add unit tests for RRIM, Multi-Hillshade RGB, and Local Dominance
TestRRIM: TIF generation, 3-band RGB output, uint8 range, no NaN
TestMultiHillshade: TIF generation, 3-band RGB output, uint8 range
TestLocalDominance: TIF generation, values in [0,1], NaN mask preservation

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
2026-05-14 01:05:55 +02:00
7f6b816ed6 Add RRIM, Multi-Hillshade RGB, and Local Dominance visualizations
Three new visualizations complementing existing SVF/openness/LRM/MSRM:

- RRIM (Red Relief Image Map): RGB composite combining positive openness
  (R), inverted slope (G), negative openness (B). Uses ray-tracing
  to compute both openness values in a single pass.

- Multi-Hillshade RGB: 3 azimuths (315°, 135°, 45°) mapped to R/G/B
  channels with slope blending. Color reveals structure orientation.

- Local Dominance: (dem - local_min) / (local_max - local_min) using
  min/max filters. Measures local height position — complements openness.

Also adds:
- _compute_openness_both() helper for shared ray-tracing (used by RRIM)
- xp_maximum_filter() in gpu.py (GPU/CPU abstraction)
- Entries in COLORMAPS, RGB_LEGENDS, VIZ_STEPS, and is_rgb detection
- All NaN handling follows existing patterns (nan_mask restoration)

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
2026-05-14 01:03:47 +02:00
1cf8e1752f Remove PMF, fix NaN in gradient visualizations, fix pos_open/neg_open shared param
- Remove PMF from ground classification options (PDAL recommends SMRF over PMF)
- Auto-detection now uses CSF for urban/complex terrain instead of PMF
- Add z_std > 30m heuristic to auto-select CSF for complex terrain
- Fix pos_open/neg_open lambda missing 'shared' parameter (NameError in workers)
- Fix NaN mask not restored in hillshade, slope, aspect, curvature
  (gradient-based products computed on filled DEM lost NaN transparency)
- Add nan_mask parameter to _save_tif for centralized NaN restoration
- DTM TIF kept by default (no longer deleted after WebP conversion)

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
2026-05-14 00:50:45 +02:00
eac482874d Fix bugs and improve pipeline flexibility
- Fix gpu_cleanup import missing in visualizations.py (NameError in workers)
- Fix t_pdf referenced before assignment when PDF is skipped
- Skip classification+DTM when DTM exists regardless of --force
- --force now only regenerates WebP/PDF, not classification/DTM
- --force-classification forces reclassification when needed
- Add laspy repair fallback for corrupt LAZ files (EVLR errors)
- Keep DTM TIF by default for reuse (--no-keep-tif to delete)
- Increase space between image and bottom cartouche (0.12→0.19)

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
2026-05-14 00:08:25 +02:00
5b74322077 Skip ground classification when DTM already exists
If the DTM .tif exists and --force is not set, skip both ground
classification and DTM generation entirely. Previously, the pipeline
would spend 3+ minutes reclassifying ground even when the DTM was
already present and would be reused anyway.

Also includes: SharedDEM cache, enhanced WebP cartouche (compass rose,
adaptive scale bar, enriched info bar), removed COG/viewer, UTF-8
fix for parallel workers, skip logic for DTM and PDF.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
2026-05-13 23:41:21 +02:00
15 changed files with 2142 additions and 1222 deletions

1
.gitignore vendored
View File

@ -45,3 +45,4 @@ htmlcov/
# Éventuels fichiers de cache matplotlib
matplotlibrc
.playwright-mcp/

View File

@ -4,7 +4,7 @@ This file provides guidance to Claude Code (claude.ai/code) when working with co
## Project Overview
LiDAR archaeological processing pipeline that generates 17 terrain visualizations from LAZ/LAS point clouds. Runs in Docker with optional NVIDIA GPU acceleration (CuPy). Designed for French LiDAR HD data in Lambert 93 (EPSG:2154).
LiDAR archaeological processing pipeline that generates 16 terrain visualizations from LAZ/LAS point clouds. Runs in Docker with optional NVIDIA GPU acceleration (CuPy). Designed for French LiDAR HD data in Lambert 93 (EPSG:2154).
## Commands
@ -14,12 +14,15 @@ All commands run inside Docker. Use `./run.sh` as the primary interface.
./run.sh -g # Standard run with GPU
./run.sh -g -w 4 # GPU + 4 parallel workers
./run.sh -g -r 0.2 # High resolution (0.2m/px)
./run.sh -g -r 0.5,0.2 # Multi-resolution (0.5m + 0.2m)
./run.sh --test # Run unit tests
./run.sh -g --file LHD_FXX_1000_6882_PTS_LAMB93_IGN69.copc # Single file
./run.sh --ground-classification pmf # Force PMF ground classification
./run.sh -g --keep-tif # Keep intermediate TIFF files
./run.sh -g --no-viewer # Skip web viewer generation
./run.sh serve # Start web map server
./run.sh --ground-classification csf # Force CSF ground classification (complex terrain)
./run.sh -g --keep-tif # Keep TIFF files (allows WebP regeneration without recalculating DTM)
./run.sh -g --only hillshade svf lrm # Only generate specific visualizations
./run.sh -g --skip ortho topo # Exclude specific visualizations
./run.sh -g --quality 90 # WebP quality 90 (default: 85)
./run.sh -g --lossless # Lossless WebP compression
./run.sh # Print help (no args)
```
@ -34,34 +37,50 @@ 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.
- **`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.
- **`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. Multi-resolution support: `self.resolutions` list, `_res_suffix()` for naming, `generate_all_visualizations()` accepts `vis_dir` override.
- **`dtm.py`** — PDAL ground classification (SMRF/CSF + auto-detection) and DTM generation via scipy `binned_statistic_2d`. `create_dtm_fast()` accepts `output_suffix` for multi-resolution DTM naming.
- **`visualizations.py`** — 13 `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. Lazy evaluation: properties computed on first access.
- **`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. Quality parameter controls WebP compression (default 85).
### 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
### Ground classification
Auto-detection in `dtm.py` `detect_ground_method()`:
- Single-return ratio > 0.6 → PMF (urban terrain)
- Single-return ratio > 0.6 → CSF (urban terrain, cloth simulation)
- Height std > 30m → CSF (complex/mountainous terrain)
- Default → SMRF (natural terrain)
Override with `--ground-classification {auto,smrf,pmf,csf}`.
Override with `--ground-classification {auto,smrf,csf}`.
### NaN handling
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.
### Multi-resolution
`-r 0.5,0.2` processes each tile at both 0.5m and 0.2m. Ground classification is shared (done once per tile). Each resolution gets its own DTM (`_dtm.tif` / `_dtm_r0p2.tif`) and visualization subdirectory (`basename/` / `basename_r0p2/`).
### Parallel processing
@ -72,7 +91,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 AVIF (quality 98 by default, best quality/size ratio). Use `--format webp` for WebP output. TIFF intermediates deleted by default. Use `--keep-tif` to keep DTM+TIF for regeneration with `--force`. No PDF reports, no COGs or viewer.
- **Compression**: TIF intermediates use `deflate` compression (faster than LZW for float32 data).
- **Tests**: Run only inside Docker via `./run.sh --test`. Synthetic DEM fixture in `tests/conftest.py`.

View File

@ -3,17 +3,21 @@ FROM nvidia/cuda:12.4.0-devel-ubuntu22.04
ENV DEBIAN_FRONTEND=noninteractive
ENV TZ=Europe/Paris
# Install PDAL and system packages
# Install system packages + Miniforge for PDAL >= 2.5 (Ubuntu 22.04 ships PDAL 2.3 which can't read COPC v1.1)
RUN apt-get update && apt-get install -y --no-install-recommends \
pdal \
liblaszip8 \
gdal-bin \
python3-gdal \
python3-pip \
python3-dev \
build-essential \
wget \
&& rm -rf /var/lib/apt/lists/*
&& rm -rf /var/lib/apt/lists/* \
&& wget -q https://github.com/conda-forge/miniforge/releases/latest/download/Miniforge3-Linux-x86_64.sh -O /tmp/miniforge.sh \
&& bash /tmp/miniforge.sh -b -p /opt/conda \
&& rm /tmp/miniforge.sh \
&& /opt/conda/bin/conda install -y -c conda-forge pdal \
&& ln -sf /opt/conda/bin/pdal /usr/local/bin/pdal \
&& /opt/conda/bin/conda clean -afy
WORKDIR /app
@ -37,7 +41,8 @@ RUN pip3 install --no-cache-dir \
titiler.core \
fastapi \
uvicorn \
piexif
piexif \
pillow-avif-plugin
# 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"

584
download_lidar.sh Executable file
View File

@ -0,0 +1,584 @@
#!/bin/bash
# download_lidar.sh — Téléchargement de fichiers LiDAR HD depuis l'IGN
#
# Usage:
# ./download_lidar.sh LHD_FXX_1049_6895_PTS_LAMB93_IGN69.copc.laz
# ./download_lidar.sh 1049_6895
# ./download_lidar.sh 1049_6895 1029_6884
# ./download_lidar.sh --list-zones
# ./download_lidar.sh --search 1049
# ./download_lidar.sh --zone RE --output input/
# ./download_lidar.sh --zone RE --dry-run
#
# L'API IGN expose les fichiers organisés par zones (RE, SE, AE, etc.)
# Le script cherche automatiquement dans quelle zone se trouve chaque tuile.
set -euo pipefail
BASE_URL="https://data.geopf.fr/telechargement/download/LiDARHD-NUALID"
API_URL="https://data.geopf.fr/telechargement/resource/LiDARHD-NUALID"
OUTPUT_DIR="${OUTPUT_DIR:-.}"
# Colors
RED='\033[0;31m'
GREEN='\033[0;32m'
YELLOW='\033[1;33m'
CYAN='\033[0;36m'
NC='\033[0m'
log_info() { echo -e "${CYAN}[INFO]${NC} $*"; }
log_ok() { echo -e "${GREEN}[OK]${NC} $*"; }
log_warn() { echo -e "${YELLOW}[WARN]${NC} $*"; }
log_err() { echo -e "${RED}[ERR]${NC} $*"; }
# Normalize a tile identifier to a full filename
# Supports: "1049_6895", "LHD_FXX_1049_6895_PTS_LAMB93_IGN69", "LHD_FXX_1049_6895_PTS_LAMB93_IGN69.copc.laz"
normalize_name() {
local input="$1"
# If it's already a full filename, return as-is
if [[ "$input" == LHD_FXX_*_PTS_LAMB93_IGN69.copc.laz ]]; then
echo "$input"
return
fi
if [[ "$input" == LHD_FXX_*_PTS_LAMB93_IGN69 ]]; then
echo "${input}.copc.laz"
return
fi
# Bare coordinates like "1049_6895"
echo "LHD_FXX_${input}_PTS_LAMB93_IGN69.copc.laz"
}
# List all available zone releases
list_zones() {
log_info "Récupération de la liste des zones..."
local page=1
local total_pages=""
while true; do
local content
content=$(curl -s "${API_URL}?page=${page}&limit=50") || true
if [ -z "$content" ]; then
break
fi
# Parse zone codes and dates
echo "$content" | python3 -c "
import sys, xml.etree.ElementTree as ET
ns = {'atom': 'http://www.w3.org/2005/Atom'}
try:
tree = ET.parse(sys.stdin)
root = tree.getroot()
for entry in root.findall('atom:entry', ns):
title = entry.find('atom:title', ns)
if title is not None and title.text:
# Extract zone code and date
parts = title.text.split('_')
# NUALHD_1-0__LAZ_LAMB93_{ZONE}_{DATE}
zone = parts[-2] if len(parts) >= 2 else '?'
date = parts[-1] if len(parts) >= 1 else '?'
print(f'{zone}\t{date}\t{title.text}')
except Exception:
pass
" 2>/dev/null || true
# Check if more pages
if [ -z "$total_pages" ]; then
total_pages=$(echo "$content" | python3 -c "
import sys, xml.etree.ElementTree as ET
ns = {'atom': 'http://www.w3.org/2005/Atom', 'gpf_dl': 'https://data.geopf.fr/annexes/ressources/xsd/gpf_dl.xsd'}
try:
tree = ET.parse(sys.stdin)
root = tree.getroot()
pc = root.get('{https://data.geopf.fr/annexes/ressources/xsd/gpf_dl.xsd}pagecount', '1')
print(pc)
except:
print('1')
" 2>/dev/null)
fi
if [ "$page" -ge "${total_pages:-1}" ]; then
break
fi
page=$((page + 1))
done | sort
}
# Search for a tile across all zones
search_tile() {
local tile_name="$1"
local filename
filename=$(normalize_name "$tile_name")
log_info "Recherche de ${filename} dans toutes les zones..."
local page=1
local found=0
local total_pages=""
while true; do
local content
content=$(curl -s "${API_URL}?page=${page}&limit=50") || true
[ -z "$content" ] && break
# Get sub-resources and search for the file
echo "$content" | python3 -c "
import sys, xml.etree.ElementTree as ET
ns = {'atom': 'http://www.w3.org/2005/Atom'}
try:
tree = ET.parse(sys.stdin)
root = tree.getroot()
for entry in root.findall('atom:entry', ns):
title = entry.find('atom:title', ns)
if title is not None and title.text:
print(title.text)
except:
pass
" 2>/dev/null | while read -r release; do
# Check if file exists in this release
local url="${BASE_URL}/${release}/${filename}"
local code
code=$(curl -sI -L -o /dev/null -w "%{http_code}" "$url" 2>/dev/null)
if [ "$code" = "200" ]; then
log_ok "Trouvé: ${filename} dans la zone ${release}"
echo " URL: ${url}"
found=1
fi
done
if [ -z "$total_pages" ]; then
total_pages=$(echo "$content" | python3 -c "
import sys, xml.etree.ElementTree as ET
ns = {'atom': 'http://www.w3.org/2005/Atom', 'gpf_dl': 'https://data.geopf.fr/annexes/ressources/xsd/gpf_dl.xsd'}
try:
tree = ET.parse(sys.stdin)
root = tree.getroot()
pc = root.get('{https://data.geopf.fr/annexes/ressources/xsd/gpf_dl.xsd}pagecount', '1')
print(pc)
except:
print('1')
" 2>/dev/null)
fi
if [ "$page" -ge "${total_pages:-1}" ]; then
break
fi
page=$((page + 1))
done
if [ "$found" -eq 0 ]; then
log_err "${filename} non trouvé dans aucune zone"
fi
}
# Download a single file with resume support
download_file() {
local url="$1"
local output_path="$2"
local filename
filename=$(basename "$output_path")
# Resume download with retries
local max_retries=15
local retry=0
while [ $retry -lt $max_retries ]; do
retry=$((retry + 1))
local result
result=$(curl -L --http1.1 -C - -o "$output_path" "$url" 2>&1)
local exit_code=$?
if [ $exit_code -eq 0 ]; then
# Verify it's a valid LAS/LAZ file
local magic
magic=$(head -c 4 "$output_path" 2>/dev/null | xxd -p 2>/dev/null || echo "")
if [ "$magic" = "4c415346" ]; then
local size
size=$(stat -c%s "$output_path" 2>/dev/null || echo "0")
log_ok "${filename} téléchargé (${size} octets)"
return 0
else
log_err "${filename} n'est pas un fichier LAZ valide (magic: ${magic})"
rm -f "$output_path"
return 1
fi
fi
# curl exit code 18 = partial download, resume
if [ $exit_code -eq 18 ]; then
local size
size=$(stat -c%s "$output_path" 2>/dev/null || echo "0")
log_warn "Téléchargement partiel (${size} octets), tentative ${retry}/${max_retries}..."
sleep 2
continue
fi
# HTTP 404
if echo "$result" | grep -q "404"; then
log_err "Fichier non trouvé (404): ${url}"
rm -f "$output_path"
return 1
fi
log_err "Erreur curl (${exit_code}): ${result}"
return 1
done
# If we exhausted retries, check if partial file is usable
local magic
magic=$(head -c 4 "$output_path" 2>/dev/null | xxd -p 2>/dev/null || echo "")
if [ "$magic" = "4c415346" ]; then
local size
size=$(stat -c%s "$output_path" 2>/dev/null || echo "0")
log_warn "Fichier partiel mais valide: ${filename} (${size} octets)"
return 0
fi
log_err "Échec après ${max_retries} tentatives"
rm -f "$output_path"
return 1
}
# Find which zone contains a given tile and download it
find_and_download() {
local tile_name="$1"
local output_dir="$2"
local filename
filename=$(normalize_name "$tile_name")
local output_path="${output_dir}/${filename}"
# If file already exists and is valid, skip
if [ -f "$output_path" ]; then
local magic
magic=$(head -c 4 "$output_path" 2>/dev/null | xxd -p 2>/dev/null || echo "")
if [ "$magic" = "4c415346" ]; then
local size
size=$(stat -c%s "$output_path" 2>/dev/null || echo "0")
if [ "$size" -gt 1000000 ]; then
log_info "${filename} déjà présent (${size} octets) — ignoré"
return 0
fi
fi
# File exists but seems invalid/truncated, re-download
log_warn "${filename} existe mais semble incomplet, re-téléchargement..."
fi
log_info "Recherche de ${filename}..."
# Try all zone releases
local page=1
local total_pages=""
while true; do
local content
content=$(curl -s "${API_URL}?page=${page}&limit=50") || true
[ -z "$content" ] && break
echo "$content" | python3 -c "
import sys, xml.etree.ElementTree as ET
ns = {'atom': 'http://www.w3.org/2005/Atom'}
try:
tree = ET.parse(sys.stdin)
root = tree.getroot()
for entry in root.findall('atom:entry', ns):
title = entry.find('atom:title', ns)
if title is not None and title.text:
print(title.text)
except:
pass
" 2>/dev/null | while read -r release; do
local url="${BASE_URL}/${release}/${filename}"
local code
code=$(curl -sI -L -o /dev/null -w "%{http_code}" "$url" 2>/dev/null || echo "000")
if [ "$code" = "200" ]; then
log_ok "Trouvé dans la zone ${release}"
download_file "$url" "$output_path"
return 0
fi
done
if [ -z "$total_pages" ]; then
total_pages=$(echo "$content" | python3 -c "
import sys, xml.etree.ElementTree as ET
ns = {'atom': 'http://www.w3.org/2005/Atom', 'gpf_dl': 'https://data.geopf.fr/annexes/ressources/xsd/gpf_dl.xsd'}
try:
tree = ET.parse(sys.stdin)
root = tree.getroot()
pc = root.get('{https://data.geopf.fr/annexes/ressources/xsd/gpf_dl.xsd}pagecount', '1')
print(pc)
except:
print('1')
" 2>/dev/null)
fi
if [ "$page" -ge "${total_pages:-1}" ]; then
break
fi
page=$((page + 1))
done
log_err "${filename} non trouvé dans aucune zone disponible"
return 1
}
# Download all tiles from a specific zone
download_zone() {
local zone="$1"
local output_dir="$2"
local dry_run="${3:-false}"
# Find the release name for this zone
local release=""
local page=1
while true; do
local content
content=$(curl -s "${API_URL}?page=${page}&limit=50") || true
[ -z "$content" ] && break
release=$(echo "$content" | python3 -c "
import sys, xml.etree.ElementTree as ET
ns = {'atom': 'http://www.w3.org/2005/Atom'}
try:
tree = ET.parse(sys.stdin)
root = tree.getroot()
for entry in root.findall('atom:entry', ns):
title = entry.find('atom:title', ns)
if title is not None and title.text:
# Extract zone code from title like NUALHD_1-0__LAZ_LAMB93_RE_2025-02-17
parts = title.text.split('_')
z = parts[-2] if len(parts) >= 2 else ''
if z == '${zone}':
print(title.text)
break
except:
pass
" 2>/dev/null)
if [ -n "$release" ]; then
break
fi
page=$((page + 1))
# Safety limit
if [ $page -gt 25 ]; then
break
fi
done
if [ -z "$release" ]; then
log_err "Zone '${zone}' non trouvée"
return 1
fi
log_info "Zone ${zone}: ${release}"
# Get list of files in this release
local sub_url="${API_URL}/${release}?limit=50"
local total_files=0
local page=1
while true; do
local content
content=$(curl -s "${API_URL}/${release}?page=${page}&limit=50") || true
[ -z "$content" ] && break
local files
files=$(echo "$content" | python3 -c "
import sys, xml.etree.ElementTree as ET
ns = {'atom': 'http://www.w3.org/2005/Atom'}
try:
tree = ET.parse(sys.stdin)
root = tree.getroot()
for entry in root.findall('atom:entry', ns):
title = entry.find('atom:title', ns)
if title is not None and title.text:
print(title.text)
except:
pass
" 2>/dev/null)
if [ -z "$files" ]; then
break
fi
local count
count=$(echo "$files" | wc -l)
total_files=$((total_files + count))
if [ "$dry_run" = "true" ]; then
echo "$files" | head -10
echo "... (${count} fichiers sur cette page)"
else
echo "$files" | while read -r fname; do
local url="${BASE_URL}/${release}/${fname}"
download_file "$url" "${output_dir}/${fname}"
done
fi
page=$((page + 1))
# Safety limit
if [ $page -gt 500 ]; then
break
fi
done
log_info "Total: ${total_files} fichiers dans la zone ${zone}"
}
# Show usage
usage() {
cat <<'EOF'
download_lidar.sh — Téléchargement de fichiers LiDAR HD depuis l'IGN
Usage:
./download_lidar.sh <tuile> [tuile2...] Télécharger une ou plusieurs tuiles
./download_lidar.sh --search <motif> Chercher une tuile dans toutes les zones
./download_lidar.sh --list-zones Lister les zones disponibles
./download_lidar.sh --zone <ZONE> [--dry-run] Télécharger une zone entière
./download_lidar.sh --check [dossier/] Vérifier l'intégrité des fichiers LAZ
Formats de tuile acceptés:
1049_6895 Coordonnées seules
LHD_FXX_1049_6895_PTS_LAMB93_IGN69 Nom complet sans extension
LHD_FXX_1049_6895_PTS_LAMB93_IGN69.copc.laz Nom complet avec extension
Options:
-o, --output DIR Répertoire de destination (défaut: .)
-n, --dry-run Afficher sans télécharger (avec --zone)
--check [DIR] Vérifier l'intégrité des fichiers LAZ (défaut: ./)
--list-zones Lister toutes les zones disponibles
--search MOTIF Chercher un motif dans les noms de tuile
--zone ZONE Télécharger tous les fichiers d'une zone
-h, --help Afficher cette aide
Exemples:
./download_lidar.sh 1049_6895
./download_lidar.sh LHD_FXX_0713_6347_PTS_LAMB93_IGN69.copc.laz -o input/
./download_lidar.sh --search 1049
./download_lidar.sh --list-zones
./download_lidar.sh --zone RE -o input/
./download_lidar.sh --check input/
EOF
}
# Check integrity of LAZ files
check_files() {
local dir="${1:-.}"
local ok=0
local fail=0
log_info "Vérification de l'intégrité des fichiers LAZ dans ${dir}/"
for f in "${dir}"/*.copc.laz "${dir}"/*.laz; do
[ -f "$f" ] || continue
local magic
magic=$(head -c 4 "$f" 2>/dev/null | xxd -p 2>/dev/null || echo "")
local size
size=$(stat -c%s "$f" 2>/dev/null || echo "0")
if [ "$magic" = "4c415346" ] && [ "$size" -gt 1000000 ]; then
log_ok "$(basename $f) (${size} octets)"
ok=$((ok + 1))
elif [ "$magic" = "7b227469" ]; then
log_err "$(basename $f) est une page HTML (404), pas un fichier LAZ"
fail=$((fail + 1))
elif [ "$size" -lt 1000 ]; then
log_err "$(basename $f) trop petit (${size} octets) — probablement corrompu"
fail=$((fail + 1))
else
log_warn "$(basename $f) magic=${magic}, taille=${size} — vérification nécessaire"
fail=$((fail + 1))
fi
done
log_info "Résultat: ${ok} OK, ${fail} problème(s)"
return 0
}
# Main
CMD=""
TILES=()
ZONE=""
DRY_RUN=false
while [ $# -gt 0 ]; do
case "$1" in
-h|--help)
usage
exit 0
;;
--list-zones)
CMD="list-zones"
shift
;;
--search)
CMD="search"
shift
TILES+=("$1")
shift
;;
--zone)
CMD="zone"
shift
ZONE="$1"
shift
;;
--check)
CMD="check"
shift
OUTPUT_DIR="${1:-.}"
[ -d "$OUTPUT_DIR" ] && shift || true
;;
-o|--output)
shift
OUTPUT_DIR="$1"
shift
;;
-n|--dry-run)
DRY_RUN=true
shift
;;
-*)
log_err "Option inconnue: $1"
usage
exit 1
;;
*)
TILES+=("$1")
shift
;;
esac
done
case "${CMD}" in
list-zones)
list_zones
;;
search)
for tile in "${TILES[@]}"; do
search_tile "$tile"
done
;;
zone)
if [ -z "$ZONE" ]; then
log_err "Zone non spécifiée. Utilisez --zone <CODE>"
exit 1
fi
mkdir -p "$OUTPUT_DIR"
download_zone "$ZONE" "$OUTPUT_DIR" "$DRY_RUN"
;;
check)
check_files "$OUTPUT_DIR"
;;
*)
# Default: download tiles
if [ ${#TILES[@]} -eq 0 ]; then
usage
exit 1
fi
mkdir -p "$OUTPUT_DIR"
for tile in "${TILES[@]}"; do
find_and_download "$tile" "$OUTPUT_DIR"
done
;;
esac

View File

@ -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"
@ -91,9 +97,9 @@ Exemples:
)
parser.add_argument(
"-r", "--resolution",
type=float,
default=0.5,
help="Résolution en mètres par pixel (défaut: 0.5)"
type=str,
default="0.5",
help="Résolution en m/px, ou multiples séparées par virgules (défaut: 0.5, ex: 0.5,0.2)"
)
parser.add_argument(
"-w", "--workers",
@ -114,18 +120,44 @@ Exemples:
parser.add_argument(
"--keep-tif",
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)"
help="Conserver les fichiers TIFF (DTM + visualisations) pour pouvoir régénérer les WebP sans recalculer"
)
parser.add_argument(
"--ground-classification",
choices=["auto", "smrf", "pmf", "csf"],
choices=["auto", "smrf", "csf"],
default="auto",
help="Méthode de classification du sol : auto (détection), smrf, pmf, csf (défaut: auto)"
help="Méthode de classification du sol : auto (détection), smrf, csf (défaut: auto)"
)
parser.add_argument(
"--quality",
type=int,
default=98,
help="Qualité image (1-100, défaut: 98). Utilisez 100 pour lossless."
)
parser.add_argument(
"--lossless",
action="store_true",
help="Forcer la compression lossless (équivalent à --quality 100)"
)
parser.add_argument(
"--format",
choices=["webp", "avif"],
default="avif",
help="Format de sortie : avif (défaut, meilleure qualité) ou webp"
)
parser.add_argument(
"--only",
nargs="+",
type=str,
default=None,
help="Générer uniquement ces visualisations (ex: --only hillshade svf lrm)"
)
parser.add_argument(
"--skip",
nargs="+",
type=str,
default=None,
help="Exclure ces visualisations (ex: --skip ortho topo)"
)
parser.add_argument(
"--file",
@ -167,6 +199,14 @@ Exemples:
log_gpu_status()
try:
quality = 100 if args.lossless else args.quality
# Parse --only and --skip: accept comma-separated values
only_viz = None
if args.only:
only_viz = [v.strip() for item in args.only for v in item.split(',')]
skip_viz = None
if args.skip:
skip_viz = [v.strip() for item in args.skip for v in item.split(',')]
pipeline = LidarArchaeoPipeline(
input_dir=args.input,
output_dir=args.output,
@ -176,7 +216,10 @@ Exemples:
ground_method=args.ground_classification,
force_classify=args.force_classification,
keep_tif=args.keep_tif,
no_viewer=args.no_viewer
quality=quality,
only_viz=only_viz,
skip_viz=skip_viz,
output_format=args.format,
)
# If --file is specified, process only matching files
@ -187,6 +230,12 @@ Exemples:
# Also supports bare name without .copc (e.g. LHD_FXX_1000_6882_PTS_LAMB93_IGN69)
selected_files = []
for pattern in args.file:
# Try exact filename first (e.g. LHD_FXX_...copc.laz)
exact_match = input_dir / pattern
if exact_match.exists() and exact_match.is_file():
matches = [exact_match]
else:
# Try with added extensions (e.g. pattern=LHD_FXX_...IGN69)
matches = (list(input_dir.glob(f"{pattern}.laz"))
+ list(input_dir.glob(f"{pattern}.las"))
+ list(input_dir.glob(f"{pattern}.copc.laz"))
@ -235,9 +284,15 @@ def _kill_orphan_pdal(signum=None, frame=None):
"""Kill orphan PDAL processes on interrupt or exit."""
import subprocess
try:
subprocess.run(["pkill", "-f", "pdal"], capture_output=True, timeout=5)
subprocess.run(["pkill", "-9", "-f", "pdal"], capture_output=True, timeout=3)
except Exception:
pass
if signum is not None:
logger.info("Interruption — nettoyage des processus PDAL")
logger.info("Interruption — nettoyage des processus")
# Force-kill all child processes immediately
try:
import os
os.killpg(os.getpgrp(), signal.SIGKILL)
except Exception:
pass
sys.exit(130)

View File

@ -1,6 +1,6 @@
"""DTM generation from classified LiDAR point clouds.
Handles ground classification via PDAL/SMRF and DTM rasterisation
Handles ground classification via PDAL (SMRF or CSF) and DTM rasterisation
using scipy binned_statistic_2d. Zones without LiDAR data remain as NaN.
"""
@ -27,13 +27,13 @@ def _create_ground_pipeline(input_laz, output_las, method):
1. Reset Classification to 0
2. ELM (Extended Local Minimum) — mark low outliers as noise (Classification=7)
3. Statistical outlier removal
4. Ground classification (SMRF/PMF/CSF)
4. Ground classification (SMRF or CSF)
5. Extract ground points (Classification=2)
Args:
input_laz: Path to input LAZ/LAS file.
output_las: Path to output classified LAS file.
method: Ground classification method ('smrf', 'pmf', or 'csf').
method: Ground classification method ('smrf' or 'csf').
Returns:
JSON string of the PDAL pipeline.
@ -84,15 +84,6 @@ def _create_ground_pipeline(input_laz, output_las, method):
"threshold": 0.5,
"scalar": 1.25
}
elif method == 'pmf':
ground_step = {
"type": "filters.pmf",
"max_window": 33,
"slope": 0.15,
"max_distance": 2.5,
"initial_distance": 0.15,
"cell_size": 1.0
}
elif method == 'csf':
ground_step = {
"type": "filters.csf",
@ -128,22 +119,118 @@ def create_smrf_pipeline(input_laz, output_las):
return _create_ground_pipeline(input_laz, output_las, 'smrf')
def create_pmf_pipeline(input_laz, output_las):
"""Create a PDAL pipeline JSON for PMF ground classification."""
return _create_ground_pipeline(input_laz, output_las, 'pmf')
def create_csf_pipeline(input_laz, output_las):
"""Create a PDAL pipeline JSON for CSF ground classification."""
return _create_ground_pipeline(input_laz, output_las, 'csf')
def validate_laz(laz_file):
"""Quick integrity check for a LAZ/LAS file.
Tries laspy first (fast header read), then PDAL as fallback for COPC files
that laspy cannot read. Also checks that the file contains points.
Returns:
True if file is readable and contains points, False otherwise.
"""
# Try laspy first (fast)
import laspy
try:
with laspy.open(str(laz_file)) as f:
header = f.header
point_count = header.point_count
if point_count == 0:
logger.error(f" ✗ Fichier vide (0 points): {laz_file.name}")
logger.error(f" → Re-télécharger depuis https://ign.fr/lidar-hd")
return False
return True
except Exception:
pass
# Fallback: try PDAL (handles COPC v1.1 that laspy can't read)
import subprocess
try:
result = subprocess.run(
["pdal", "info", str(laz_file), "--summary"],
capture_output=True, text=True, timeout=30
)
if result.returncode == 0:
# Check point count from PDAL info output
import json as _json
try:
info = _json.loads(result.stdout)
count = info.get('summary', {}).get('num_points', 0)
if count == 0:
logger.error(f" ✗ Fichier vide (0 points PDAL): {laz_file.name}")
logger.error(f" → Re-télécharger depuis https://ign.fr/lidar-hd")
return False
except Exception:
pass # Can't parse — assume valid
return True
logger.error(f" ✗ Fichier illisible: {laz_file.name}")
logger.error(f" PDAL: {result.stderr.strip()[:200]}")
except (subprocess.TimeoutExpired, FileNotFoundError):
logger.error(f" ✗ Impossible de vérifier le fichier: {laz_file.name}")
logger.error(f" → Re-télécharger depuis https://ign.fr/lidar-hd")
return False
def _read_with_pdal(laz_file):
"""Read a LAZ/LAS file via PDAL when laspy fails (e.g. COPC v1.1).
Returns a laspy.LasData object, or None on failure.
"""
import subprocess
import tempfile
import os
try:
# Convert COPC to LAS via PDAL, then read with laspy
with tempfile.NamedTemporaryFile(suffix='.las', delete=False) as tmp:
tmp_path = tmp.name
pipeline = json.dumps({
"pipeline": [
str(laz_file),
{"type": "writers.las", "filename": tmp_path}
]
})
result = subprocess.run(
["pdal", "pipeline", "--stdin"],
input=pipeline, capture_output=True, text=True, timeout=300
)
if result.returncode != 0:
logger.warning(f" PDAL conversion échouée: {result.stderr[:200]}")
try:
os.unlink(tmp_path)
except Exception:
pass
return None
import laspy
las = laspy.read(tmp_path)
try:
os.unlink(tmp_path)
except Exception:
pass
if len(las.points) == 0:
logger.warning(f" PDAL: conversion réussie mais 0 points")
return None
return las
except Exception as e:
logger.warning(f" PDAL fallback échoué: {e}")
return None
def detect_ground_method(laz_file):
"""Detect the best ground classification method based on point cloud statistics.
Auto-selects between SMRF (natural terrain) and PMF (urban) only.
CSF is available only via --ground-classification csf (slow but robust
on complex terrain).
Auto-selects between SMRF and CSF:
- SMRF: fast, robust for most natural terrain (PDAL recommended default)
- CSF: cloth simulation, better for complex/urban terrain
Falls back to SMRF if the file cannot be read or attributes are missing.
@ -151,18 +238,28 @@ def detect_ground_method(laz_file):
laz_file: Path to input LAZ/LAS file.
Returns:
String: 'smrf', 'pmf', or 'csf'
String: 'smrf' or 'csf'
"""
import laspy
# Try laspy first, then PDAL for COPC files
las = None
try:
las = laspy.read(str(laz_file))
except Exception as e:
logger.warning(f" Impossible de lire le nuage pour auto-détection: {e}")
logger.warning(f" laspy: {e}")
logger.info(f" → Lecture via PDAL pour auto-détection...")
las = _read_with_pdal(laz_file)
if las is None:
logger.info(f" → Méthode: SMRF (défaut — lecture impossible)")
return 'smrf'
total_points = len(las.points)
if total_points == 0:
logger.warning(f" Nuage vide (0 points) — méthode par défaut: SMRF")
return 'smrf'
z = np.array(las.z)
# Height variance (always available)
@ -182,13 +279,16 @@ def detect_ground_method(laz_file):
f"ratio_retours_uniques={single_return_ratio:.2f}, "
f"écart_Z={z_std:.1f}m, amplitude_Z={z_range:.1f}m")
# Decision logic (auto selects between SMRF and PMF only):
# - High single-return ratio (>0.6) → urban (buildings, roads) → PMF
# Decision logic:
# - High single-return ratio (>0.6) → urban (buildings, roads) → CSF (cloth simulation)
# - High elevation variance (>30m) → complex/mountainous terrain → CSF
# - Default → SMRF (fast, robust for most natural terrain)
# Note: CSF is available only via --ground-classification csf (slow but robust on complex terrain)
if single_return_ratio > 0.6:
method = 'pmf'
method = 'csf'
reason = f"ratio retours uniques={single_return_ratio:.2f} > 0.6 → milieu urbain"
elif z_std > 30:
method = 'csf'
reason = f"écart_Z={z_std:.1f}m > 30m → terrain complexe"
else:
method = 'smrf'
reason = f"terrain naturel standard"
@ -203,7 +303,7 @@ def classify_ground(laz_file, temp_dir, method='auto', force=False):
Args:
laz_file: Path to input LAZ/LAS file.
temp_dir: Directory for temporary files (pipeline.json, ground.las).
method: Ground classification method ('auto', 'smrf', 'pmf', or 'csf').
method: Ground classification method ('auto', 'smrf', or 'csf').
force: If True, reclassify even if output file already exists.
Returns:
@ -243,14 +343,66 @@ def classify_ground(laz_file, temp_dir, method='auto', force=False):
["pdal", "pipeline", str(pipeline_file)],
capture_output=True, check=True
)
# Verify that ground file has points
if output_las.exists() and output_las.stat().st_size < 100:
logger.error(f" ✗ Fichier ground vide (taille < 100 octets)")
output_las.unlink(missing_ok=True)
return None
logger.info(f" ✓ Classification sol {method.upper()} terminée")
return output_las
except subprocess.CalledProcessError as e:
logger.error(f" ✗ Erreur classification PDAL ({method.upper()}): {e.stderr.decode()}")
error_msg = e.stderr.decode() if e.stderr else str(e)
logger.warning(f" ✗ Erreur classification PDAL ({method.upper()}): {error_msg}")
# Try repairing file with laspy if PDAL fails on EVLR/VLR
if 'VLR' in error_msg or 'Invalid' in error_msg:
logger.info(f" → Tentative de réparation du fichier avec laspy...")
repaired_las = temp_dir / f"{laz_base}_repaired.las"
if _repair_laz_with_laspy(laz_file, repaired_las):
# Retry PDAL pipeline with repaired file
pipeline_json = _create_ground_pipeline(repaired_las, output_las, method)
with open(pipeline_file, 'w') as f:
f.write(pipeline_json)
try:
subprocess.run(
["pdal", "pipeline", str(pipeline_file)],
capture_output=True, check=True
)
logger.info(f" ✓ Classification sol {method.upper()} terminée (fichier réparé)")
return output_las
except subprocess.CalledProcessError as e2:
error_msg2 = e2.stderr.decode() if e2.stderr else str(e2)
logger.error(f" ✗ Échec classification même après réparation: {error_msg2}")
else:
logger.error(f" ✗ Impossible de réparer le fichier")
return None
def create_dtm_fast(las_file, basename, dtm_dir, resolution):
def _repair_laz_with_laspy(input_laz, output_las):
"""Try to repair a corrupt LAZ file by re-reading with laspy and saving as LAS.
Works around PDAL errors like 'Invalid Extended VLR size' by stripping
problematic VLR/EVLR metadata during re-save.
Args:
input_laz: Path to corrupt LAZ/LAS file.
output_las: Path for repaired LAS output.
Returns:
True if repair succeeded, False otherwise.
"""
import laspy
try:
las = laspy.read(str(input_laz))
las.write(str(output_las))
logger.info(f" ✓ Fichier réparé via laspy ({len(las.points):,} points)")
return True
except Exception as e:
logger.warning(f" ✗ Réparation laspy échouée: {e}")
return False
def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, output_suffix=""):
"""Create DTM using fast binning method with gap filling.
Args:
@ -258,16 +410,38 @@ 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.
output_suffix: Suffix for output filename (e.g. '_r0p2' for additional resolutions).
Returns:
Path to output DTM GeoTIFF, or None on failure.
"""
output_tif = dtm_dir / f"{basename}_dtm{output_suffix}.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...")
try:
las = laspy.read(str(las_file))
except Exception as e:
# laspy can't read COPC v1.1 — try PDAL conversion
logger.warning(f" laspy: {e}")
logger.info(f" → Conversion via PDAL pour lecture COPC...")
las = _read_with_pdal(las_file)
if las is None:
logger.error(f" ✗ Impossible de lire {las_file.name}")
return None
if len(las.points) == 0:
logger.error(f" ✗ Fichier vide (0 points): {las_file.name}")
return None
try:
min_x, max_x = float(las.header.min[0]), float(las.header.max[0])
min_y, max_y = float(las.header.min[1]), float(las.header.max[1])
@ -309,7 +483,7 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution):
logger.info(f" {remaining:,} pixels restent sans données (grands écarts)")
# Save as GeoTIFF
output_tif = dtm_dir / f"{basename}_dtm.tif"
output_tif = dtm_dir / f"{basename}_dtm{output_suffix}.tif"
transform = from_bounds(min_x, min_y, max_x, max_y, width, height)
with rasterio.open(

View File

@ -110,6 +110,16 @@ def xp_minimum_filter(arr, footprint=None, size=None):
return ndimage.minimum_filter(arr, footprint=footprint, size=size)
def xp_maximum_filter(arr, footprint=None, size=None):
"""Maximum filter — uses GPU if array is on GPU, CPU otherwise."""
if _cp is not None and isinstance(arr, _cp.ndarray):
try:
return _cp_ndimage.maximum_filter(arr, footprint=footprint, size=size)
except Exception:
arr = to_cpu(arr)
return ndimage.maximum_filter(arr, footprint=footprint, size=size)
def gpu_cleanup():
"""Free GPU memory. Call between visualizations to prevent OOM."""
if _cp is not None:

View File

@ -76,6 +76,8 @@ def _lat_lon_to_px(lat, lon, zoom, tile_size=256):
def download_ign_tiles(min_x, max_x, min_y, max_y, layer, zoom_level=15):
"""Download IGN WMTS tiles for the given bounds using Web Mercator (PM).
If the first tile returns 404, automatically retries at lower zoom levels.
Args:
min_x, max_x, min_y, max_y: Bounds in Lambert 93.
layer: IGN WMTS layer name.
@ -106,25 +108,32 @@ def download_ign_tiles(min_x, max_x, min_y, max_y, layer, zoom_level=15):
tile_matrix_set = "PM"
tile_size = 256
col_min, row_min = _lat_lon_to_tile(nw_lat, nw_lon, zoom_level)
col_max, row_max = _lat_lon_to_tile(se_lat, se_lon, zoom_level)
# Try downloading at the requested zoom level; fall back to lower zooms on 404
min_zoom = 15
for zoom in range(zoom_level, min_zoom - 1, -1):
col_min, row_min = _lat_lon_to_tile(nw_lat, nw_lon, zoom)
col_max, row_max = _lat_lon_to_tile(se_lat, se_lon, zoom)
nw_px_x, nw_px_y = _lat_lon_to_px(nw_lat, nw_lon, zoom_level)
se_px_x, se_px_y = _lat_lon_to_px(se_lat, se_lon, zoom_level)
nw_px_x, nw_px_y = _lat_lon_to_px(nw_lat, nw_lon, zoom)
se_px_x, se_px_y = _lat_lon_to_px(se_lat, se_lon, zoom)
out_width = int(se_px_x - nw_px_x)
out_height = int(se_px_y - nw_px_y)
if out_width <= 0 or out_height <= 0 or out_width > 10000 or out_height > 10000:
logger.warning(f" Image IGN trop grande: {out_width}x{out_height}px — abandon")
return None
logger.warning(f" Image IGN trop grande: {out_width}x{out_height}px — zoom {zoom} abandon")
continue
total_tiles = (col_max - col_min + 1) * (row_max - row_min + 1)
logger.info(f" Zoom {zoom_level}: {total_tiles} tuiles à télécharger ({out_width}x{out_height}px)")
if zoom < zoom_level:
logger.info(f" ↓ Retry zoom {zoom_level}{zoom}: {total_tiles} tuiles ({out_width}x{out_height}px)")
else:
logger.info(f" Zoom {zoom}: {total_tiles} tuiles à télécharger ({out_width}x{out_height}px)")
composite = np.full((out_height, out_width, 3), 255, dtype=np.uint8)
tiles_downloaded = 0
tiles_404 = 0
fmt = "image/png" if 'PLAN' in layer else "image/jpeg"
for col in range(col_min, col_max + 1):
@ -133,7 +142,7 @@ def download_ign_tiles(min_x, max_x, min_y, max_y, layer, zoom_level=15):
f"{wmts_url}?SERVICE=WMTS&VERSION=1.0.0&REQUEST=GetTile"
f"&LAYER={layer}&STYLE=normal"
f"&TILEMATRIXSET={tile_matrix_set}"
f"&TILEMATRIX={zoom_level}&TILECOL={col}&TILEROW={row}"
f"&TILEMATRIX={zoom}&TILECOL={col}&TILEROW={row}"
f"&FORMAT={fmt}"
)
@ -167,16 +176,33 @@ def download_ign_tiles(min_x, max_x, min_y, max_y, layer, zoom_level=15):
tile_arr[src_y:src_y+src_h, src_x:src_x+src_w]
tiles_downloaded += 1
except Exception as e:
if tiles_downloaded == 0 and col == col_min and row == row_min:
logger.error(f" ✗ Erreur tuile IGN: {e}")
except urllib.error.HTTPError as e:
if e.code == 404:
tiles_404 += 1
# If the very first tile is 404, this zoom level is unavailable
if col == col_min and row == row_min:
logger.info(f" Zoom {zoom} non disponible (404) — essai zoom inférieur")
break
continue
except Exception:
continue
else:
continue
# Only reach here if inner loop broke (first tile 404)
break
if tiles_404 > 0 and tiles_downloaded == 0:
# No tiles at this zoom, try lower
continue
logger.info(f"{tiles_downloaded} tuiles IGN téléchargées ({layer})")
if tiles_downloaded == 0:
return None
continue
return composite
logger.error(" ✗ Aucun zoom disponible pour cette zone")
return None
def generate_ign_overlay(dem_file, basename, vis_dir, resolution, layer, title, legend_label, description, out_suffix):
"""Generate an IGN basemap overlay visualization.

View File

@ -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,16 +56,16 @@ _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_lrm, generate_openness,
generate_mslrm, generate_tpi, generate_sailore,
generate_roughness, generate_anomalies, generate_wavelet,
generate_flow,
)
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
# Ordered list of visualization steps.
@ -75,10 +76,9 @@ VIZ_STEPS = [
('slope', generate_slope),
('aspect', generate_aspect),
('curvature', generate_curvature),
('svf', generate_svf),
('lrm', generate_lrm),
('pos_open', lambda d, b, v, r: generate_openness(d, b, v, r, positive=True)),
('neg_open', lambda d, b, v, r: generate_openness(d, b, v, r, positive=False)),
('pos_open', lambda d, b, v, r, shared=None: generate_openness(d, b, v, r, positive=True, shared=shared)),
('neg_open', lambda d, b, v, r, shared=None: generate_openness(d, b, v, r, positive=False, shared=shared)),
('mslrm', generate_mslrm),
('tpi', generate_tpi),
('sailore', generate_sailore),
@ -106,16 +106,26 @@ 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, quality=98, only_viz=None, skip_viz=None, output_format='avif'):
self.input_dir = Path(input_dir)
self.output_dir = Path(output_dir)
self.resolution = resolution
# Accept single float or comma-separated string for multi-resolution
if isinstance(resolution, str):
self.resolutions = [float(r.strip()) for r in resolution.split(',')]
elif isinstance(resolution, (list, tuple)):
self.resolutions = [float(r) for r in resolution]
else:
self.resolutions = [float(resolution)]
self.resolution = self.resolutions[0] # Primary resolution (backward compat)
self.workers = workers
self.force = force
self.ground_method = ground_method
self.force_classify = force_classify
self.keep_tif = keep_tif
self.no_viewer = no_viewer
self.quality = quality
self.only_viz = only_viz
self.skip_viz = skip_viz
self.output_format = output_format
self.temp_dir = self.output_dir / "temp"
if not self.input_dir.exists():
@ -126,21 +136,43 @@ class LidarArchaeoPipeline:
self.dtm_dir = self.output_dir / "DTM"
self.vis_dir = self.output_dir / "visualisations"
self.pdf_dir = self.output_dir / "rapports"
for d in [self.dtm_dir, self.vis_dir, self.pdf_dir]:
for d in [self.dtm_dir, self.vis_dir]:
d.mkdir(exist_ok=True)
# Filter visualizations based on --only / --skip
all_viz_names = [name for name, _ in VIZ_STEPS]
if only_viz:
invalid = set(only_viz) - set(all_viz_names)
if invalid:
raise ValueError(f"Visualisations inconnues: {', '.join(invalid)}. Disponibles: {', '.join(all_viz_names)}")
self.viz_steps = [(n, f) for n, f in VIZ_STEPS if n in only_viz]
elif skip_viz:
invalid = set(skip_viz) - set(all_viz_names)
if invalid:
raise ValueError(f"Visualisations inconnues: {', '.join(invalid)}. Disponibles: {', '.join(all_viz_names)}")
self.viz_steps = [(n, f) for n, f in VIZ_STEPS if n not in skip_viz]
else:
self.viz_steps = VIZ_STEPS
logger.info("Pipeline initialisé")
logger.info(f" Entrée : {self.input_dir}")
logger.info(f" Sortie : {self.output_dir}")
logger.info(f" Résolution : {resolution}m/px")
if len(self.resolutions) > 1:
logger.info(f" Résolutions : {', '.join(f'{r}m/px' for r in self.resolutions)}")
else:
logger.info(f" Résolution : {self.resolution}m/px")
logger.info(f" Workers : {workers}")
logger.info(f" Force : {'OUI' if self.force else 'non (skip existing)'}")
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'}")
logger.info(f" Qualité {self.output_format.upper()}: {self.quality if self.quality < 100 else 'lossless'}")
if only_viz:
logger.info(f" Visualisations: uniquement {', '.join(only_viz)}")
elif skip_viz:
logger.info(f" Visualisations: tout sauf {', '.join(skip_viz)}")
logger.info(f" Visualisations: {len(self.viz_steps)}/{len(VIZ_STEPS)}")
def find_laz_files(self):
"""Find all LAZ/LAS files in input directory."""
@ -162,27 +194,76 @@ class LidarArchaeoPipeline:
return False
return True
def generate_all_visualizations(self, dtm_file, basename):
@staticmethod
def _expected_output_path(name, basename, file_vis_dir, output_format='avif'):
"""Return the expected output filename for a visualization step."""
ext = 'avif' if output_format == 'avif' else 'webp'
if name == 'pos_open':
return file_vis_dir / f"{basename}_positive_openness.{ext}"
elif name == 'neg_open':
return file_vis_dir / f"{basename}_negative_openness.{ext}"
elif name == 'hillshade':
return file_vis_dir / f"{basename}_hillshade_multi.{ext}"
else:
return file_vis_dir / f"{basename}_{name}.{ext}"
def generate_all_visualizations(self, dtm_file, basename, resolution=None, vis_dir=None):
"""Generate all archaeological visualizations for one DTM file.
Returns a dict of {name: tif_path} for successful generations.
Optimisation: SharedDEM is only computed if at least one visualization
needs to be generated. When all WebP outputs exist, SharedDEM is
skipped entirely (saves ~2min per file on re-runs).
"""
if resolution is None:
resolution = self.resolution
logger.info(" Génération visualisations:")
# Create per-file subdirectory
file_vis_dir = self.vis_dir / basename
# Use provided vis_dir (for multi-resolution subdirectories) or default
file_vis_dir = vis_dir if vis_dir else (self.vis_dir / basename)
file_vis_dir.mkdir(exist_ok=True)
total = len(self.viz_steps)
# Phase 1: determine which visualizations need generation
needs_generation = {} # name -> True/False
for name, func in self.viz_steps:
if self.force:
needs_generation[name] = True
else:
expected_webp = self._expected_output_path(name, basename, file_vis_dir, self.output_format)
needs_generation[name] = not expected_webp.exists()
to_generate = [n for n, needed in needs_generation.items() if needed]
ign_only = all(name in ('ortho', 'topo') for name in to_generate)
needs_shared = any(name not in ('ortho', 'topo') for name in to_generate)
if not to_generate:
logger.info(" Toutes les visualisations déjà existantes — ignorées")
# Still need to return results dict for PDF check
vis_results = {}
total = len(VIZ_STEPS)
elapsed_times = []
for name, func in self.viz_steps:
vis_results[name] = self._expected_output_path(name, basename, file_vis_dir, self.output_format)
return vis_results
# Phase 2: compute SharedDEM only if needed
shared = None
if needs_shared:
logger.info(" Pré-calcul données partagées (gradient, LRM)...")
t_shared = time.time()
shared = SharedDEM(dtm_file, resolution)
logger.info(f" ✓ Données partagées prêtes ({time.time()-t_shared:.1f}s)")
# Phase 3: generate visualizations
vis_results = {}
for idx, (name, func) in enumerate(self.viz_steps, 1):
if not needs_generation[name]:
logger.info(f" [{idx}/{total}] {name}: déjà existant, ignoré")
vis_results[name] = self._expected_output_path(name, basename, file_vis_dir, self.output_format)
continue
for idx, (name, func) in enumerate(VIZ_STEPS, 1):
# When --force, delete existing TIF to ensure clean regeneration
if self.force:
for tif in file_vis_dir.glob(f"{basename}_{name}.tif"):
tif.unlink(missing_ok=True)
# Special cases for differently-named TIFs
if name == 'pos_open':
for tif in file_vis_dir.glob(f"{basename}_positive_openness.tif"):
tif.unlink(missing_ok=True)
@ -193,40 +274,18 @@ class LidarArchaeoPipeline:
for tif in file_vis_dir.glob(f"{basename}_hillshade_multi.tif"):
tif.unlink(missing_ok=True)
# Check if output WebP already exists (skip unless --force)
if not self.force:
# Determine expected WebP filename from the viz name
# Special cases for openness and IGN overlays
if name == 'pos_open':
expected_webp = file_vis_dir / f"{basename}_positive_openness.webp"
elif name == 'neg_open':
expected_webp = file_vis_dir / f"{basename}_negative_openness.webp"
elif name == 'hillshade':
expected_webp = file_vis_dir / f"{basename}_hillshade_multi.webp"
elif name in ('ortho', 'topo'):
expected_webp = file_vis_dir / f"{basename}_{name}.webp"
else:
expected_webp = file_vis_dir / f"{basename}_{name}.webp"
if expected_webp.exists():
logger.info(f" [{idx}/{total}] {name}: déjà existant, ignoré")
vis_results[name] = expected_webp # Track as existing file
continue
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, resolution)
else:
result = func(dtm_file, basename, file_vis_dir, resolution, shared=shared)
vis_results[name] = result
elapsed = time.time() - t0
elapsed_times.append(elapsed)
if result:
eta = ""
if len(elapsed_times) > 1:
avg_time = sum(elapsed_times) / len(elapsed_times)
remaining = (total - idx) * avg_time
eta = f" — ETA: {remaining:.0f}s"
logger.info(f" [{idx}/{total}] ✓ {name} ({elapsed:.1f}s){eta}")
logger.info(f" [{idx}/{total}] ✓ {name} ({elapsed:.1f}s)")
else:
logger.warning(f" [{idx}/{total}] ✗ {name} — no output ({elapsed:.1f}s)")
except Exception as e:
@ -236,25 +295,43 @@ class LidarArchaeoPipeline:
# Free GPU memory between visualizations to prevent OOM
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)
# Convert to output format (only newly generated TIFs, not skipped ones)
fmt_label = self.output_format.upper()
logger.info(f" Conversion images {fmt_label}:")
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)
if webp_file:
logger.info(f"{webp_file.name}")
img_file = tif_to_png(tif_file, file_vis_dir, resolution, keep_tif=self.keep_tif, source_info=source_info, quality=self.quality, output_format=self.output_format)
if img_file:
logger.info(f"{img_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
@staticmethod
def _res_suffix(resolution):
"""Return suffix for additional resolutions (empty string for primary)."""
if resolution == 0.5:
return "" # Default resolution — no suffix
res_str = f"{resolution}".replace('.', 'p')
return f"_r{res_str}"
def process_file(self, laz_file):
"""Process a single LAZ file through the full pipeline."""
"""Process a single LAZ file through the full pipeline.
If self.resolutions has multiple entries, processes each resolution:
- Primary resolution uses current naming (no suffix)
- Additional resolutions use _r0p2 suffix in directories/filenames
- Ground classification is done once and shared across resolutions
"""
basename = _file_basename(laz_file)
_file_filter.basename = basename
t_start = time.time()
@ -263,8 +340,40 @@ class LidarArchaeoPipeline:
logger.info(f"FICHIER : {basename}")
logger.info("=" * 60)
# Step 1: Ground classification
logger.info("[1/6] Classification du sol...")
# Validate file integrity before any processing
from .dtm import validate_laz
if not validate_laz(laz_file):
return False
# Step 1: Ground classification (shared across all resolutions)
las_file = None
t_classif = 0
for i, res in enumerate(self.resolutions):
res_suffix = self._res_suffix(res)
dtm_path = self.dtm_dir / f"{basename}_dtm{res_suffix}.tif"
if dtm_path.exists():
import rasterio
try:
with rasterio.open(dtm_path) as src:
existing_res = abs(src.transform.a)
if abs(existing_res - res) > 0.01:
logger.info(f" DTM{res_suffix} existant à {existing_res}m/px — résolution demandée {res}m/px → régénération")
dtm_path.unlink()
else:
if i == 0:
logger.info(f"[1/5] Classification du sol — sautée (DTM existant)")
logger.info(f"[2/5] Génération DTM {res}m/px — sautée (DTM existant)")
else:
logger.info(f" DTM {res}m/px déjà existant — ignoré")
continue
except Exception:
logger.warning(f"Impossible de lire le DTM existant — régénération")
dtm_path.unlink()
# Need to classify/generate DTM for this resolution
if las_file is None:
# First time: do 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
@ -273,54 +382,47 @@ class LidarArchaeoPipeline:
return False
logger.info(f" ✓ Classification terminée ({t_classif:.1f}s)")
# Step 2: Generate DTM
logger.info("[2/6] Génération DTM...")
# Generate DTM at this resolution
logger.info(f"{'[2/5]' if i == 0 else ' '} Génération DTM {res}m/px...")
t2 = time.time()
dtm_file = create_dtm_fast(las_file, basename, self.dtm_dir, self.resolution)
dtm_file = create_dtm_fast(las_file, basename, self.dtm_dir, res, force=self.force, output_suffix=res_suffix)
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)")
logger.error(f" ✗ Échec DTM {res}m/px ({t_dtm:.1f}s)")
if i == 0:
return False # Primary resolution failure is fatal
continue # Additional resolution failure is non-fatal
logger.info(f" ✓ DTM {res}m/px terminé ({t_dtm:.1f}s)")
# Step 3: Visualizations
logger.info("[3/6] Visualisations archéologiques...")
self.generate_all_visualizations(dtm_file, basename)
# Process each resolution: visualizations + PDF
all_vis_results = {}
for res in self.resolutions:
res_suffix = self._res_suffix(res)
dtm_path = self.dtm_dir / f"{basename}_dtm{res_suffix}.tif"
# 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)")
if not dtm_path.exists():
logger.warning(f" DTM {res}m/px manquant — visualisations ignorées")
continue
# 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)")
import rasterio
with rasterio.open(dtm_path) as src:
actual_res = abs(src.transform.a)
# 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)")
if len(self.resolutions) > 1:
logger.info(f" --- Résolution {res}m/px ---")
# For additional resolutions, use suffixed subdirectory
if res_suffix:
vis_dir = self.vis_dir / f"{basename}{res_suffix}"
else:
logger.info("[6/6] Viewer web: ignoré (--no-viewer)")
vis_dir = self.vis_dir / basename
vis_dir.mkdir(exist_ok=True)
self.generate_all_visualizations(dtm_path, basename, actual_res, vis_dir=vis_dir)
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")
_file_filter.basename = None
return True
@ -348,11 +450,14 @@ class LidarArchaeoPipeline:
logger.info(f"Traitement parallèle avec {self.workers} workers...")
logger.info(f"Fichiers: {len(files)}")
with ProcessPoolExecutor(max_workers=self.workers) as executor:
# Pass resolutions as comma-separated string for multiprocessing serialization
resolutions_str = ','.join(str(r) for r in self.resolutions)
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), resolutions_str, self.force, self.ground_method, self.force_classify, self.keep_tif, self.quality, self.only_viz, self.skip_viz, self.output_format): laz_file
for laz_file in files
}
done = 0
try:
for future in as_completed(future_to_file):
laz_file = future_to_file[future]
done += 1
@ -365,12 +470,22 @@ class LidarArchaeoPipeline:
logger.error(f" [{done}/{len(files)}] ✗ {laz_file.name}: {e}")
logger.debug(f" Traceback:", exc_info=True)
results[laz_file.name] = False
except KeyboardInterrupt:
logger.info("Interruption — annulation des travaux en cours...")
for f in future_to_file:
f.cancel()
executor.shutdown(wait=False, cancel_futures=True)
logger.info("Travaux annulés.")
return
else:
total = len(files)
for idx, laz_file in enumerate(files, 1):
logger.info(f"--- Fichier {idx}/{total} ---")
try:
results[laz_file.name] = self.process_file(laz_file)
except KeyboardInterrupt:
logger.info("Interruption — arrêt immédiat.")
return
except Exception as e:
logger.error(f"✗ Erreur traitement {laz_file.name}: {e}")
logger.debug("Traceback:", exc_info=True)
@ -384,18 +499,18 @@ class LidarArchaeoPipeline:
logger.info("=" * 60)
logger.info("RÉSUMÉ")
logger.info("=" * 60)
for name, ok in results.items():
status = "" if ok else ""
logger.info(f" {status} {name}")
logger.info("-" * 60)
logger.info(f" Succès : {success_count}/{len(results)}")
if fail_count:
logger.info(f" Échecs : {fail_count}/{len(results)}")
for name, ok in results.items():
if not ok:
logger.info(f"{name}")
logger.info(f" Durée totale : {t_pipeline_total:.1f}s ({t_pipeline_total/60:.1f}min)")
logger.info(f"\nRésultats dans: {self.output_dir}")
logger.info(f" • DTM : {self.dtm_dir}")
logger.info(f" • Visualisations: {self.vis_dir}")
logger.info(f" • Rapports PDF : {self.pdf_dir}")
# Clean up temporary files
logger.info("Nettoyage des fichiers temporaires...")
@ -411,7 +526,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, quality=98, only_viz=None, skip_viz=None, output_format='avif'):
"""Standalone function for multiprocessing — creates its own pipeline instance.
Each worker gets its own temp directory to avoid file conflicts.
@ -419,6 +534,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 +547,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, quality=quality, only_viz=only_viz, skip_viz=skip_viz, output_format=output_format)
basename = _file_basename(laz_file_str)
pipeline.temp_dir = pipeline.output_dir / "temp" / basename
pipeline.temp_dir.mkdir(exist_ok=True)

View File

@ -1,13 +1,14 @@
"""Rendering module: colormap registry, GeoTIFF-to-WebP conversion, and PDF report generation.
"""Rendering module: colormap registry, GeoTIFF-to-image conversion, and PDF report generation.
Contains:
- COLORMAPS: registry mapping filename keywords to (cmap, title, legend, description)
- tif_to_png(): convert a GeoTIFF to a WebP visualization with legend, scale bar, north arrow
- tif_to_png(): convert a GeoTIFF to a WebP/AVIF visualization with legend, scale bar, north arrow
- generate_pdf_report(): generate an A3 PDF report with all visualizations
"""
import logging
import time
from datetime import datetime
from pathlib import Path
import numpy as np
@ -73,18 +74,10 @@ COLORMAPS = {
'description': 'Détecte les ruptures de pente — utile pour bords de terrasses et levées',
'vmin_mode': 'symmetric', 'sym_pct': (5, 95),
},
'svf': {
'cmap': 'viridis',
'title': 'Sky-View Factor (Ray-tracing 16 directions)',
'legend': 'Fraction de ciel visible depuis chaque point\nSombre = Encaissé (fossés, vallées, rues)\nClair = Dégagé (sommets, plateformes, plateaux)',
'description': 'Ray-tracing sur 16 azimuts — élimine l\'éclairage, détecte structures linéaires et enclos',
'vmin_mode': 'percentile', 'vmin_pct': 5,
'vmax_mode': 'percentile', 'vmax_pct': 95,
},
'mslrm': {
'cmap': 'RdBu_r',
'title': 'MSRM - Multi-Scale Relief Model (5 échelles)',
'legend': 'Relief combiné à 5 échelles (5m à 100m)\nRouge = Surélévation (mur, tumulus, levée)\nBleu = Dépression (fossé, douve, fossé)\n\nDifférence avec LRM:\nLRM = 1 échelle (15m)\nMSRM = 5 échelles combinées\nMSRM détecte du micro au macro',
'title': 'MSRM - Multi-Scale Relief Model (échelles adaptatives)',
'legend': 'Relief combiné multi-échelles (adapté à la résolution)\nRouge = Surélévation (mur, tumulus, levée)\nBleu = Dépression (fossé, douve)\n\nDifférence avec LRM:\nLRM = 1 échelle (15m)\nMSRM = échelles combinées pondérées\nMSRM détecte du micro au macro',
'description': 'Combine LRM à 5 échelles — détecte structures de 5m à 100m simultanément',
'vmin_mode': 'symmetric', 'sym_pct': (2, 98),
},
@ -113,8 +106,8 @@ COLORMAPS = {
},
'tpi': {
'cmap': 'BrBG',
'title': 'TPI - Topographic Position Index (2 échelles)',
'legend': 'Position dans le paysage\nBrun/Sombre = Plus bas que le voisinage (fossé, vallée)\nVert/Clair = Plus haut que le voisinage (crête, plateau)\nCombine échelle fine 5m + large 100m',
'title': 'TPI - Topographic Position Index (4 échelles)',
'legend': 'Position dans le paysage\nBrun/Sombre = Plus bas que le voisinage (fossé, vallée)\nVert/Clair = Plus haut que le voisinage (crête, plateau)\nCombine 4 échelles : 3m, 15m, 50m, 200m',
'description': 'Identifie la position topographique — utile pour repérer crêtes vs vallées à grande échelle',
'vmin_mode': 'symmetric', 'sym_pct': (2, 98),
},
@ -127,23 +120,23 @@ COLORMAPS = {
},
'roughness': {
'cmap': 'magma',
'title': 'Rugosité de Surface (Écart-type local 5m)',
'legend': 'Irrégularité du terrain dans un voisinage de 5m\nSombre = Surface lisse (route, mur, sol plat)\nClair = Surface rugueuse (végétation, ruines, pierres)\nMax: {vmax:.2f}m',
'title': 'Rugosité Multi-Échelle (3m + 15m)',
'legend': 'Irrégularité du terrain combinée fine + large\nSombre = Surface lisse (route, mur, sol plat)\nClair = Surface rugueuse (végétation, ruines, pierres)\nCombine rugosité fine 3m (70%) + large 15m (30%)',
'description': 'Mesure la variabilité locale — surfaces anthropiques lisses vs naturelles rugueuses',
'vmin_mode': 'fixed', 'vmin_val': 0,
'vmax_mode': 'percentile', 'vmax_pct': 97,
},
'anomalies': {
'cmap': 'coolwarm',
'title': 'Anomalies Statistiques (Z-score x Moran\'s I)',
'legend': 'Anomalies topographiques significatives\nRouge vif = Surélévation anormale (mur, tumulus)\nBleu vif = Dépression anormale (fossé, doline)\nBlanc/gris = Normal\n\nCombine z-score (intensité) et\nMoran\'s I (regroupement spatial)',
'title': 'Anomalies Statistiques (MSRM multi-échelle + Moran\'s I)',
'legend': 'Anomalies topographiques significatives\nRouge vif = Surélévation anormale (mur, tumulus)\nBleu vif = Dépression anormale (fossé, doline)\nBlanc/gris = Normal\n\nCombine MSRM normalisé (intensité) et\nMoran\'s I (regroupement spatial)',
'description': 'Détecte uniquement les anomalies statistiquement significatives — filtre le bruit de fond',
'vmin_mode': 'symmetric', 'sym_pct': (5, 95),
},
'wavelet': {
'cmap': 'cividis',
'title': 'Ondelette Mexican Hat (CWT multi-échelle)',
'legend': 'Réponse de la transformée en ondelette à 5 échelles\nÉchelles: 2m, 5m, 10m, 20m, 50m\n\nClair = Structure détectée à cette échelle\nSombre = Pas de structure\n\nOptimisé pour formes circulaires:\ntumulus, enclos, fossés annulaires',
'legend': 'Réponse de la transformée en ondelette\nÉchelles adaptées à la résolution\n\nClair = Structure détectée à cette échelle\nSombre = Pas de structure\n\nOptimisé pour formes circulaires:\ntumulus, enclos, fossés annulaires',
'description': 'Transformée en ondelette 2D — excellente pour détecter structures circulaires',
'vmin_mode': 'symmetric', 'sym_pct': (2, 98),
},
@ -238,26 +231,54 @@ 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):
"""Convert GeoTIFF to visualization WebP with GPS coordinates, legend, and scale bar.
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, quality=98, output_format='avif'):
"""Convert GeoTIFF to visualization image (WebP or AVIF) with GPS coordinates, legend, and scale bar.
Args:
tif_file: Path to input GeoTIFF.
vis_dir: Output directory for the WebP file.
vis_dir: Output directory for the image file.
resolution: Grid resolution in m/px.
keep_tif: If True, keep the source TIFF after conversion.
source_info: Dict with method/date/basename for metadata.
quality: Image quality (1-100). Use 100 for lossless. Default 95.
output_format: Output format ('webp' or 'avif'). Default 'webp'.
Returns:
Path to output WebP file, or None on failure.
Path to output image file, or None on failure.
"""
if not tif_file or not tif_file.exists():
return None
webp_file = vis_dir / f"{tif_file.stem}.webp"
ext = 'avif' if output_format == 'avif' else 'webp'
output_file = vis_dir / f"{tif_file.stem}.{ext}"
try:
with rasterio.open(tif_file) as src:
is_rgb = src.count >= 3 and ('ortho' in str(tif_file) or 'topo' in str(tif_file))
is_rgb = src.count >= 3 and any(k in str(tif_file) for k in ('ortho', 'topo'))
if is_rgb:
data = src.read([1, 2, 3])
@ -362,9 +383,9 @@ def tif_to_png(tif_file, vis_dir, resolution, keep_tif=False):
# Fixed data area position — identical for ALL visualization types
# This ensures overlay/superposition works across all WebP images
data_left = 0.08
data_bottom = 0.12
data_bottom = 0.19
data_width_frac = 0.74
data_height_frac = 0.78
data_height_frac = 0.71
ax = fig.add_axes([data_left, data_bottom, data_width_frac, data_height_frac])
if is_rgba or is_rgb:
@ -440,20 +461,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,46 +498,77 @@ 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,
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(bar_left + bar_width_frac / 2, bar_bottom + bar_height + 0.12,
f"{scale_m} m", ha='center', va='bottom', fontsize=9, fontweight='bold',
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')
# Save as PNG then convert to WebP — fixed layout, no bbox_inches='tight'
# Save as PNG then convert to final format — fixed layout, no bbox_inches='tight'
save_dpi = 200 if width > 3000 else 150
png_temp = vis_dir / f"{tif_file.stem}_temp.png"
try:
@ -513,159 +577,24 @@ def tif_to_png(tif_file, vis_dir, resolution, keep_tif=False):
plt.close()
img = PILImage.open(str(png_temp))
img.save(str(webp_file), format='WEBP', lossless=True)
pil_format = 'AVIF' if output_format == 'avif' else 'WEBP'
if quality >= 100:
img.save(str(output_file), format=pil_format, lossless=True)
else:
img.save(str(output_file), format=pil_format, quality=quality)
png_temp.unlink(missing_ok=True)
# Delete source TIFF (unless --keep-tif)
if not keep_tif:
tif_file.unlink(missing_ok=True)
return webp_file
return output_file
except Exception as e:
logger.error(f" Erreur conversion WebP: {e}", exc_info=True)
logger.error(f" Erreur conversion {ext.upper()}: {e}", exc_info=True)
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.

View File

@ -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(
'<html><head><title>LiDAR Server</title>'
'<style>body{font-family:sans-serif;max-width:600px;margin:40px auto;padding:0 20px}'
'h1{color:#1a1a2e}a{color:#4fc3f7}</style></head>'
'<body>'
'<h1>Serveur LiDAR Archéologique</h1>'
'<p><a href="/viewer">Visualisations</a></p>'
'<p>TiTiler API: <a href="/cog/">/cog/</a></p>'
'</body></html>'
)
@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'<li><a href="/viewer/{f.stem}">{f.stem}</a></li>'
for f in viewers
)
return HTMLResponse(
f'<html><head><title>LiDAR Viewer</title>'
f'<style>body{{font-family:sans-serif;max-width:600px;margin:40px auto;padding:0 20px}}'
f'h1{{color:#1a1a2e}}li{{margin:8px 0}}a{{color:#4fc3f7}}</style></head>'
f'<body><h1>Visualisations LiDAR</h1><ul>{links}</ul></body></html>'
)
return HTMLResponse('<html><body><p>Aucun viewer disponible</p></body></html>')
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()

View File

@ -199,3 +199,34 @@ class TestFlow:
# log1p(x) >= 0 for x >= 0
valid = data[~np.isnan(data)]
assert np.nanmin(valid) >= 0
class TestLocalDominance:
def test_generates_tif(self, synthetic_dem, tmp_output_dir):
from lidar_pipeline.visualizations import generate_local_dominance
result = generate_local_dominance(synthetic_dem, "test", tmp_output_dir, 5.0)
assert result is not None
assert result.exists()
assert result.suffix == ".tif"
def test_dominance_values_0_1(self, synthetic_dem, tmp_output_dir):
import rasterio
from lidar_pipeline.visualizations import generate_local_dominance
result = generate_local_dominance(synthetic_dem, "test", tmp_output_dir, 5.0)
with rasterio.open(result) as src:
data = src.read(1)
valid = data[~np.isnan(data)]
assert np.nanmin(valid) >= 0, "Local dominance should be >= 0"
assert np.nanmax(valid) <= 1, "Local dominance should be <= 1"
def test_dominance_nan_mask_preserved(self, synthetic_dem, tmp_output_dir):
"""Check that NaN zones from original DEM are preserved."""
import rasterio
from lidar_pipeline.visualizations import generate_local_dominance
result = generate_local_dominance(synthetic_dem, "test", tmp_output_dir, 5.0)
# The synthetic DEM has no NaN, so this just verifies the output is valid
with rasterio.open(result) as src:
data = src.read(1)
# Shape should match input
assert data.shape[0] > 0
assert data.shape[1] > 0

View File

@ -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' <div class="layer-control">\n'
f' <label>\n'
f' <input type="checkbox" data-layer="{layer["name"]}" {checked}>\n'
f' <span class="layer-title">{layer["title"]}</span>\n'
f' <span class="layer-type">{layer_type}</span>\n'
f' </label>\n'
f' <input type="range" class="opacity-slider" data-layer="{layer["name"]}" '
f'min="0" max="100" value="{initial_opacity}">\n'
f' </div>'
)
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'''<!DOCTYPE html>
<html lang="fr">
<head>
<meta charset="utf-8">
<meta name="viewport" content="width=device-width, initial-scale=1">
<title>LiDAR Archéologique — {basename}</title>
<link rel="stylesheet" href="https://unpkg.com/maplibre-gl@4.7.1/dist/maplibre-gl.css">
<script src="https://unpkg.com/maplibre-gl@4.7.1/dist/maplibre-gl.js"></script>
<style>
* {{ margin: 0; padding: 0; box-sizing: border-box; }}
body {{ font-family: -apple-system, BlinkMacSystemFont, 'Segoe UI', Roboto, sans-serif; }}
#map {{ position: absolute; top: 0; left: 280px; right: 0; bottom: 0; }}
#panel {{
position: absolute; top: 0; left: 0; width: 280px; bottom: 0;
background: #1a1a2e; color: #e0e0e0; overflow-y: auto;
z-index: 10; box-shadow: 2px 0 8px rgba(0,0,0,0.3);
}}
#panel h1 {{
font-size: 14px; padding: 12px 14px; background: #16213e;
border-bottom: 1px solid #0f3460; letter-spacing: 0.5px;
}}
#panel h2 {{
font-size: 11px; padding: 8px 14px; color: #a0a0a0;
text-transform: uppercase; letter-spacing: 1px;
border-bottom: 1px solid #2a2a4a;
}}
.layer-group {{ padding: 4px 0; border-bottom: 1px solid #2a2a4a; }}
.layer-control {{
padding: 6px 14px;
transition: background 0.15s;
}}
.layer-control:hover {{ background: #252550; }}
.layer-control label {{
display: flex; align-items: center; gap: 8px;
font-size: 13px; cursor: pointer; margin-bottom: 4px;
}}
.layer-control input[type="checkbox"] {{
width: 16px; height: 16px; accent-color: #4fc3f7;
}}
.layer-control .layer-title {{ flex: 1; }}
.layer-control .layer-type {{
font-size: 10px; color: #888; text-transform: uppercase;
}}
.opacity-slider {{
width: 100%; height: 4px; margin: 2px 0 0 24px;
-webkit-appearance: none; appearance: none;
background: #333; border-radius: 2px; outline: none;
}}
.opacity-slider::-webkit-slider-thumb {{
-webkit-appearance: none; appearance: none;
width: 14px; height: 14px; border-radius: 50%;
background: #4fc3f7; cursor: pointer;
}}
.basemap-section {{
padding: 10px 14px; border-bottom: 1px solid #2a2a4a;
}}
.basemap-section h3 {{
font-size: 11px; color: #a0a0a0; text-transform: uppercase;
letter-spacing: 1px; margin-bottom: 8px;
}}
.basemap-btn {{
display: inline-block; padding: 5px 10px; margin: 2px;
background: #2a2a4a; color: #ccc; border: 1px solid #444;
border-radius: 4px; cursor: pointer; font-size: 12px;
transition: all 0.15s;
}}
.basemap-btn:hover {{ background: #3a3a6a; }}
.basemap-btn.active {{ background: #4fc3f7; color: #000; border-color: #4fc3f7; }}
#coords {{
position: absolute; bottom: 8px; left: 288px;
background: rgba(26,26,46,0.9); color: #4fc3f7;
padding: 4px 10px; border-radius: 4px; font-size: 12px;
font-family: monospace; z-index: 5;
}}
</style>
</head>
<body>
<div id="panel">
<h1>LiDAR Archéologique</h1>
<div style="padding: 6px 14px; font-size: 12px; color: #888;">
{basename}<br>
<span style="font-size: 10px;">Res: {resolution}m/px | EPSG:2154</span>
</div>
<div class="basemap-section">
<h3>Fond de carte</h3>
<button class="basemap-btn active" onclick="setBasemap('osm')">OSM</button>
<button class="basemap-btn" onclick="setBasemap('ign-photo')">Photo IGN</button>
<button class="basemap-btn" onclick="setBasemap('ign-map')">Carte IGN</button>
</div>
<h2>Visualisations</h2>
<div class="layer-group" id="layers">
{layer_controls_html}
</div>
</div>
<div id="map"></div>
<div id="coords"></div>
<script>
const TITILER_URL = '{titiler_url}';
const BASENAME = '{basename}';
const BOUNDS = {bounds_json};
const LAYERS = {layers_json};
const COG_PATH_PREFIX = '{cog_path_prefix}';
const LAYER_COLORMAPS = {colormaps_json};''')
# JavaScript code — use regular string to avoid f-string issues with {z}/{x}/{y}
html_parts.append(r'''
// Map initialization
const map = new maplibregl.Map({
container: 'map',
style: {
version: 8,
sources: {},
layers: []
},
center: [' + f'{center_lng}, {center_lat}' + r'],
zoom: 13,
maxZoom: 19,
minZoom: 8,
});
// Basemap layers
const basemapStyles = {
'osm': {
type: 'raster',
tiles: ['https://tile.openstreetmap.org/{z}/{x}/{y}.png'],
tileSize: 256,
attribution: '&copy; OpenStreetMap',
maxzoom: 19
},
'ign-photo': {
type: 'raster',
tiles: ['https://data.geopf.fr/wmts?SERVICE=WMTS&VERSION=1.0.0&REQUEST=GetTile&LAYER=ORTHOIMAGERY.ORTHOPHOTOS&STYLE=normal&TILEMATRIXSET=PM&TILEMATRIX={z}&TILECOL={x}&TILEROW={y}&FORMAT=image/jpeg'],
tileSize: 256,
attribution: '&copy; IGN',
maxzoom: 20
},
'ign-map': {
type: 'raster',
tiles: ['https://data.geopf.fr/wmts?SERVICE=WMTS&VERSION=1.0.0&REQUEST=GetTile&LAYER=GEOGRAPHICALGRIDSYSTEMS.PLANIGNV2&STYLE=normal&TILEMATRIXSET=PM&TILEMATRIX={z}&TILECOL={x}&TILEROW={y}&FORMAT=image/png'],
tileSize: 256,
attribution: '&copy; IGN',
maxzoom: 19
}
};
let currentBasemap = 'osm';
function setBasemap(name) {
if (map.getLayer('basemap')) map.removeLayer('basemap');
if (map.getSource('basemap')) map.removeSource('basemap');
currentBasemap = name;
const style = basemapStyles[name];
map.addSource('basemap', style);
map.addLayer({ id: 'basemap', type: 'raster', source: 'basemap' });
document.querySelectorAll('.basemap-btn').forEach(btn => {
const names = {'osm': 'OSM', 'ign-photo': 'Photo IGN', 'ign-map': 'Carte IGN'};
btn.classList.toggle('active', btn.textContent.trim() === names[name]);
});
// Re-add data layers on top of basemap
LAYERS.forEach(layer => {
const layerId = 'layer-' + layer.name;
if (map.getLayer(layerId)) {
map.moveLayer(layerId, 'basemap');
}
});
}
// Add visualization layers
function addVisualizationLayer(layerInfo) {
const name = layerInfo.name;
const isRgb = layerInfo.is_rgb;
const cogUrl = COG_PATH_PREFIX + '/' + layerInfo.file;
let tileUrl;
if (isRgb) {
tileUrl = TITILER_URL + '/cog/tiles/{z}/{x}/{y}.png?url=' +
encodeURIComponent(cogUrl) + '&bidx=1&bidx=2&bidx=3';
} else {
const cmap = LAYER_COLORMAPS[name] || 'viridis';
tileUrl = TITILER_URL + '/cog/tiles/{z}/{x}/{y}.png?url=' +
encodeURIComponent(cogUrl) + '&colormap_name=' + cmap + '&rescale=-1,1';
}
const sourceId = 'source-' + name;
const layerId = 'layer-' + name;
map.addSource(sourceId, {
type: 'raster',
tiles: [tileUrl],
tileSize: 256,
minzoom: 8,
maxzoom: 19,
});
map.addLayer({
id: layerId,
type: 'raster',
source: sourceId,
paint: {
'raster-opacity': name === 'hillshade_multi' ? 1.0 : 0.0,
'raster-resampling': 'bilinear',
}
});
}
// Layer control event handlers
document.querySelectorAll('.layer-control input[type="checkbox"]').forEach(cb => {
cb.addEventListener('change', () => {
const layerId = 'layer-' + cb.dataset.layer;
if (map.getLayer(layerId)) {
map.setLayoutProperty(layerId, 'visibility',
cb.checked ? 'visible' : 'none');
}
// When checking a layer, set its opacity to 70% if it was 0
const slider = document.querySelector('.opacity-slider[data-layer="' + cb.dataset.layer + '"]');
if (cb.checked && slider && parseInt(slider.value) === 0) {
slider.value = 70;
const layerId = 'layer-' + cb.dataset.layer;
if (map.getLayer(layerId)) {
map.setPaintProperty(layerId, 'raster-opacity', 0.7);
}}
});
});
document.querySelectorAll('.opacity-slider').forEach(slider => {
slider.addEventListener('input', () => {
const layerId = 'layer-' + slider.dataset.layer;
if (map.getLayer(layerId)) {
map.setPaintProperty(layerId, 'raster-opacity', slider.value / 100);
}
});
});
// Coordinate display
map.on('mousemove', (e) => {
const { lng, lat } = e.lngLat;
document.getElementById('coords').textContent =
lat.toFixed(6) + 'N ' + lng.toFixed(6) + 'E';
});
// Fit map to bounds on load
map.on('load', () => {
setBasemap('osm');
// Add all visualization layers
LAYERS.forEach(layer => {
addVisualizationLayer(layer);
});
// Initially hide all except hillshade
LAYERS.forEach(layer => {
if (layer.name !== 'hillshade_multi') {
const layerId = 'layer-' + layer.name;
if (map.getLayer(layerId)) {
map.setLayoutProperty(layerId, 'visibility', 'none');
}
}
});
// Fit bounds
map.fitBounds([[BOUNDS.west, BOUNDS.south], [BOUNDS.east, BOUNDS.north]], { padding: 20 });
});
</script>
</body>
</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

File diff suppressed because it is too large Load Diff

119
run.sh
View File

@ -3,21 +3,22 @@
# Utilisation: ./run.sh [options]
#
# Options:
# -r RESOLUTION Résolution en m/px (défaut: 0.5)
# -r RESOLUTION Résolution en m/px, ou multiples séparées par virgules (défaut: 0.5, ex: 0.5,0.2)
# -w WORKERS Nombre de workers parallèles (défaut: 1)
# -g Activer l'accélération GPU
# -v Mode verbeux (timestamps + niveaux)
# --debug Mode debug (détails internes fichier:ligne)
# -f / --force Régénérer tous les fichiers même si existants
# --keep-tif Conserver les fichiers TIFF intermédiaires
# --no-viewer Ne pas générer le viewer web (COGs + HTML)
# -v Mode verbeux
# --debug Mode debug
# -f / --force Régénérer tous les fichiers
# --keep-tif Conserver les fichiers TIFF
# --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
# --ground-classification {auto,smrf,csf}
# --quality N Qualité image 1-100 (défaut: 98)
# --lossless Compression lossless
# --format FMT Format de sortie : avif (défaut) ou webp
# --only VIZ... Générer uniquement ces visualisations
# --skip VIZ... Exclure ces visualisations
# --file NOM... Traiter un ou plusieurs fichiers LAZ
# --test Exécuter les tests unitaires
# serve Démarrer le serveur cartographique
# -h Afficher l'aide complète
set -e
@ -27,36 +28,14 @@ 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)"
echo " -r RESOLUTION Résolution en m/px, ou multiples (défaut: 0.5, ex: 0.5,0.2)"
echo " -w WORKERS Nombre de workers CPU parallèles (défaut: 1)"
echo " -g Activer l'accélération GPU NVIDIA"
echo " -v Mode verbeux (timestamps + niveaux)"
@ -64,13 +43,11 @@ if [ $# -eq 0 ]; then
echo " -f / --force Régénérer tous les fichiers même si les WebP existent"
echo " --force-classification"
echo " Reclassifier le sol même si le fichier .las existe"
echo " --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 " --keep-tif Conserver les TIFF pour régénérer les WebP"
echo " --ground-classification {auto,smrf,csf}"
echo " Méthode de classification du sol (défaut: auto)"
echo " --file NOM... Traiter un ou plusieurs fichiers LAZ (nom complet sans .laz/.las)"
echo " --test Exécuter les tests unitaires"
echo " serve Démarrer le serveur cartographique"
echo " -h Afficher cette aide"
echo ""
echo "Exemples:"
@ -78,11 +55,11 @@ if [ $# -eq 0 ]; then
echo " $0 -g -w 4 # GPU + 4 workers"
echo " $0 -g -v # GPU + verbeux"
echo " $0 -g -r 0.2 # Haute résolution"
echo " $0 -g --force # Tout régénérer (WebP + classification)"
echo " $0 -g -r 0.5,0.2 # Multi-résolution (0.5m + 0.2m)"
echo " $0 -g --force # Régénérer WebP (DTM conservé si --keep-tif)"
echo " $0 -g --force-classification # Reclassifier le sol seulement"
echo " $0 -g --ground-classification pmf # Forcer PMF"
echo " $0 -g --ground-classification csf # Forcer CSF (terrain complexe)"
echo " $0 -g --file LHD_...IGN69.copc # Un fichier"
echo " $0 serve # Démarrer le serveur web"
exit 0
fi
@ -95,7 +72,10 @@ FILE_ARGS=""
GROUND_METHOD=""
FORCE_CLASSIFY_FLAG=""
KEEP_TIF_FLAG=""
NO_VIEWER_FLAG=""
QUALITY=""
FORMAT_FLAG=""
ONLY_FLAG=""
SKIP_FLAG=""
# Parse arguments manually (more robust than getopts for mixed short/long options)
while [ $# -gt 0 ]; do
@ -109,16 +89,19 @@ 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 ;;
--quality) QUALITY="--quality $2"; shift 2 ;;
--lossless) QUALITY="--lossless"; shift ;;
--format) FORMAT_FLAG="--format $2"; shift 2 ;;
--only) shift; ONLY_FLAG="--only"; while [ $# -gt 0 ] && [[ ! "$1" =~ ^- ]]; do ONLY_FLAG="$ONLY_FLAG $1"; shift; done ;;
--skip) shift; SKIP_FLAG="--skip"; while [ $# -gt 0 ] && [[ ! "$1" =~ ^- ]]; do SKIP_FLAG="$SKIP_FLAG $1"; shift; done ;;
--file) shift; while [ $# -gt 0 ] && [[ ! "$1" =~ ^- ]]; do FILE_ARGS="$FILE_ARGS $1"; shift; done ;;
--test) ;; # Handled below
-h|--help|-help)
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)"
@ -129,25 +112,35 @@ while [ $# -gt 0 ]; do
echo " -f / --force Régénérer tous les fichiers même si les WebP existent"
echo " --force-classification"
echo " Reclassifier le sol même si le fichier .las existe"
echo " --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 " --keep-tif Conserver les TIFF pour régénérer les WebP"
echo " --ground-classification {auto,smrf,csf}"
echo " Méthode de classification du sol (défaut: auto)"
echo " --file NOM... Traiter un ou plusieurs fichiers LAZ (nom complet sans .laz/.las)"
echo " --quality N Qualité image 1-100 (défaut: 98, 100=lossless)"
echo " --lossless Compression lossless (équivalent à --quality 100)"
echo " --format FMT Format de sortie : avif (défaut) ou webp"
echo " --only VIZ... Générer uniquement ces visualisations"
echo " --skip VIZ... Exclure ces visualisations"
echo " --file NOM... Traiter un ou plusieurs fichiers LAZ"
echo " --test Exécuter les tests unitaires"
echo " serve Démarrer le serveur cartographique"
echo " -h Afficher cette aide"
echo ""
echo "Visualisations disponibles:"
echo " hillshade slope aspect curvature lrm pos_open neg_open"
echo " mslrm tpi sailore roughness anomalies wavelet flow ortho topo"
echo ""
echo "Exemples:"
echo " $0 -g # GPU, auto"
echo " $0 -g -w 4 # GPU + 4 workers"
echo " $0 -g -v # GPU + verbeux"
echo " $0 -g -r 0.2 # Haute résolution"
echo " $0 -g --force # Tout régénérer (WebP + classification)"
echo " $0 -g --force-classification # Reclassifier le sol seulement"
echo " $0 -g --ground-classification pmf # Forcer PMF"
echo " $0 -g -r 0.5,0.2 # Multi-résolution (0.5m + 0.2m)"
echo " $0 -g --force # Régénérer WebP"
echo " $0 -g --only hillshade svf lrm # Seulement 3 visualisations"
echo " $0 -g --skip ortho topo # Sans les overlays IGN"
echo " $0 -g --lossless # WebP lossless"
echo " $0 -g --quality 90 # WebP qualité 90"
echo " $0 -g --ground-classification csf # Forcer CSF (terrain complexe)"
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 ;;
@ -167,7 +160,7 @@ if [[ " $* " == *" --test "* ]]; then
echo "============================================"
echo " Tests unitaires LiDAR Pipeline"
echo "============================================"
docker run --rm $GPU_FLAG \
docker run --rm --init $GPU_FLAG \
"$IMAGE_NAME" \
python3 -m pytest -v --pyargs lidar_pipeline.tests
exit $?
@ -193,22 +186,34 @@ 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 " Qualité image : $([ -n "$QUALITY" ] && echo "$QUALITY" || echo '98')"
echo " Format : $([ -n "$FORMAT_FLAG" ] && echo "${FORMAT_FLAG#--format }" || echo 'avif')"
echo " Classification sol : $([ -n "$GROUND_METHOD" ] && echo "$GROUND_METHOD" || echo 'auto')"
if [ -n "$ONLY_FLAG" ]; then
echo " Visualisations: uniquement${ONLY_FLAG#--only}"
elif [ -n "$SKIP_FLAG" ]; then
echo " Visualisations: tout sauf${SKIP_FLAG#--skip}"
fi
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 $QUALITY $FORMAT_FLAG"
if [ -n "$GROUND_METHOD" ]; then
CMD_ARGS="$CMD_ARGS --ground-classification $GROUND_METHOD"
fi
if [ -n "$ONLY_FLAG" ]; then
CMD_ARGS="$CMD_ARGS $ONLY_FLAG"
fi
if [ -n "$SKIP_FLAG" ]; then
CMD_ARGS="$CMD_ARGS $SKIP_FLAG"
fi
if [ -n "$FILE_ARGS" ]; then
CMD_ARGS="$CMD_ARGS --file $FILE_ARGS"
fi
docker run --rm $GPU_FLAG \
docker run --rm --init $GPU_FLAG \
--user 1000:1000 \
-v "${INPUT_DIR}:/data/input:ro" \
-v "${OUTPUT_DIR}:/data/output" \