Add web map with zone generation API, side job queue, and restore historical DTM hole rendering

- webapp.py: FastAPI serving the continuous map (port 8973) with
  /api/preview, /api/generate and /api/status; tiles are downloaded
  from IGN and processed in a logged subprocess, tracked live in a
  side "File de génération" panel that survives page reloads
- fetch_ign.py: download missing 1 km LiDAR HD tiles from the IGN
  geoplateforme before processing
- index.py: tile thumbnails and 500 m subtiles are now invalidated by
  mtime so regenerating a tile refreshes its cached images; progress
  logging per tile
- dtm.py: back to the historical gap handling (small gaps filled by
  fillnodata only, larger holes left as nodata rendered black);
  lowest-return floor only via --bare-earth, IGN class selection via
  --ign-classes
- cli.py: positional input now optional (--rebuild-index works alone)
- docker-compose.yml: serve (GPU, port 8973) and process services;
  launch via docker compose only (documented in AGENTS.md/AGENTS.md)
- tests: 131 passing, incl. regressions for thumbnail staleness,
  --rebuild-index without input, and nodata rendering
This commit is contained in:
Antoine Jacquin
2026-08-31 18:07:14 +02:00
parent 35bd827790
commit 8ca65155db
19 changed files with 3676 additions and 771 deletions

View File

@ -8,6 +8,9 @@
- lint: not configured - lint: not configured
- format: not configured - format: not configured
- after every edit: `./run.sh --test` - after every edit: `./run.sh --test`
- **RÈGLE 1 — toujours lancer via docker compose** (jamais `docker run` direct) : carte/API → `docker compose up -d serve` (port 8973) ; traitement ponctuel → `docker compose run --rm process [options]` ; logs → `docker compose logs -f serve` ; arrêt → `docker compose down`.
- **RÈGLE 2 — après chaque édition de code : rebuild de l'image puis relance du conteneur.** Le code est baké dans l'image (jamais monté) : sans `docker compose build` suivi d'un `docker compose up -d serve` (recrée le conteneur), l'ANCIEN code continue de tourner. Toujours reconstruire avant de faire tester/valider une modif par l'utilisateur.
- test rapide sans rebuild (code monté par-dessus l'image): `docker run --rm -e PYTHONPATH=/app -v $(pwd)/lidar_pipeline:/app/lidar_pipeline lidar-lidar python3 -m pytest --pyargs lidar_pipeline.tests`
- debug: `./run.sh --debug` (file:line logging); container shell: `docker run --rm -it -v $(pwd)/input:/data/input -v $(pwd)/output:/data/output --entrypoint bash lidar-lidar` - debug: `./run.sh --debug` (file:line logging); container shell: `docker run --rm -it -v $(pwd)/input:/data/input -v $(pwd)/output:/data/output --entrypoint bash lidar-lidar`
## Conventions ## Conventions
@ -21,6 +24,7 @@
- **Default output is AVIF**, not WebP. Use `--format webp` for WebP. Quality default is 98. - **Default output is AVIF**, not WebP. Use `--format webp` for WebP. Quality default is 98.
- **Tests use lazy imports inside each test function**, never at module top, to avoid importing CuPy/GDAL at import time. - **Tests use lazy imports inside each test function**, never at module top, to avoid importing CuPy/GDAL at import time.
- **`_`-prefixed names are critical private**: `_create_ground_pipeline`, `_fallback_to_smrf`, `_fill_nans`, `_init_gpu`, `_process_file_standalone` — do not call from outside their module. - **`_`-prefixed names are critical private**: `_create_ground_pipeline`, `_fallback_to_smrf`, `_fill_nans`, `_init_gpu`, `_process_file_standalone` — do not call from outside their module.
- **`build_index()` writes 3 files**: `output/index.html` (data shell, `const TILES` embedded), `output/assets/app.css` and `output/assets/app.js` (source: `_APP_CSS`/`_APP_JS` constants in `index.py`). `webapp.py` serves `/assets` with no-cache headers. Each tile carries `meta` — ground method read from `DTM/*_dtm{_rXpY}_method.txt` (falls back to the primary-resolution sidecar) + per-viz dates/sizes.
## Commit & Pull Request Guidelines ## Commit & Pull Request Guidelines

View File

@ -1,38 +1,61 @@
version: '3.8' # Lancement du pipeline LiDAR — TOUJOURS via docker compose :
# docker compose build # après chaque édition de code (code baké dans l'image)
# docker compose up -d serve # carte interactive + API sur http://localhost:8973
# docker compose logs -f serve # journal du serveur
# docker compose down # arrêt
# Traitement ponctuel (sans serveur) :
# docker compose run --rm process [-r 0.5,0.2 | --force | --file ...]
# Jupyter (opt-in) :
# docker compose --profile interactive up -d jupyter
services: services:
lidar: # Carte interactive + API de génération (mode ./run.sh --serve)
serve:
build: . build: .
container_name: lidar-archeo image: lidar-lidar
container_name: lidar-serve
init: true
user: "1000:1000" user: "1000:1000"
gpus: all
ports:
- "8973:8973"
volumes: volumes:
# Mount your LAZ files directory here # input/ en écriture : l'API y télécharge les dalles IGN manquantes
- ./input:/data/input:ro - ./input:/data/input
# Output directory
- ./output:/data/output - ./output:/data/output
# Optional: Mount a large data directory
# - /path/to/your/laz/files:/data/input:ro
environment: environment:
- TZ=Europe/Paris - TZ=Europe/Paris
# Processing parameters - LIDAR_INPUT_DIR=/data/input
- RESOLUTION=0.5 - LIDAR_OUTPUT_DIR=/data/output
- WHITEBOX_THREADS=4 # Les générations lancées depuis la carte utilisent le GPU
# Resource limits (adjust based on your system) - LIDAR_GPU=1
deploy: - LIDAR_WORKERS=2
resources: command: python3 -m uvicorn lidar_pipeline.webapp:app --host 0.0.0.0 --port 8973
limits: restart: unless-stopped
cpus: '4'
memory: 8G
reservations:
cpus: '2'
memory: 4G
# Override default command
command: ["process_lidar.py", "/data/input", "-o", "/data/output", "-r", "0.5"]
# Optional: Jupyter notebook for interactive exploration # Traitement ponctuel des dalles input/ (une passe puis arrêt)
process:
build: .
image: lidar-lidar
container_name: lidar-process
init: true
user: "1000:1000"
gpus: all
volumes:
- ./input:/data/input
- ./output:/data/output
environment:
- TZ=Europe/Paris
command: ["python3", "-m", "lidar_pipeline", "/data/input", "-o", "/data/output", "-r", "0.5,0.2", "-g", "all"]
profiles:
- process
# Exploration interactive (opt-in : --profile interactive)
jupyter: jupyter:
build: . build: .
image: lidar-lidar
container_name: lidar-jupyter container_name: lidar-jupyter
init: true
ports: ports:
- "8888:8888" - "8888:8888"
volumes: volumes:

View File

@ -0,0 +1,141 @@
# Classification du sol — options et références
Note de synthèse pour choisir/améliorer l'algorithme de détection du sol.
Contexte : tiles LiDAR HD IGN (Lambert 93), zones de relief fort / rochers /
forêt dense où le sol est sous-classifié et le MNT présente de grands trous.
Tile de référence : `LHD_FXX_0999_6778_PTS_LAMB93_IGN69`
- 65 047 919 points, ~1 km², résolutions 0.5 m et 0.2 m.
- Répartition classes (pré-classification fournisseur) :
classe 2 (sol) **25.79 %**, classe 5 (veg haute) **63.2 %**,
classe 3 (veg basse) 7.85 %, classe 4 (veg moyenne) 2.66 %,
classe 1 0.47 %, classe 6 0.01 %. Aucun point en classe 0.
- Trous dans le MNT existant (avant correction) : **43.5 %** à 0.5 m,
**44.1 %** à 0.2 m (surface du tile, bornes du header).
## Benchmark mesuré (PDAL, 1 km²)
| Méthode | Temps | Points sol | Surface sol* | Trous* |
|---|---|---|---|---|
| **IGN** (pré-classif.) | **9.4 s** | 25.8 % | 84.6 % | 15.4 % |
| **SMRF** | 326.3 s | 36.1 % | 90.4 % | 9.6 % |
| **CSF** | 355.2 s | 17.6 % | 46.9 % | 53.1 % |
\* « Surface sol » calculée sur l'emprise des points (bounding box du nuage),
à 0.5 m. Les % de trous du MNT final (bornes du header, plus grandes) sont
supérieurs : voir le tile de référence ci-dessus.
## A. Filtres géométriques (stack PDAL actuelle)
- **IGN** (pré-classification fournisseur, classe 2) : le plus rapide (~9 s).
Fiable là où le fournisseur a confiance ; trous sous forêt dense / relief.
Aucun paramètre à régler.
- **SMRF** — Pingel, Clarke & McBride 2013, *ISPRS J. Photogramm. Remote
Sens.* 77:21-30. Filtre **raster** (opère sur un DSM, pas sur les points),
donc plus rapide que les filtres point-based ; **minimise les erreurs de
type I** (omission de sol) → bien adapté quand le sol est rare (forêt).
Meilleure couverture des trois ici (90.4 %) mais ~5.4 min/tile.
- **CSF** — Zhang et al. 2016, *Remote Sensing* 8(6):501. Toile inversée
drapée sur le nuage ; simple, précis, mais **la toile ne touche plus le sol
en terrain raide/vallonné** → mauvaise classification. Le plus lent ici et
le pire sur ce tile. À réserver aux zones urbaines.
- **PTD/PTIN** (Progressive TIN Densification) — Axelsson 2000, ISPRS
Congress. **Gagnant de la littérature** : le plus robuste sur terrain
complexe + forêt (Moudrý et al. 2020, *Measurement* 150:107047 ; Cai et al.
2019, *Remote Sensing* 11(9):1037) et le plus rapide (benchmark lidR :
PTD ~20 s vs CSF ~156 s vs PMF ~1800 s). **NON disponible dans la version
PDAL de cette image** (`filters.ground` / TIN absents) — à ajouter pour
l'utiliser (ou via lidR / une implémentation maison).
## B. Hybride rapide (choisi pour implémentation)
Principe **PTD / Wack & Wimmer** (Wack & Wimmer 2002, *ISPRS Archives*
XXXIV/3A:293-296 : MNT par retour le plus bas, en excluant le 1 % le plus bas
par cellule pour écarter les outliers) :
1. **Base = pré-classification IGN** (classe 2, ~9 s, fiable et officielle).
2. **Comblement mesuré des trous** : pour chaque cellule sans point sol,
prendre le **retour le plus bas robuste** (min du 99 % des points de la
cellule) → ajoute du sol *mesuré* là où le fournisseur a échoué
(rochers, clairières, sol forestier).
3. **Inpainting topographique** des derniers vides (interpolation
terrain-aware déjà implémentée dans `dtm.py:_interpolate_holes`).
Attendu : MNT **continu** (0 % de trous), robuste en forêt/relief,
**~10-15 s/tile** au lieu de 326-355 s. Aucune dépendance GPU, aucun
entraînement.
## C. Modèles IA / ML (supervisés — nécessitent des labels)
Avertissement (Qin et al. 2023, *ISPRS J. Photogramm. Remote Sens.*
202:246-261) : **tout est supervisé** ; le principal risque est la
**généralisation** — un modèle entraîné sur une région dégrade ailleurs.
Aucun filtre DL entièrement non-supervisé publié à date.
**Basés sur les points (3D) :**
| Modèle | Année | Archi | Précision | Vitesse (~/km², GPU) |
|---|---|---|---|---|
| KPConv / RandLA-Net (Qin, OpenGF) | 2021 | KPConv / RandLA-Net | 97.8 % OA, RMSE DTM 0.20 m, IoU sol 95 % | 0.5-2.5 min |
| PFCN (Jin, *IEEE JSTARS* 13:3958) | 2020 | point-FCN | Te 1.73 %, Kappa 93.9 % | ~1/3 du coût PointNet++ |
| Terrain-Net (Li, *Remote Sensing* 14(22):5798) | 2022 | KPConv + self-attention | OA 98 %, mIoU 0.933 | param-free au transfert |
| MSVC (Štroner, *Remote Sensing* 17(4):615) | 2025 | DNN voxel 9x9x9 | bat CSF en F-score | — |
**Rasterisés (sortent directement le MNT — le plus proche du besoin) :**
| Modèle | Année | Archi | Résultat |
|---|---|---|---|
| Precursor (Rizaldy, *ISPRS Annals* IV-2:231) | 2018 | 2D FCN | Te 5.22 %, 78x plus rapide |
| DeepTerRa / ALS2DTM (Lê, *IEEE JSTARS* 15:2778) | 2022 | GAN pix2pix (U-Net) | RMSE MNT < 1 m, filtre + interp en 1 passe |
| DSM2DTM (Bittner, *ISPRS Annals* X-1/W1-2023:925) | 2023 | U-Net (EfficientNet) | masque non-sol + hauteur sol/pixel |
**Jeu de données d'entraînement** : OpenGF (Qin et al., CVPRW 2021,
arXiv:2101.09641 — 47.7 km², 542 M pts) ; ALS2DTM (Lê et al. 2022,
arXiv:2206.03778 — 52 km², 1.66 Md pts, urbain/forêt/montagne).
**Coûts / obstacles pour notre cas** : (1) labels → à générer en
pseudo-labels (sortie SMRF/PTD haute qualité sur un échantillon représentatif
de nos tiles) ou pré-entraînement OpenGF/ALS2DTM ; (2) généralisation sur le
terrain divers de LiDAR HD (plaine/forêt/montagne/urbain) ; (3) infra :
checkpoint + chemin d'inférence GPU dans l'image Docker.
**Meilleur fit si on part sur l'IA** : un **U-Net rasterisé (style
DSM2DTM)** — rasteriser le nuage en grilles multi-canaux (altitude, pente,
courbure, densité, stats de retours), sortir masque sol + hauteur sol.
2D = très rapide et trivial à déployer sur GPU, fusionne filtrage +
interpolation. Le KPConv/RandLA-Net est plus précis en 3D pur mais plus lourd
à déployer.
## Synthèse / décision
- « Rapide » contrainte dure + faible maintenance → **hybride (B)**
(~10-40 s/tile, zéro entraînement, zéro GPU). ← **choix retenu, IMPLÉMENTÉ**
- Base = pré-classification IGN (rapide, ~10 s). `auto` la préfère dès que
≥ 20 % des points sont classés sol (seuil abaissé de 30 % à 20 %, car le
MNT est ensuite complété — voir ci-dessous).
- Le MNT est **toujours** complété dans `create_dtm_fast` : comblement par le
**retour le plus bas par cellule** (`_min_return_grid`, Wack & Wimmer 2002)
pour les trous, puis **interpolation terrain-aware** (`_interpolate_holes`).
Résultat : MNT continu (0 % de trous) pour n'importe quelle base.
- Vérifié sur le tile 0999_6778 : base CSF + comblement → 1,9 M trous
comblés par retour le plus bas, 5,3 % interpolés, MNT 0 % de trous.
- Qualité max dans les cas durs (raide + dense), ~1-2 min/tile + GPU +
entraînement acceptés → **U-Net rasterisé (C)**. (non implémenté)
- Meilleur filtre géométrique disponible dans PDAL → **SMRF (A)** (meilleure
couverture 90.4 % mais 5.4 min/tile). Sélectionnable via `--ground-classification smrf`.
- « Gagnant » absolu de la littérature (rapide + robuste) → **PTD/PTIN
(A)** : à intégrer (pas dans la stack PDAL actuelle).
## Références
- Axelsson (2000), PTIN/PTD, ISPRS Congress.
- Pingel, Clarke & McBride (2013), SMRF, ISPRS J. P&RS 77:21-30.
- Zhang et al. (2016), CSF, Remote Sensing 8(6):501.
- Wack & Wimmer (2002), DMT par retour le plus bas, ISPRS Archives XXXIV/3A.
- Moudrý et al. (2020), comparaison CSF/PTIN/PMF/SMRF, Measurement 150:107047.
- Cai et al. (2019), CS+PTD, Remote Sensing 11(9):1037.
- Qin et al. (2021), OpenGF, CVPR Workshops (arXiv:2101.09641).
- Qin et al. (2023), dataset + evaluation + survey, ISPRS J. P&RS 202:246-261.
- Lê et al. (2022), DeepTerRa/ALS2DTM, IEEE JSTARS 15:2778 (arXiv:2206.03778).
- Bittner et al. (2023), DSM2DTM, ISPRS Annals X-1/W1-2023:925.
- lidR book (comparaison PTD/CSF/PMF) : https://r-lidar.github.io/lidRbook/gnd.html

View File

@ -97,7 +97,10 @@ def main():
) )
parser.add_argument( parser.add_argument(
"input", "input",
help="Dossier contenant les fichiers LAZ/LAS" nargs="?",
default="/data/input",
help="Dossier contenant les fichiers LAZ/LAS (défaut: /data/input ; "
"optionnel pour --rebuild-index)"
) )
parser.add_argument( parser.add_argument(
"-o", "--output", "-o", "--output",
@ -134,7 +137,16 @@ def main():
parser.add_argument( parser.add_argument(
"--force-classification", "--force-classification",
action="store_true", action="store_true",
help="Reclassifier le sol même si le fichier .las existe déjà" help="Reclassifier le sol même si la méthode est inchangée (régénère aussi le DTM "
"et les images). Sans ce flag, changer --ground-classification suffit : la "
"méthode enregistrée est comparée et un changement déclenche la reclassification."
)
parser.add_argument(
"--bare-earth",
action="store_true",
help="Sol nu : ramener le DTM au retour le plus bas de chaque cellule. "
"Requalifie le point le plus bas de chaque colonne en terrain — utile sous "
"végétation dense ou en relief raide où la classification du sol sous-couvre le terrain."
) )
parser.add_argument( parser.add_argument(
"--keep-tif", "--keep-tif",
@ -143,9 +155,22 @@ def main():
) )
parser.add_argument( parser.add_argument(
"--ground-classification", "--ground-classification",
choices=["auto", "smrf", "csf"], choices=["auto", "ign", "smrf", "csf"],
default="auto", default="auto",
help="Méthode de classification du sol : auto (détection), smrf, csf (défaut: auto)" help="Méthode de classification du sol : auto (préfère la pré-classification IGN si "
"présente — base rapide — sinon détection SMRF/CSF), ign, smrf, csf. "
"Avec ign, le MNT est la rasterisation pure des classes choisies "
"(--ign-classes) sans aucune retouche ; avec smrf/csf, il est complété "
"par le retour le plus bas par cellule + interpolation des trous. (défaut: auto)"
)
parser.add_argument(
"--ign-classes",
default="sol",
help="Classes LAS extraites pour le MNT avec la classification IGN (méthode "
"ign/auto) : liste noms ou codes séparés par virgules — "
"sol(2), unclassified(1), eau(9), virtuel(66), pont(17), sursol(64). "
"Ex: --ign-classes sol,unclassified. Changer la liste reclassifie les "
"dalles concernées. (défaut: sol)"
) )
parser.add_argument( parser.add_argument(
"--quality", "--quality",
@ -185,6 +210,14 @@ def main():
default=None, default=None,
help="Traiter un ou plusieurs fichiers LAZ/LAS (nom complet sans extension, ex: LHD_FXX_1000_6882_PTS_LAMB93_IGN69.copc)" help="Traiter un ou plusieurs fichiers LAZ/LAS (nom complet sans extension, ex: LHD_FXX_1000_6882_PTS_LAMB93_IGN69.copc)"
) )
parser.add_argument(
"--fetch-tiles",
nargs="+",
default=None,
metavar="COL,ROW",
help="Télécharger ces dalles LiDAR HD depuis l'IGN avant traitement "
"(tuiles non encore générées, ex: --fetch-tiles 1055,6882 1056,6883)"
)
parser.add_argument( parser.add_argument(
"-v", "--verbose", "-v", "--verbose",
action="store_true", action="store_true",
@ -253,6 +286,21 @@ def main():
logger.warning("Aucune tuile traitée trouvée — carte globale non générée") logger.warning("Aucune tuile traitée trouvée — carte globale non générée")
return return
# Téléchargement des dalles IGN manquantes avant le traitement
if args.fetch_tiles:
from .fetch_ign import fetch_tiles, parse_tile_specs
try:
specs = parse_tile_specs(args.fetch_tiles)
except ValueError as e:
logger.error(str(e))
return
logger.info(f"Téléchargement de {len(specs)} dalle(s) LiDAR HD depuis l'IGN...")
fetched = fetch_tiles(args.input, specs, args.output)
if fetched:
logger.info(f"{len(fetched)} dalle(s) téléchargée(s) — traitement...")
else:
logger.warning("Aucune dalle téléchargée (déjà présentes ou introuvables)")
quality = 100 if args.lossless else args.quality quality = 100 if args.lossless else args.quality
# Parse --only and --skip: accept comma-separated values # Parse --only and --skip: accept comma-separated values
only_viz = None only_viz = None
@ -268,8 +316,10 @@ def main():
workers=args.workers, workers=args.workers,
force=args.force, force=args.force,
ground_method=args.ground_classification, ground_method=args.ground_classification,
ign_classes=args.ign_classes,
force_classify=args.force_classification, force_classify=args.force_classification,
keep_tif=args.keep_tif, keep_tif=args.keep_tif,
bare_earth=args.bare_earth,
quality=quality, quality=quality,
only_viz=only_viz, only_viz=only_viz,
skip_viz=skip_viz, skip_viz=skip_viz,
@ -315,30 +365,8 @@ def main():
logger.info(f"Traitement de {len(unique_files)} fichier(s) sélectionné(s)") logger.info(f"Traitement de {len(unique_files)} fichier(s) sélectionné(s)")
for laz_file in unique_files: for laz_file in unique_files:
logger.info(f" → {laz_file.name}") logger.info(f" → {laz_file.name}")
for laz_file in unique_files: # Réutilise process_all : workers parallèles, résumé, index, nettoyage
pipeline.process_file(laz_file) pipeline.process_all(files=unique_files)
# Clean up temporary files
logger.info("Nettoyage des fichiers temporaires...")
try:
if pipeline.temp_dir.exists():
shutil.rmtree(pipeline.temp_dir)
temp_base = pipeline.output_dir / "temp"
if temp_base.exists():
shutil.rmtree(temp_base)
logger.info(" ✓ Fichiers temporaires supprimés")
except Exception as e:
logger.warning(f" Note: Impossible de supprimer les fichiers temporaires: {e}")
# Génère la carte globale après traitement --file
if not args.no_index:
try:
from .index import build_index
index_path = build_index(pipeline.output_dir, pipeline.output_format)
if index_path:
logger.info(f"Carte globale générée : {index_path}")
except Exception as e:
logger.warning(f"Index global non généré: {e}")
else: else:
pipeline.process_all() pipeline.process_all()
except Exception as e: except Exception as e:

View File

@ -1,7 +1,10 @@
"""DTM generation from classified LiDAR point clouds. """DTM generation from classified LiDAR point clouds.
Handles ground classification via PDAL (SMRF or CSF) and DTM rasterisation Handles ground classification via PDAL (IGN supplier pre-classification,
using scipy binned_statistic_2d. Zones without LiDAR data remain as NaN. SMRF or CSF) and DTM rasterisation
using scipy binned_statistic_2d. Gaps without LiDAR data (common in
complex/rocky terrain) are filled with a terrain-aware interpolation so the
DTM stays continuous.
""" """
import json import json
@ -16,8 +19,67 @@ from scipy.stats import binned_statistic_2d
logger = logging.getLogger("lidar") logger = logging.getLogger("lidar")
# Classes LAS exploitables de la pré-classification LiDAR HD (noms → codes)
IGN_CLASS_NAMES = {
"sol": 2,
"unclassified": 1,
"non-classe": 1,
"eau": 9,
"virtuel": 66,
"pont": 17,
"sursol": 64,
}
def _create_ground_pipeline(input_laz, output_las, method):
def parse_ign_classes(spec):
"""Convertit une liste de classes IGN (noms ou codes) en codes LAS triés.
Args:
spec: Chaîne séparée par virgules, ex. "sol,unclassified" ou "2,1".
Returns:
Liste triée de codes LAS uniques.
Raises:
ValueError: Si un élément n'est ni un nom connu ni un code LAS 0-255,
ou si la liste est vide.
"""
codes = set()
for token in str(spec).split(","):
token = token.strip().lower()
if not token:
continue
if token in IGN_CLASS_NAMES:
codes.add(IGN_CLASS_NAMES[token])
else:
try:
code = int(token)
except ValueError:
raise ValueError(
f"Classe IGN inconnue: '{token}' "
f"(noms: {', '.join(sorted(IGN_CLASS_NAMES))} ou code LAS 0-255)")
if not 0 <= code <= 255:
raise ValueError(f"Code LAS hors bornes (0-255): {code}")
codes.add(code)
if not codes:
raise ValueError("Aucune classe IGN fournie")
return sorted(codes)
def ign_method_label(codes):
"""Étiquette de méthode encodant les classes IGN (ex. 'ign_1_2').
'ign' seul = sol uniquement (code 2), rétrocompatible avec les fichiers de
classification existants. Toute autre combinaison est encodée dans le nom
pour invalider le cache et déclencher la reclassification.
"""
codes = sorted(codes)
if codes == [2]:
return "ign"
return "ign_" + "_".join(str(c) for c in codes)
def _create_ground_pipeline(input_laz, output_las, method, ign_codes=None):
"""Create a PDAL pipeline JSON for ground classification. """Create a PDAL pipeline JSON for ground classification.
All methods include a ReturnNumber/NumberOfReturns >= 1 filter to handle All methods include a ReturnNumber/NumberOfReturns >= 1 filter to handle
@ -33,7 +95,10 @@ def _create_ground_pipeline(input_laz, output_las, method):
Args: Args:
input_laz: Path to input LAZ/LAS file. input_laz: Path to input LAZ/LAS file.
output_las: Path to output classified LAS file. output_las: Path to output classified LAS file.
method: Ground classification method ('smrf' or 'csf'). method: Ground classification method ('ign', 'smrf' or 'csf').
ign_codes: LAS class codes to extract with the 'ign' method
(default: [2] = sol). Multiple ranges on Classification are
logically ORed by filters.range (documented PDAL semantics).
Returns: Returns:
JSON string of the PDAL pipeline. JSON string of the PDAL pipeline.
@ -44,6 +109,38 @@ def _create_ground_pipeline(input_laz, output_las, method):
"limits": "ReturnNumber[1:],NumberOfReturns[1:]" "limits": "ReturnNumber[1:],NumberOfReturns[1:]"
} }
# Classification filter (ground points only)
ground_filter = {
"type": "filters.range",
"limits": "Classification[2:2]"
}
# LiDAR HD IGN : le fichier est pré-classifié par le fournisseur.
# On réutilise la classification telle quelle (mode pur) : les classes
# extraites sont paramétrables — par défaut le sol seul (2), mais on peut
# ajouter p.ex. unclassified (1) pour combler les trous sans retouche.
# Les plages multiples sur Classification sont combinées en OU logique
# par filters.range (sémantique PDAL documentée).
if method == 'ign':
codes = sorted(ign_codes) if ign_codes else [2]
ign_filter = {
"type": "filters.range",
"limits": ",".join(f"Classification[{c}:{c}]" for c in codes)
}
pipeline = {
"pipeline": [
str(input_laz),
return_filter,
ign_filter,
{
"type": "writers.las",
"filename": str(output_las),
"extra_dims": "all"
}
]
}
return json.dumps(pipeline)
# Reset Classification to 0 before preprocessing # Reset Classification to 0 before preprocessing
reset_classification = { reset_classification = {
"type": "filters.assign", "type": "filters.assign",
@ -68,12 +165,6 @@ def _create_ground_pipeline(input_laz, output_las, method):
"multiplier": 3.0 "multiplier": 3.0
} }
# Classification filter (ground points only)
ground_filter = {
"type": "filters.range",
"limits": "Classification[2:2]"
}
# Method-specific ground classification filter # Method-specific ground classification filter
if method == 'smrf': if method == 'smrf':
ground_step = { ground_step = {
@ -85,9 +176,12 @@ def _create_ground_pipeline(input_laz, output_las, method):
"scalar": 1.25 "scalar": 1.25
} }
elif method == 'csf': elif method == 'csf':
# resolution 1.0 m : un cloth à 0.5 m (= 4 M particules pour 1 km²)
# rend la classification ~4× plus lente sans gain visible sur le MNT
# (la résolution finale du MNT vient de la rasterisation, pas du cloth).
ground_step = { ground_step = {
"type": "filters.csf", "type": "filters.csf",
"resolution": 0.5, "resolution": 1.0,
"rigidness": 3, "rigidness": 3,
"smooth": True, "smooth": True,
"threshold": 0.5 "threshold": 0.5
@ -119,6 +213,11 @@ def create_smrf_pipeline(input_laz, output_las):
return _create_ground_pipeline(input_laz, output_las, 'smrf') return _create_ground_pipeline(input_laz, output_las, 'smrf')
def create_ign_pipeline(input_laz, output_las):
"""Create a PDAL pipeline JSON using the IGN supplier pre-classification."""
return _create_ground_pipeline(input_laz, output_las, 'ign')
def create_csf_pipeline(input_laz, output_las): def create_csf_pipeline(input_laz, output_las):
"""Create a PDAL pipeline JSON for CSF ground classification.""" """Create a PDAL pipeline JSON for CSF ground classification."""
return _create_ground_pipeline(input_laz, output_las, 'csf') return _create_ground_pipeline(input_laz, output_las, 'csf')
@ -258,7 +357,7 @@ def detect_ground_method(laz_file):
laz_file: Path to input LAZ/LAS file. laz_file: Path to input LAZ/LAS file.
Returns: Returns:
String: 'smrf' or 'csf' String: 'ign', 'smrf' or 'csf'
""" """
import laspy import laspy
@ -280,6 +379,22 @@ def detect_ground_method(laz_file):
logger.warning(f" Nuage vide (0 points) — méthode par défaut: SMRF") logger.warning(f" Nuage vide (0 points) — méthode par défaut: SMRF")
return 'smrf' return 'smrf'
# LiDAR HD IGN : les données livrées sont pré-classifiées par le fournisseur
# (classe 2 = sol). C'est la base la plus rapide (~10 s) et de référence.
# Le MNT est ensuite complété par le retour le plus bas par cellule +
# interpolation (voir create_dtm_fast), ce qui « rattrape » les trous de la
# pré-classification (forêt dense / relief). On la préfère donc dès qu'une
# part raisonnable des points est classée sol, plutôt que de refiltrer.
try:
cls = np.asarray(las.classification, dtype=np.int32)
ground_ratio = float(np.mean(cls == 2))
except Exception:
ground_ratio = 0.0
if ground_ratio >= 0.2:
logger.info(f" → Méthode: IGN (pré-classification fournisseur — "
f"{ground_ratio * 100:.1f}% de points classe 2)")
return 'ign'
z = np.array(las.z) z = np.array(las.z)
# Height variance (always available) # Height variance (always available)
@ -317,14 +432,17 @@ def detect_ground_method(laz_file):
return method return method
def classify_ground(laz_file, temp_dir, method='auto', force=False): def classify_ground(laz_file, temp_dir, method='auto', force=False, ign_classes="sol"):
"""Classify ground points using PDAL ground classification filter. """Classify ground points using PDAL ground classification filter.
Args: Args:
laz_file: Path to input LAZ/LAS file. laz_file: Path to input LAZ/LAS file.
temp_dir: Directory for temporary files (pipeline.json, ground.las). temp_dir: Directory for temporary files (pipeline.json, ground.las).
method: Ground classification method ('auto', 'smrf', or 'csf'). method: Ground classification method ('auto', 'ign', 'smrf' or 'csf').
force: If True, reclassify even if output file already exists. force: If True, reclassify even if output file already exists.
ign_classes: Classes LAS extraites par la méthode IGN (noms ou codes
séparés par virgules, ex. "sol,unclassified"). Ignoré pour les
autres méthodes.
Returns: Returns:
Path to classified ground LAS file, or None on failure. Path to classified ground LAS file, or None on failure.
@ -338,11 +456,17 @@ def classify_ground(laz_file, temp_dir, method='auto', force=False):
else: else:
logger.info(f" Classification sol: {method.upper()} (forcé)") logger.info(f" Classification sol: {method.upper()} (forcé)")
# Les classes IGN sont encodées dans le nom de fichier (ex. ign_1_2)
# pour qu'un changement de classes invalide le cache et déclenche la
# reclassification.
ign_codes = parse_ign_classes(ign_classes) if method == 'ign' else None
method_label = ign_method_label(ign_codes) if ign_codes else method
# Use shared basename extraction function # Use shared basename extraction function
from .pipeline import _file_basename from .pipeline import _file_basename
laz_base = _file_basename(laz_file) laz_base = _file_basename(laz_file)
output_las = temp_dir / f"{laz_base}_ground_{method}.las" output_las = temp_dir / f"{laz_base}_ground_{method_label}.las"
if output_las.exists() and not force: if output_las.exists() and not force:
logger.info(f" Classification {method.upper()} déjà effectuée — fichier existant réutilisé") logger.info(f" Classification {method.upper()} déjà effectuée — fichier existant réutilisé")
@ -352,8 +476,8 @@ def classify_ground(laz_file, temp_dir, method='auto', force=False):
logger.info(f" Reclassification forcée — suppression de {output_las.name}") logger.info(f" Reclassification forcée — suppression de {output_las.name}")
output_las.unlink() output_las.unlink()
pipeline_json = _create_ground_pipeline(laz_file, output_las, method) pipeline_json = _create_ground_pipeline(laz_file, output_las, method, ign_codes=ign_codes)
pipeline_file = temp_dir / f"pipeline_{method}.json" pipeline_file = temp_dir / f"pipeline_{method_label}.json"
with open(pipeline_file, 'w') as f: with open(pipeline_file, 'w') as f:
f.write(pipeline_json) f.write(pipeline_json)
@ -367,9 +491,9 @@ def classify_ground(laz_file, temp_dir, method='auto', force=False):
if output_las.exists() and output_las.stat().st_size < 100: if output_las.exists() and output_las.stat().st_size < 100:
logger.error(f" ✗ Fichier ground vide (taille < 100 octets)") logger.error(f" ✗ Fichier ground vide (taille < 100 octets)")
output_las.unlink(missing_ok=True) output_las.unlink(missing_ok=True)
# Fallback: if CSF produced no ground points, retry with SMRF # Fallback: si la méthode ne produit aucun point sol, réessayer avec SMRF
if method == 'csf': if method in ('csf', 'ign'):
return _fallback_to_smrf(laz_file, temp_dir, laz_base, force) return _fallback_to_smrf(laz_file, temp_dir, laz_base, force, source=method_label)
return None return None
logger.info(f" ✓ Classification sol {method.upper()} terminée") logger.info(f" ✓ Classification sol {method.upper()} terminée")
return output_las return output_las
@ -377,9 +501,9 @@ def classify_ground(laz_file, temp_dir, method='auto', force=False):
error_msg = e.stderr.decode() if e.stderr else str(e) error_msg = e.stderr.decode() if e.stderr else str(e)
logger.warning(f" ✗ Erreur classification PDAL ({method.upper()}): {error_msg}") logger.warning(f" ✗ Erreur classification PDAL ({method.upper()}): {error_msg}")
# Fallback: if CSF failed, retry with SMRF # Fallback: si CSF ou la pré-classification échouent, réessayer avec SMRF
if method == 'csf': if method in ('csf', 'ign'):
return _fallback_to_smrf(laz_file, temp_dir, laz_base, force) return _fallback_to_smrf(laz_file, temp_dir, laz_base, force, source=method_label)
# Try repairing file with laspy if PDAL fails on EVLR/VLR # Try repairing file with laspy if PDAL fails on EVLR/VLR
if 'VLR' in error_msg or 'Invalid' in error_msg: if 'VLR' in error_msg or 'Invalid' in error_msg:
@ -405,28 +529,30 @@ def classify_ground(laz_file, temp_dir, method='auto', force=False):
return None return None
def _fallback_to_smrf(laz_file, temp_dir, laz_base, force=False): def _fallback_to_smrf(laz_file, temp_dir, laz_base, force=False, source='csf'):
"""Retry ground classification with SMRF when CSF fails. """Retry ground classification with SMRF when CSF/IGN fails.
CSF (Cloth Simulation Filter) can fail on certain terrain types where CSF (Cloth Simulation Filter) can fail on certain terrain types where
SMRF (Simple Morphological Filter) succeeds. This fallback ensures SMRF (Simple Morphological Filter) succeeds, and a file without usable
processing continues even when auto-detection selects CSF incorrectly. pre-classification produces an empty ground extract. This fallback ensures
processing continues even when the selected method fails.
Args: Args:
laz_file: Path to input LAZ/LAS file. laz_file: Path to input LAZ/LAS file.
temp_dir: Directory for temporary files. temp_dir: Directory for temporary files.
laz_base: Base name for the file. laz_base: Base name for the file.
force: If True, reclassify even if output exists. force: If True, reclassify even if output exists.
source: Method that failed ('csf' or 'ign').
Returns: Returns:
Path to classified ground LAS file, or None on failure. Path to classified ground LAS file, or None on failure.
""" """
logger.info(f" → Basculement CSF → SMRF (fallback)") logger.info(f" → Basculement {source.upper()} → SMRF (fallback)")
# Clean up failed CSF output if it exists # Clean up failed output if it exists
csf_output = temp_dir / f"{laz_base}_ground_csf.las" failed_output = temp_dir / f"{laz_base}_ground_{source}.las"
if csf_output.exists(): if failed_output.exists():
csf_output.unlink(missing_ok=True) failed_output.unlink(missing_ok=True)
output_las = temp_dir / f"{laz_base}_ground_smrf.las" output_las = temp_dir / f"{laz_base}_ground_smrf.las"
@ -481,7 +607,118 @@ def _repair_laz_with_laspy(input_laz, output_las):
return False return False
def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, output_suffix=""): def _interpolate_holes(dtm, downsample=8):
"""Fill remaining NaN holes with a terrain-aware surface interpolation.
In complex / rocky terrain the ground under-classification leaves interior
holes far too large for a 1 m gap fill, which otherwise become flat
nearest-neighbor patches in the downstream layers. This helper triangulates
the valid cells on a downsampled grid (linear, nearest as a fallback for
cells outside the data hull) and bilinearly upsamples the result, keeping
the operation fast even for large rasters.
Args:
dtm: 2-D float array (may contain NaN holes).
downsample: Coarsening factor for the interpolation grid.
Returns:
Tuple (filled_array, filled_count).
"""
holes = np.isnan(dtm)
if not holes.any():
return dtm, 0
valid = ~holes
if not valid.any():
return dtm, 0
height, width = dtm.shape
step = max(1, downsample)
coarse = dtm[::step, ::step].astype(np.float64)
c_valid = ~np.isnan(coarse)
c_holes = np.isnan(coarse)
if not c_holes.any() or not c_valid.any():
return dtm, 0
from scipy.interpolate import griddata
from scipy.ndimage import map_coordinates
cy, cx = np.where(c_valid)
c_coords = np.column_stack([cx, cy]).astype(np.float64)
c_vals = coarse[c_valid]
hy, hx = np.where(c_holes)
h_coords = np.column_stack([hx, hy]).astype(np.float64)
interp = griddata(c_coords, c_vals, h_coords, method='linear')
bad = np.isnan(interp)
if bad.any():
interp[bad] = griddata(c_coords, c_vals, h_coords[bad], method='nearest')
coarse_filled = coarse.copy()
coarse_filled[c_holes] = interp
# Coarse cell i represents fine column/row i*step, so fine index c maps to
# coarse coordinate c/step (no half-cell offset).
rows = np.arange(height) / step
cols = np.arange(width) / step
grid_y, grid_x = np.meshgrid(rows, cols, indexing='ij')
upsampled = map_coordinates(coarse_filled, [grid_y, grid_x], order=1)
filled = dtm.copy()
filled[holes] = upsampled[holes]
return filled, int(holes.sum())
def _min_return_grid(laz_file, width, height, bounds, chunk_size=2_000_000):
"""Rasterize the per-cell minimum z (lowest return) of the full point cloud.
In complex/forested terrain the ground is under-classified, leaving DTM
holes. Filling them with the *lowest measured return* of the cell (Wack &
Wimmer 2002) recovers a real ground surface (forest floor, rock, clearing)
instead of a pure interpolation. The read is streamed in chunks so memory
stays bounded to the output grid regardless of the point count.
Args:
laz_file: Path to the full (unclassified) LAZ/LAS file.
width, height: Output grid dimensions (pixels).
bounds: (min_x, min_y, max_x, max_y) the grid covers.
chunk_size: Points per streaming chunk.
Returns:
(height, width) float32 array of per-cell min z (NaN where no point).
"""
import laspy
min_x, min_y, max_x, max_y = bounds
grid = np.full((height, width), np.nan, dtype=np.float32)
rng = [[min_x, max_x], [min_y, max_y]]
def process(points):
if len(points) == 0:
return
x = np.asarray(points.x, dtype=np.float64)
y = np.asarray(points.y, dtype=np.float64)
z = np.asarray(points.z, dtype=np.float64)
st = binned_statistic_2d(x, y, z, statistic='min',
bins=[width, height], range=rng)
# Match the DTM convention: .T then flip Y (north at top).
cell_min = st.statistic.T[::-1, :].astype(np.float32)
# fmin ignores NaN so cells without a point in this chunk stay NaN.
np.fmin(grid, cell_min, out=grid)
try:
with laspy.open(str(laz_file)) as las:
for chunk in las.chunk_iterator(chunk_size):
process(chunk)
except Exception as e:
logger.warning(f" Lecture streaming impossible ({e}) — lecture complète")
las = _read_with_pdal(laz_file)
if las is None:
return grid
process(las)
return grid
def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False,
output_suffix="", source_laz=None, bare_earth=False,
pure=False):
"""Create DTM using fast binning method with gap filling. """Create DTM using fast binning method with gap filling.
Args: Args:
@ -491,6 +728,15 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, output
resolution: Grid resolution in meters per pixel. resolution: Grid resolution in meters per pixel.
force: If True, regenerate even if DTM already exists. force: If True, regenerate even if DTM already exists.
output_suffix: Suffix for output filename (e.g. '_r0p2' for additional resolutions). output_suffix: Suffix for output filename (e.g. '_r0p2' for additional resolutions).
source_laz: Optionnel : chemin du LAZ complet (non classé). Utilisé
uniquement avec bare_earth (plancher au retour le plus bas).
bare_earth: If True, pull the DTM down to the lowest measured return of
each cell (bare-earth floor). This requalifies the lowest point of
every column as terrain, recovering the ground under dense
vegetation / steep relief that the ground classifier rejected.
pure: Sans effet (conservé pour compatibilité). Fonctionnement
historique rétabli : petits trous comblés par fillnodata, grands
trous laissés en nodata (rendus en noir dans les rendus).
Returns: Returns:
Path to output DTM GeoTIFF, or None on failure. Path to output DTM GeoTIFF, or None on failure.
@ -542,7 +788,19 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, output
dtm = stat.statistic.T dtm = stat.statistic.T
dtm = dtm[::-1, :] # Flip Y so north is at top dtm = dtm[::-1, :] # Flip Y so north is at top
# Fill small gaps (< 1m from existing data) while keeping large gaps as NaN # Comblement « historique » (fonctionnement d'origine, rétabli) :
# seuls les petits trous proches des données sont remplis ; les grands
# trous restent en nodata et apparaissent en noir dans les rendus.
# Le plancher au retour le plus bas n'est appliqué qu'à la demande
# explicite (--bare-earth).
if bare_earth and source_laz is not None:
min_grid = _min_return_grid(source_laz, width, height,
(min_x, min_y, max_x, max_y))
lower = ~np.isnan(min_grid) & (min_grid < dtm)
dtm = np.where(lower, min_grid, dtm)
logger.info(f" Sol nu : {int(lower.sum()):,} cellules raménées au retour le plus bas")
# Fill small gaps (< 1 m from data) precisely — comme avant
nan_count = np.count_nonzero(np.isnan(dtm)) nan_count = np.count_nonzero(np.isnan(dtm))
if nan_count > 0: if nan_count > 0:
total = dtm.size total = dtm.size
@ -558,8 +816,6 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, output
if filled_count > 0: if filled_count > 0:
dtm = np.where(small_gap_mask, dtm_filled, dtm) dtm = np.where(small_gap_mask, dtm_filled, dtm)
logger.info(f" {filled_count:,} petits trous comblés (< {max_gap_pixels}px)") logger.info(f" {filled_count:,} petits trous comblés (< {max_gap_pixels}px)")
remaining = np.count_nonzero(np.isnan(dtm))
logger.info(f" {remaining:,} pixels restent sans données (grands écarts)")
# Save as GeoTIFF # Save as GeoTIFF
output_tif = dtm_dir / f"{basename}_dtm{output_suffix}.tif" output_tif = dtm_dir / f"{basename}_dtm{output_suffix}.tif"

162
lidar_pipeline/fetch_ign.py Normal file
View File

@ -0,0 +1,162 @@
"""Téléchargement des dalles LiDAR HD de l'IGN pour les tuiles non générées.
Catalogue STAC (à jour) : https://browser.stac.teledetection.fr/collections/lidarhd
API : https://api.stac.teledetection.fr/collections/lidarhd/items
Fichiers (géoplateforme): https://data.geopf.fr/telechargement/download/...
Chaque dalle couvre 1 km × 1 km en Lambert 93 et est nommée par son coin
nord-ouest : LHD_FXX_{col}_{row}_PTS_LAMB93_IGN69.copc.laz
(col = X ouest en km, row = Y nord en km, cf. propriété STAC
"lidarhd:coordonnees_NW" au format "0816-6847").
"""
import json
import logging
import time
import urllib.parse
import urllib.request
from pathlib import Path
logger = logging.getLogger("lidar")
_STAC_ITEMS_URL = "https://api.stac.teledetection.fr/collections/lidarhd/items"
_HEADERS = {"User-Agent": "Mozilla/5.0 (lidar-archeo-pipeline)"}
def parse_tile_specs(args):
"""Convertit des spécifications "col,row" ou "col:row" en liste de tuples.
Args:
args: liste de chaînes (ex: ["1055,6882", "1056:6883"]).
Returns:
Liste de tuples (col, row).
Raises:
ValueError: si une spécification est mal formée.
"""
specs = []
for raw in args:
text = raw.strip().replace(":", ",").replace(";", ",")
parts = [p.strip() for p in text.split(",") if p.strip()]
if len(parts) != 2:
raise ValueError(f"Spécification de tuile invalide: {raw!r} (attendu: col,row)")
try:
col, row = int(parts[0]), int(parts[1])
except ValueError:
raise ValueError(f"Spécification de tuile invalide: {raw!r} (col et row doivent être des entiers)")
specs.append((col, row))
return specs
def tile_filename(col, row):
"""Nom de fichier LAZ standard d'une dalle (col, row)."""
return f"LHD_FXX_{col:04d}_{row:04d}_PTS_LAMB93_IGN69.copc.laz"
def _bbox_wgs84(col, row):
"""Bbox WGS84 de la dalle (col,row) pour la requête STAC (peut être élargie)."""
try:
from rasterio.warp import transform as warp_transform
xs = [col * 1000, (col + 1) * 1000, col * 1000, (col + 1) * 1000]
ys = [(row - 1) * 1000] * 2 + [row * 1000] * 2
lons, lats = warp_transform('EPSG:2154', 'EPSG:4326', xs, ys)
except Exception:
from .index import _approx_l93_to_wgs84
pts = [_approx_l93_to_wgs84(x, y)
for x in (col * 1000, (col + 1) * 1000)
for y in ((row - 1) * 1000, row * 1000)]
lons = [p[0] for p in pts]
lats = [p[1] for p in pts]
pad = 0.005 # ~500 m de marge pour éviter les erreurs d'arrondi aux bords
return (min(lons) - pad, min(lats) - pad, max(lons) + pad, max(lats) + pad)
def match_feature(features, col, row):
"""Retourne l'item STAC correspondant à la dalle (col,row), sinon None."""
want = f"{col:04d}-{row:04d}"
for feature in features:
props = feature.get("properties", {})
if props.get("lidarhd:coordonnees_NW") == want:
return feature
return None
def find_tile_url(col, row, timeout=20):
"""Cherche l'URL de téléchargement de la dalle (col,row) dans le catalogue STAC.
Returns:
URL (str) ou None si la dalle n'est pas (encore) publiée par l'IGN.
"""
w, s, e, n = _bbox_wgs84(col, row)
query = urllib.parse.urlencode({"bbox": f"{w:.6f},{s:.6f},{e:.6f},{n:.6f}", "limit": 50})
req = urllib.request.Request(f"{_STAC_ITEMS_URL}?{query}", headers=_HEADERS)
with urllib.request.urlopen(req, timeout=timeout) as response:
data = json.loads(response.read().decode("utf-8"))
feature = match_feature(data.get("features", []), col, row)
if not feature:
return None
return feature.get("assets", {}).get("data", {}).get("href")
def download_file(url, dest_path, timeout=120, chunk=1024 * 1024):
"""Télécharge url vers dest_path en streaming. Retourne la taille en octets."""
req = urllib.request.Request(url, headers=_HEADERS)
t0 = time.time()
with urllib.request.urlopen(req, timeout=timeout) as response, open(dest_path, "wb") as out:
done = 0
while True:
block = response.read(chunk)
if not block:
break
out.write(block)
done += len(block)
elapsed = time.time() - t0
logger.info(f" {done / 1e6:.0f} Mo en {elapsed:.0f}s"
f" ({done / 1e6 / max(elapsed, 0.1):.1f} Mo/s)")
return done
def fetch_tiles(input_dir, specs, output_dir=None):
"""Télécharge les dalles IGN spécifiées, sauf celles déjà présentes/générées.
Args:
input_dir: dossier des fichiers LAZ (écriture autorisée requise).
specs: liste de tuples (col, row).
output_dir: dossier de sortie (optionnel) — permet d'ignorer les
tuiles dont les visualisations existent déjà.
Returns:
Liste des chemins téléchargés.
"""
input_dir = Path(input_dir)
downloaded = []
for col, row in specs:
name = tile_filename(col, row)
dest = input_dir / name
if dest.exists():
logger.info(f" {name} : déjà présent dans input/ — aucun téléchargement")
continue
if output_dir is not None:
vis_dir = Path(output_dir) / "visualisations"
if list(vis_dir.glob(f"LHD_FXX_{col:04d}_{row:04d}_PTS*")):
logger.info(f" {name} : visualisations déjà générées — ignorée")
continue
logger.info(f" {name} : recherche dans le catalogue IGN...")
try:
url = find_tile_url(col, row)
except Exception as e:
logger.warning(f" ✗ {name} : erreur catalogue ({e})")
continue
if not url:
logger.warning(f" ✗ {name} : introuvable dans le catalogue IGN (zone non publiée ?)")
continue
logger.info(f" {name} : téléchargement depuis la géoplateforme...")
try:
download_file(url, dest)
logger.info(f" ✓ {name} téléchargée")
downloaded.append(dest)
except Exception as e:
dest.unlink(missing_ok=True)
logger.warning(f" ✗ {name} : échec du téléchargement ({e})")
return downloaded

File diff suppressed because it is too large Load Diff

View File

@ -109,7 +109,7 @@ VIZ_STEPS = [
class LidarArchaeoPipeline: class LidarArchaeoPipeline:
"""Orchestrates the LiDAR archaeological analysis pipeline.""" """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, quality=98, only_viz=None, skip_viz=None, output_format='avif', gpu_ids=None, no_index=False): def __init__(self, input_dir, output_dir, resolution=0.5, workers=1, force=False, ground_method='auto', ign_classes="sol", force_classify=False, keep_tif=False, bare_earth=False, quality=98, only_viz=None, skip_viz=None, output_format='avif', gpu_ids=None, no_index=False):
self.input_dir = Path(input_dir) self.input_dir = Path(input_dir)
self.output_dir = Path(output_dir) self.output_dir = Path(output_dir)
# Accept single float or comma-separated string for multi-resolution # Accept single float or comma-separated string for multi-resolution
@ -123,7 +123,9 @@ class LidarArchaeoPipeline:
self.workers = workers self.workers = workers
self.force = force self.force = force
self.ground_method = ground_method self.ground_method = ground_method
self.ign_classes = ign_classes
self.force_classify = force_classify self.force_classify = force_classify
self.bare_earth = bare_earth
self.keep_tif = keep_tif self.keep_tif = keep_tif
self.quality = quality self.quality = quality
self.only_viz = only_viz self.only_viz = only_viz
@ -212,7 +214,7 @@ class LidarArchaeoPipeline:
else: else:
return file_vis_dir / f"{basename}_{name}.{ext}" return file_vis_dir / f"{basename}_{name}.{ext}"
def generate_all_visualizations(self, dtm_file, basename, resolution=None, vis_dir=None): def generate_all_visualizations(self, dtm_file, basename, resolution=None, vis_dir=None, force_images=None):
"""Generate all archaeological visualizations for one DTM file. """Generate all archaeological visualizations for one DTM file.
Optimisation: SharedDEM is only computed if at least one visualization Optimisation: SharedDEM is only computed if at least one visualization
@ -229,9 +231,10 @@ class LidarArchaeoPipeline:
total = len(self.viz_steps) total = len(self.viz_steps)
# Phase 1: determine which visualizations need generation # Phase 1: determine which visualizations need generation
force_viz = self.force if force_images is None else force_images
needs_generation = {} # name -> True/False needs_generation = {} # name -> True/False
for name, func in self.viz_steps: for name, func in self.viz_steps:
if self.force: if force_viz:
needs_generation[name] = True needs_generation[name] = True
else: else:
expected_webp = self._expected_output_path(name, basename, file_vis_dir, self.output_format) expected_webp = self._expected_output_path(name, basename, file_vis_dir, self.output_format)
@ -265,8 +268,8 @@ class LidarArchaeoPipeline:
vis_results[name] = self._expected_output_path(name, basename, file_vis_dir, self.output_format) vis_results[name] = self._expected_output_path(name, basename, file_vis_dir, self.output_format)
continue continue
# When --force, delete existing TIF to ensure clean regeneration # When regenerating, delete existing TIF to ensure clean regeneration
if self.force: if force_viz:
for tif in file_vis_dir.glob(f"{basename}_{name}.tif"): for tif in file_vis_dir.glob(f"{basename}_{name}.tif"):
tif.unlink(missing_ok=True) tif.unlink(missing_ok=True)
if name == 'pos_open': if name == 'pos_open':
@ -325,6 +328,48 @@ class LidarArchaeoPipeline:
res_str = f"{resolution}".replace('.', 'p') res_str = f"{resolution}".replace('.', 'p')
return f"_r{res_str}" return f"_r{res_str}"
def _dtm_method_path(self, basename, res_suffix):
"""Sidecar path storing which ground method produced a DTM."""
return self.dtm_dir / f"{basename}_dtm{res_suffix}_method.txt"
def _dtm_method_name(self, basename, res_suffix):
"""Read the recorded ground method for a DTM, or None if unknown."""
p = self._dtm_method_path(basename, res_suffix)
if p.exists():
try:
return p.read_text(encoding="utf-8").strip() or None
except Exception:
return None
return None
def _effective_ground_method(self):
"""Méthode effective pour le suivi de cache, classes IGN incluses.
La méthode 'ign' est étiquetée avec les classes choisies (ex. 'ign_1_2')
pour qu'un changement de --ign-classes déclenche la reclassification.
"""
if self.ground_method == 'ign':
from .dtm import parse_ign_classes, ign_method_label
return ign_method_label(parse_ign_classes(self.ign_classes))
return self.ground_method
def _dtm_method_matches(self, basename, res_suffix):
"""True if the recorded ground method matches the requested one.
A DTM without a recorded method is treated as matching so the existing
cache is preserved; the method is adopted on its next reclassification.
"""
recorded = self._dtm_method_name(basename, res_suffix)
return recorded is None or recorded == self._effective_ground_method()
def _write_dtm_method(self, basename, res_suffix):
"""Record the ground classification method used to build a DTM."""
try:
self._dtm_method_path(basename, res_suffix).write_text(
self._effective_ground_method(), encoding="utf-8")
except Exception:
pass
def process_file(self, laz_file): def process_file(self, laz_file):
"""Process a single LAZ file through the full pipeline. """Process a single LAZ file through the full pipeline.
@ -349,10 +394,18 @@ class LidarArchaeoPipeline:
# Step 1: Ground classification (shared across all resolutions) # Step 1: Ground classification (shared across all resolutions)
las_file = None las_file = None
t_classif = 0 t_classif = 0
dtm_rebuilt = False
# The ground method is shared across resolutions. It is recorded per DTM
# in a sidecar so that changing --ground-classification invalidates the
# cache; otherwise the cached DTM would be reused and the new method
# never applied (nothing would change visually).
primary_suffix = self._res_suffix(self.resolutions[0])
method_matches = self._dtm_method_matches(basename, primary_suffix)
for i, res in enumerate(self.resolutions): for i, res in enumerate(self.resolutions):
res_suffix = self._res_suffix(res) res_suffix = self._res_suffix(res)
dtm_path = self.dtm_dir / f"{basename}_dtm{res_suffix}.tif" dtm_path = self.dtm_dir / f"{basename}_dtm{res_suffix}.tif"
if dtm_path.exists(): if dtm_path.exists() and not self.force_classify:
if method_matches:
import rasterio import rasterio
try: try:
with rasterio.open(dtm_path) as src: with rasterio.open(dtm_path) as src:
@ -370,13 +423,16 @@ class LidarArchaeoPipeline:
except Exception: except Exception:
logger.warning(f"Impossible de lire le DTM existant — régénération") logger.warning(f"Impossible de lire le DTM existant — régénération")
dtm_path.unlink() dtm_path.unlink()
else:
logger.info(f" DTM{res_suffix} produit par {self._dtm_method_name(basename, primary_suffix) or '?'} ≠ {self._effective_ground_method()} → reclassification")
dtm_path.unlink()
# Need to classify/generate DTM for this resolution # Need to classify/generate DTM for this resolution
if las_file is None: if las_file is None:
# First time: do ground classification # First time: do ground classification
logger.info("[1/5] Classification du sol...") logger.info("[1/5] Classification du sol...")
t1 = time.time() t1 = time.time()
las_file = classify_ground(laz_file, self.temp_dir, method=self.ground_method, force=self.force_classify) las_file = classify_ground(laz_file, self.temp_dir, method=self.ground_method, force=self.force_classify, ign_classes=self.ign_classes)
t_classif = time.time() - t1 t_classif = time.time() - t1
if not las_file: if not las_file:
logger.error(f" ✗ Échec classification ({t_classif:.1f}s)") logger.error(f" ✗ Échec classification ({t_classif:.1f}s)")
@ -386,7 +442,15 @@ class LidarArchaeoPipeline:
# Generate DTM at this resolution # Generate DTM at this resolution
logger.info(f"{'[2/5]' if i == 0 else ' '} Génération DTM {res}m/px...") logger.info(f"{'[2/5]' if i == 0 else ' '} Génération DTM {res}m/px...")
t2 = time.time() t2 = time.time()
dtm_file = create_dtm_fast(las_file, basename, self.dtm_dir, res, force=self.force, output_suffix=res_suffix) # Classification IGN → mode pur : DTM = rasterisation brute des
# classes choisies, sans plancher ni comblement (sauf --bare-earth)
pure_ign = "_ground_ign" in Path(las_file).name
dtm_file = create_dtm_fast(las_file, basename, self.dtm_dir, res,
force=self.force or self.force_classify,
output_suffix=res_suffix,
source_laz=laz_file,
bare_earth=self.bare_earth,
pure=pure_ign)
t_dtm = time.time() - t2 t_dtm = time.time() - t2
if not dtm_file: if not dtm_file:
logger.error(f" ✗ Échec DTM {res}m/px ({t_dtm:.1f}s)") logger.error(f" ✗ Échec DTM {res}m/px ({t_dtm:.1f}s)")
@ -394,6 +458,8 @@ class LidarArchaeoPipeline:
return False # Primary resolution failure is fatal return False # Primary resolution failure is fatal
continue # Additional resolution failure is non-fatal continue # Additional resolution failure is non-fatal
logger.info(f" ✓ DTM {res}m/px terminé ({t_dtm:.1f}s)") logger.info(f" ✓ DTM {res}m/px terminé ({t_dtm:.1f}s)")
dtm_rebuilt = True
self._write_dtm_method(basename, res_suffix)
# Process each resolution: visualizations + PDF # Process each resolution: visualizations + PDF
all_vis_results = {} all_vis_results = {}
@ -420,16 +486,18 @@ class LidarArchaeoPipeline:
vis_dir.mkdir(exist_ok=True) vis_dir.mkdir(exist_ok=True)
self.generate_all_visualizations(dtm_path, basename, actual_res, vis_dir=vis_dir) self.generate_all_visualizations(
dtm_path, basename, actual_res, vis_dir=vis_dir,
force_images=self.force or self.force_classify or dtm_rebuilt)
t_total = time.time() - t_start t_total = time.time() - t_start
logger.info(f"✓ {basename} terminé en {t_total:.1f}s") logger.info(f"✓ {basename} terminé en {t_total:.1f}s")
_file_filter.basename = None _file_filter.basename = None
return True return True
def process_all(self): def process_all(self, files=None):
"""Process all LAZ files in input directory.""" """Process all LAZ files in input directory (or an explicit list)."""
files = self.find_laz_files() files = files if files is not None else self.find_laz_files()
if not files: if not files:
logger.error("Aucun fichier LAZ/LAS trouvé !") logger.error("Aucun fichier LAZ/LAS trouvé !")
@ -466,7 +534,7 @@ class LidarArchaeoPipeline:
active_ids = self.gpu_ids if self.gpu_ids else available_gpu_ids() active_ids = self.gpu_ids if self.gpu_ids else available_gpu_ids()
resolutions_str = ','.join(str(r) for r in self.resolutions) resolutions_str = ','.join(str(r) for r in self.resolutions)
future_to_file = { future_to_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, active_ids[file_idx % len(active_ids)] if active_ids else None): 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.ign_classes, self.force_classify, self.keep_tif, self.bare_earth, self.quality, self.only_viz, self.skip_viz, self.output_format, active_ids[file_idx % len(active_ids)] if active_ids else None): laz_file
for file_idx, laz_file in enumerate(files) for file_idx, laz_file in enumerate(files)
} }
done = 0 done = 0
@ -545,7 +613,7 @@ class LidarArchaeoPipeline:
logger.warning(f" Note: Impossible de supprimer les fichiers temporaires: {e}") 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, quality=98, only_viz=None, skip_viz=None, output_format='avif', gpu_id=None): def _process_file_standalone(laz_file_str, input_dir, output_dir, resolution, force=False, ground_method='auto', ign_classes="sol", force_classify=False, keep_tif=False, bare_earth=False, quality=98, only_viz=None, skip_viz=None, output_format='avif', gpu_id=None):
"""Standalone function for multiprocessing — creates its own pipeline instance. """Standalone function for multiprocessing — creates its own pipeline instance.
Each worker gets its own temp directory to avoid file conflicts. Each worker gets its own temp directory to avoid file conflicts.
@ -572,7 +640,7 @@ def _process_file_standalone(laz_file_str, input_dir, output_dir, resolution, fo
worker_logger.addHandler(handler) worker_logger.addHandler(handler)
worker_logger.addFilter(_file_filter) 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, quality=quality, only_viz=only_viz, skip_viz=skip_viz, output_format=output_format) pipeline = LidarArchaeoPipeline(input_dir, output_dir, resolution=resolution, workers=1, force=force, ground_method=ground_method, ign_classes=ign_classes, force_classify=force_classify, keep_tif=keep_tif, bare_earth=bare_earth, quality=quality, only_viz=only_viz, skip_viz=skip_viz, output_format=output_format)
basename = _file_basename(laz_file_str) basename = _file_basename(laz_file_str)
pipeline.temp_dir = pipeline.output_dir / "temp" / basename pipeline.temp_dir = pipeline.output_dir / "temp" / basename
pipeline.temp_dir.mkdir(exist_ok=True) pipeline.temp_dir.mkdir(exist_ok=True)

View File

@ -818,8 +818,11 @@ def tif_to_crop(tif_file, vis_dir, resolution, keep_tif=False, quality=98, outpu
# Convert to RGB using colormap # Convert to RGB using colormap
if is_rgb_result: if is_rgb_result:
# RGB images are already in RGB # RGB images are already in RGB (uint8 depuis le TIF IGN, ou float 0-1)
rgb_data = (data * 255).astype(np.uint8) if data.dtype == np.uint8:
rgb_data = data
else:
rgb_data = (np.clip(data, 0, 1) * 255).astype(np.uint8)
else: else:
# Normalize data to 0-1 range for colormap # Normalize data to 0-1 range for colormap
cmap = plt.get_cmap(cmap_name) cmap = plt.get_cmap(cmap_name)

View File

@ -86,3 +86,20 @@ class TestSetupLogging:
fmt = logger.handlers[0].formatter._fmt fmt = logger.handlers[0].formatter._fmt
assert "%(filename)s" in fmt assert "%(filename)s" in fmt
assert "%(lineno)d" in fmt assert "%(lineno)d" in fmt
def test_rebuild_index_without_input_arg(tmp_path):
"""--rebuild-index fonctionne sans l'argument positionnel input.
Régression : input était obligatoire alors que --rebuild-index ne
l'utilise pas (erreur argparse « the following arguments are required »).
"""
import subprocess
r = subprocess.run(
[sys.executable, "-m", "lidar_pipeline", "--rebuild-index",
"-o", str(tmp_path)],
capture_output=True, text=True, timeout=180,
)
assert r.returncode == 0, r.stderr
assert "the following arguments are required" not in r.stderr

View File

@ -94,12 +94,127 @@ class TestCSFPipeline:
pipeline = json.loads(result) pipeline = json.loads(result)
csf_stage = [s for s in pipeline["pipeline"] if isinstance(s, dict) and s.get("type") == "filters.csf"][0] csf_stage = [s for s in pipeline["pipeline"] if isinstance(s, dict) and s.get("type") == "filters.csf"][0]
assert csf_stage["resolution"] == 0.5 assert csf_stage["resolution"] == 1.0 # cloth 1 m : ~4× plus rapide, MNT inchangé
assert csf_stage["rigidness"] == 3 assert csf_stage["rigidness"] == 3
assert csf_stage["smooth"] is True assert csf_stage["smooth"] is True
assert "hdiff" not in csf_stage # hdiff is not a valid PDAL CSF parameter assert "hdiff" not in csf_stage # hdiff is not a valid PDAL CSF parameter
class TestInterpolateHoles:
def test_fills_interior_hole_with_surface(self):
"""Large interior NaN hole is filled (no NaN left, value is plausible)."""
from lidar_pipeline.dtm import _interpolate_holes
# Linear-in-column surface z = 0.02 * x, with a large square hole in the middle.
x = np.arange(40, dtype=float) * 0.02
dtm = np.tile(x, (40, 1))
dtm[16:24, 16:24] = np.nan
filled, count = _interpolate_holes(dtm)
assert count == 64
assert not np.isnan(filled).any()
# Filled values stay within the surrounding z range (no wild extrapolation).
zmin, zmax = np.nanmin(dtm), np.nanmax(dtm)
hole_vals = filled[16:24, 16:24]
assert np.all(hole_vals >= zmin - 1e-6)
assert np.all(hole_vals <= zmax + 1e-6)
# A linear surface is interpolated near-exactly in the interior.
expected = np.tile(x[16:24], (8, 1))
assert np.allclose(hole_vals, expected, atol=0.02)
# Original valid cells are untouched.
valid = ~np.isnan(dtm)
assert np.allclose(filled[valid], dtm[valid])
def test_no_holes_returns_unchanged(self):
"""No NaN → returns same array and zero count."""
from lidar_pipeline.dtm import _interpolate_holes
dtm = np.arange(64, dtype=float).reshape(8, 8)
filled, count = _interpolate_holes(dtm)
assert count == 0
assert np.shares_memory(filled, dtm)
def test_all_nan_returns_unchanged(self):
"""No valid data → cannot interpolate, returns zeros-free NaN array."""
from lidar_pipeline.dtm import _interpolate_holes
dtm = np.full((8, 8), np.nan)
filled, count = _interpolate_holes(dtm)
assert count == 0
assert np.isnan(filled).all()
class TestMinReturnGrid:
def test_takes_lowest_return_per_cell(self, tmp_output_dir):
"""_min_return_grid rasterise le point le plus bas par cellule (pas la moyenne)."""
import laspy
from lidar_pipeline.dtm import _min_return_grid
out = tmp_output_dir / "pts.las"
hdr = laspy.LasHeader(version='1.2', point_format=0)
las = laspy.LasData(hdr)
# Grille 2x2 sur [0,2]x[0,2]. La cellule (0,0) porte deux points z=5 et
# z=2 (min=2, moyenne=3.5) ; (1,0) z=3 ; (0,1) z=4 ; (1,1) vide.
las.x = [0.2, 0.5, 1.2, 0.3]
las.y = [0.2, 0.3, 0.4, 1.5]
las.z = [5.0, 2.0, 3.0, 4.0]
las.write(str(out))
grid = _min_return_grid(out, 2, 2, (0.0, 0.0, 2.0, 2.0))
assert grid.shape == (2, 2)
assert int(np.isnan(grid).sum()) == 1
# Une seule valeur par cellule, et la cellule (0,0) vaut le MIN (2.0).
vals = sorted(float(v) for v in grid[~np.isnan(grid)])
assert vals == [2.0, 3.0, 4.0]
assert 3.5 not in vals
class TestBareEarth:
"""Le plancher « sol nu » ramène le DTM au retour le plus bas de chaque cellule."""
def _write_las(self, path, points):
"""points: list of (x, y, z). Écrit un LAS 1.2 format 0 aux bornes [0,2]x[0,2]."""
import laspy
hdr = laspy.LasHeader(version='1.2', point_format=0)
las = laspy.LasData(hdr)
las.x = [p[0] for p in points]
las.y = [p[1] for p in points]
las.z = [p[2] for p in points]
las.write(str(path))
return path
def _make_clouds(self, tmp_output_dir):
"""Grille 2x2 (res=1.0). Les points « coin » à 0.05/1.95 imposent l'étendue
[0.05,1.95] (laspy re-déduit les bornes de l'en-tête depuis les points).
Le sol (las_file) vaut z=10 partout. Le nuage complet (source_laz) a un
retour plus bas dans les cellules (0,0) -> 2 et (1,0) -> 5 ; les deux autres
cellules n'ont que z=10."""
ground = [(0.5, 0.5, 10.0), (1.5, 0.5, 10.0), (0.5, 1.5, 10.0), (1.5, 1.5, 10.0),
(0.05, 0.05, 10.0), (1.95, 1.95, 10.0)]
source = list(ground) + [(0.3, 0.3, 2.0), (1.3, 0.3, 5.0)]
las_file = self._write_las(tmp_output_dir / "ground.las", ground)
source_laz = self._write_las(tmp_output_dir / "source.las", source)
return las_file, source_laz
def _dtm_values(self, tmp_output_dir, bare_earth):
from lidar_pipeline.dtm import create_dtm_fast
import rasterio
las_file, source_laz = self._make_clouds(tmp_output_dir)
out = create_dtm_fast(las_file, "tile", tmp_output_dir, 1.0,
force=True, source_laz=source_laz, bare_earth=bare_earth)
assert out is not None
with rasterio.open(str(out)) as src:
arr = src.read(1).astype("float64")
return arr[~np.isnan(arr)]
def test_bare_earth_pulls_dtm_to_lowest_return(self, tmp_output_dir):
"""Avec bare_earth, le DTM descend aux retours les plus bas (2 et 5)."""
vals = sorted(float(v) for v in self._dtm_values(tmp_output_dir, bare_earth=True))
# Les cellules sans retour plus bas restent à 10 ; les deux autres descendent.
assert vals == [2.0, 5.0, 10.0, 10.0]
assert vals[0] == 2.0 and vals[1] == 5.0
def test_no_bare_earth_keeps_mean(self, tmp_output_dir):
"""Sans bare_earth, le DTM garde la moyenne des points sol (10 partout)."""
vals = sorted(float(v) for v in self._dtm_values(tmp_output_dir, bare_earth=False))
assert all(v == 10.0 for v in vals)
class TestDetectGroundMethod: class TestDetectGroundMethod:
def _make_mock_las(self, num_returns, z_values): def _make_mock_las(self, num_returns, z_values):
"""Create a mock laspy object with specified NumberOfReturns and z.""" """Create a mock laspy object with specified NumberOfReturns and z."""
@ -162,6 +277,73 @@ class TestDetectGroundMethod:
assert result == 'csf' assert result == 'csf'
class TestIGNPipeline:
def test_pipeline_keeps_supplier_classification(self):
"""create_ign_pipeline réutilise la pré-classification (classe 2) sans refiltrer."""
from lidar_pipeline.dtm import create_ign_pipeline
result = create_ign_pipeline("/input/a.laz", "/output/a_ground.las")
pipeline = json.loads(result)
stages = pipeline["pipeline"]
stage_types = [s.get("type") if isinstance(s, dict) else None for s in stages]
# Aucun algorithme de classification, pas de remise à zéro, pas de filtres de bruit
assert "filters.smrf" not in stage_types
assert "filters.csf" not in stage_types
assert "filters.assign" not in stage_types
assert "filters.elm" not in stage_types
assert "filters.outlier" not in stage_types
# Filtre ReturnNumber conservé + extraction des points classe 2
range_stages = [s for s in stages if isinstance(s, dict) and s.get("type") == "filters.range"]
assert any("ReturnNumber" in str(s.get("limits", "")) for s in range_stages)
assert any(s.get("limits") == "Classification[2:2]" for s in range_stages)
writer = [s for s in stages if isinstance(s, dict) and s.get("type") == "writers.las"][0]
assert writer["filename"] == "/output/a_ground.las"
class TestDetectIGN:
def _make_mock_las(self, classification, num_returns, z_values):
mock_las = MagicMock()
mock_las.classification = classification
mock_las.NumberOfReturns = np.array(num_returns)
mock_las.z = np.array(z_values)
mock_las.points = MagicMock()
mock_las.points.__len__ = lambda self: len(num_returns)
return mock_las
@patch('lidar_pipeline.dtm._read_with_pdal')
@patch('laspy.read')
def test_preclassified_returns_ign(self, mock_read, mock_pdal):
"""Fichier pré-classifié (majorité classe 2) → méthode IGN."""
from lidar_pipeline.dtm import detect_ground_method
n = 10000
num_returns = np.ones(n, dtype=int)
cls = np.zeros(n, dtype=np.uint8)
cls[int(n * 0.15):] = 2 # 85 % de points classe 2
z_values = np.random.normal(100, 5, n)
mock_read.return_value = self._make_mock_las(cls, num_returns, z_values)
assert detect_ground_method(Path("/data/input/test.laz")) == 'ign'
@patch('lidar_pipeline.dtm._read_with_pdal')
@patch('laspy.read')
def test_unclassified_falls_back_to_smrf_or_csf(self, mock_read, mock_pdal):
"""Sans classification exploitable → détection SMRF/CSF classique."""
from lidar_pipeline.dtm import detect_ground_method
n = 10000
num_returns = np.ones(n, dtype=int)
num_returns[:int(n * 0.6)] = 2 # 60 % multi-retours (forêt) → non urbain
cls = np.zeros(n, dtype=np.uint8) # aucun point classe 2
z_values = np.random.normal(100, 5, n)
mock_read.return_value = self._make_mock_las(cls, num_returns, z_values)
assert detect_ground_method(Path("/data/input/test.laz")) == 'smrf'
class TestClassifyGroundMethod: class TestClassifyGroundMethod:
@patch('lidar_pipeline.dtm.subprocess') @patch('lidar_pipeline.dtm.subprocess')
def test_classify_ground_auto_calls_detect(self, mock_subprocess): def test_classify_ground_auto_calls_detect(self, mock_subprocess):
@ -212,3 +394,157 @@ class TestClassifyGroundMethod:
pipeline = json.loads(pipeline_file.read_text()) pipeline = json.loads(pipeline_file.read_text())
stage_types = [s.get("type") if isinstance(s, dict) else None for s in pipeline["pipeline"]] stage_types = [s.get("type") if isinstance(s, dict) else None for s in pipeline["pipeline"]]
assert "filters.csf" in stage_types assert "filters.csf" in stage_types
class TestParseIgnClasses:
def test_default_sol(self):
"""'sol' → code 2 seul."""
from lidar_pipeline.dtm import parse_ign_classes
assert parse_ign_classes("sol") == [2]
def test_names_sorted_dedup(self):
"""Noms acceptés (EN/FR), triés et dédupliqués."""
from lidar_pipeline.dtm import parse_ign_classes
assert parse_ign_classes("sol,unclassified") == [1, 2]
assert parse_ign_classes("unclassified,sol") == [1, 2]
assert parse_ign_classes("non-classe") == [1]
assert parse_ign_classes("sol,2") == [2]
def test_numeric_codes(self):
"""Codes LAS directs, triés."""
from lidar_pipeline.dtm import parse_ign_classes
assert parse_ign_classes("2,1") == [1, 2]
assert parse_ign_classes("66") == [66]
def test_invalid_raises(self):
"""Nom inconnu, code hors bornes ou liste vide → ValueError."""
from lidar_pipeline.dtm import parse_ign_classes
with pytest.raises(ValueError):
parse_ign_classes("foo")
with pytest.raises(ValueError):
parse_ign_classes("300")
with pytest.raises(ValueError):
parse_ign_classes("")
def test_method_label(self):
"""'ign' seul pour le sol, combinaison encodée sinon (invalidation cache)."""
from lidar_pipeline.dtm import ign_method_label
assert ign_method_label([2]) == "ign"
assert ign_method_label([1, 2]) == "ign_1_2"
assert ign_method_label([2, 1]) == "ign_1_2"
class TestIGNPipelineMultiClasses:
def test_multi_class_limits(self):
"""Plusieurs classes → plages OU logiques sur Classification."""
from lidar_pipeline.dtm import _create_ground_pipeline
result = _create_ground_pipeline("/input/a.laz", "/output/a_ground.las",
'ign', ign_codes=[1, 2])
pipeline = json.loads(result)
range_stages = [s for s in pipeline["pipeline"]
if isinstance(s, dict) and s.get("type") == "filters.range"]
limits = [str(s.get("limits", "")) for s in range_stages]
assert any("Classification[1:1]" in l and "Classification[2:2]" in l
for l in limits)
def test_default_sol_only(self):
"""Sans ign_codes, la voie IGN reste sol seul (2) — rétrocompatible."""
from lidar_pipeline.dtm import _create_ground_pipeline
result = _create_ground_pipeline("/input/a.laz", "/output/a_ground.las", 'ign')
pipeline = json.loads(result)
range_stages = [s for s in pipeline["pipeline"]
if isinstance(s, dict) and s.get("type") == "filters.range"]
limits = [str(s.get("limits", "")) for s in range_stages]
assert any("Classification[2:2]" in l and "Classification[1:1]" not in l
for l in limits)
class TestClassifyGroundIgnClasses:
@patch('lidar_pipeline.dtm.subprocess')
def test_ign_classes_encoded_in_filenames(self, mock_subprocess):
"""--ign-classes sol,unclassified → fichiers ign_1_2 + filtre multi-classes."""
import tempfile
from lidar_pipeline.dtm import classify_ground
mock_subprocess.run.return_value = MagicMock(returncode=0)
with tempfile.TemporaryDirectory() as tmpdir:
tmpdir = Path(tmpdir)
classify_ground(Path("/data/input/test.laz"), tmpdir,
method='ign', ign_classes="sol,unclassified")
pipeline_file = tmpdir / "pipeline_ign_1_2.json"
assert pipeline_file.exists()
pipeline = json.loads(pipeline_file.read_text())
limits = [str(s.get("limits", "")) for s in pipeline["pipeline"]
if isinstance(s, dict) and s.get("type") == "filters.range"]
assert any("Classification[1:1]" in l and "Classification[2:2]" in l
for l in limits)
@patch('lidar_pipeline.dtm.subprocess')
def test_ign_default_label_unchanged(self, mock_subprocess):
"""--ign-classes sol (défaut) → noms 'ign' inchangés (cache préservé)."""
import tempfile
from lidar_pipeline.dtm import classify_ground
mock_subprocess.run.return_value = MagicMock(returncode=0)
with tempfile.TemporaryDirectory() as tmpdir:
tmpdir = Path(tmpdir)
classify_ground(Path("/data/input/test.laz"), tmpdir, method='ign')
assert (tmpdir / "pipeline_ign.json").exists()
assert not (tmpdir / "pipeline_ign_1_2.json").exists()
class TestPureDtm:
"""Mode pur (classification IGN) : aucune retouche, trous en nodata."""
def _write_las(self, path, points):
import laspy
hdr = laspy.LasHeader(version='1.2', point_format=0)
las = laspy.LasData(hdr)
las.x = [p[0] for p in points]
las.y = [p[1] for p in points]
las.z = [p[2] for p in points]
las.write(str(path))
return path
def _make_clouds(self, tmp_output_dir):
"""Grille 2x2 (res=1.0). Sol sur 3 cellules (z=10), trou en (1,1).
Le nuage complet a un retour plus bas (z=7) dans le trou."""
corners = [(0.05, 0.05, 10.0), (1.95, 0.05, 10.0), (0.05, 1.95, 10.0)]
ground = [(0.5, 0.5, 10.0), (1.5, 0.5, 10.0), (0.5, 1.5, 10.0)] + corners
source = list(ground) + [(1.5, 1.5, 7.0)]
las_file = self._write_las(tmp_output_dir / "ground_pure.las", ground)
source_laz = self._write_las(tmp_output_dir / "source_pure.las", source)
return las_file, source_laz
def _dtm_array(self, tmp_output_dir, pure):
from lidar_pipeline.dtm import create_dtm_fast
import rasterio
las_file, source_laz = self._make_clouds(tmp_output_dir)
out = create_dtm_fast(las_file, "tile_pure", tmp_output_dir, 1.0,
force=True, source_laz=source_laz, pure=pure)
assert out is not None
with rasterio.open(str(out)) as src:
return src.read(1).astype("float64")
def test_pure_fills_holes_without_floor(self, tmp_output_dir):
"""pur=True : trous comblés par interpolation, sans plancher à 7.
Le comblement est actif dans tous les modes (comportement
historique) ; « pur » ne désactive que l'abaissement au retour
le plus bas.
"""
arr = self._dtm_array(tmp_output_dir, pure=True)
assert int(np.isnan(arr).sum()) == 0
vals = sorted(float(v) for v in arr.flatten())
assert vals == [10.0, 10.0, 10.0, 10.0]
def test_not_pure_fills_holes(self, tmp_output_dir):
"""pur=False : le trou est comblé (comportement historique conservé)."""
arr = self._dtm_array(tmp_output_dir, pure=False)
assert int(np.isnan(arr).sum()) == 0
vals = sorted(float(v) for v in arr.flatten())
assert len(vals) == 4
assert vals[-1] == 10.0

View File

@ -0,0 +1,98 @@
"""Tests du téléchargement des dalles LiDAR HD de l'IGN (fetch_ign)."""
def test_parse_tile_specs():
"""Accepte 'col,row', 'col:row' et ignore les espaces."""
from lidar_pipeline.fetch_ign import parse_tile_specs
assert parse_tile_specs(["1055,6882"]) == [(1055, 6882)]
assert parse_tile_specs(["1055:6883", " 651 , 6630 "]) == [(1055, 6883), (651, 6630)]
def test_parse_tile_specs_rejects_invalid():
"""Les spécifications mal formées lèvent une erreur explicite."""
import pytest
from lidar_pipeline.fetch_ign import parse_tile_specs
with pytest.raises(ValueError):
parse_tile_specs(["1055"])
with pytest.raises(ValueError):
parse_tile_specs(["abc,def"])
with pytest.raises(ValueError):
parse_tile_specs(["1055,6882,9999"])
def test_tile_filename_pads_coordinates():
"""Le nom de fichier DALLE utilise des coordonnées à 4 chiffres."""
from lidar_pipeline.fetch_ign import tile_filename
assert tile_filename(1055, 6882) == "LHD_FXX_1055_6882_PTS_LAMB93_IGN69.copc.laz"
assert tile_filename(651, 6630) == "LHD_FXX_0651_6630_PTS_LAMB93_IGN69.copc.laz"
def test_match_feature_uses_coordonnees_nw():
"""La correspondance se fait sur lidarhd:coordonnees_NW (format 0816-6847)."""
from lidar_pipeline.fetch_ign import match_feature
features = [
{"id": "LHD_FXX_1054_6882_PTS_LAMB93_IGN69_NE",
"properties": {"lidarhd:coordonnees_NW": "1054-6882"},
"assets": {"data": {"href": "https://example.org/a.copc.laz"}}},
{"id": "LHD_FXX_1055_6882_PTS_LAMB93_IGN69_NE",
"properties": {"lidarhd:coordonnees_NW": "1055-6882"},
"assets": {"data": {"href": "https://example.org/b.copc.laz"}}},
]
assert match_feature(features, 1055, 6882) is features[1]
assert match_feature(features, 9999, 9999) is None
def test_find_tile_url_matches_and_returns_href(monkeypatch):
"""find_tile_url interroge le STAC et retourne l'href de l'asset data."""
from lidar_pipeline import fetch_ign
class FakeResponse:
def __init__(self, payload):
self._payload = payload.encode("utf-8")
def read(self):
return self._payload
def __enter__(self):
return self
def __exit__(self, *args):
return False
payload = ('{"features": [{"properties": {"lidarhd:coordonnees_NW": "1055-6882"},'
'"assets": {"data": {"href": "https://data.geopf.fr/x.copc.laz"}}}]}')
captured = {}
def fake_urlopen(req, timeout=None):
captured["url"] = req.full_url
return FakeResponse(payload)
monkeypatch.setattr(fetch_ign.urllib.request, "urlopen", fake_urlopen)
url = fetch_ign.find_tile_url(1055, 6882)
assert url == "https://data.geopf.fr/x.copc.laz"
assert "api.stac.teledetection.fr" in captured["url"]
assert "bbox=" in captured["url"]
def test_fetch_tiles_skips_existing_and_generated(tmp_path, monkeypatch):
"""Pas de téléchargement si le LAZ existe ou si les visualisations existent."""
from lidar_pipeline import fetch_ign
input_dir = tmp_path / "input"
input_dir.mkdir()
output_dir = tmp_path / "output"
vis_dir = output_dir / "visualisations" / "LHD_FXX_1054_6882_PTS_LAMB93_IGN69"
vis_dir.mkdir(parents=True)
# 1054,6882 : visualisations déjà générées
# 1055,6882 : LAZ déjà présent
existing = input_dir / fetch_ign.tile_filename(1055, 6882)
existing.write_bytes(b"naze")
def fail_download(*args, **kwargs):
raise AssertionError("ne doit pas être appelé")
monkeypatch.setattr(fetch_ign, "find_tile_url", fail_download)
result = fetch_ign.fetch_tiles(input_dir, [(1054, 6882), (1055, 6882)],
output_dir=output_dir)
assert result == []

View File

@ -60,7 +60,11 @@ def test_compute_bbox_empty():
def _make_fake_viz_dir(vis_dir, basename, col, row, viz_keys=('hillshade_multi', 'svf'), ext='webp', res_suffix=''): def _make_fake_viz_dir(vis_dir, basename, col, row, viz_keys=('hillshade_multi', 'svf'), ext='webp', res_suffix=''):
"""Crée un faux dossier de visualisations avec de petites images.""" """Crée un faux dossier de visualisations avec de petites images.
Le suffixe de résolution apparaît seulement dans le nom du dossier (miroir
du pipeline : les fichiers restent préfixés par le basename nu).
"""
from PIL import Image as PILImage from PIL import Image as PILImage
import numpy as np import numpy as np
@ -70,7 +74,7 @@ def _make_fake_viz_dir(vis_dir, basename, col, row, viz_keys=('hillshade_multi',
for v in viz_keys: for v in viz_keys:
arr = np.random.randint(0, 255, (50, 50, 3), dtype=np.uint8) arr = np.random.randint(0, 255, (50, 50, 3), dtype=np.uint8)
img = PILImage.fromarray(arr) img = PILImage.fromarray(arr)
fname = f"LHD_FXX_{col}_{row}_PTS_LAMB93_IGN69{res_suffix}_{v}.{ext}" fname = f"LHD_FXX_{col}_{row}_PTS_LAMB93_IGN69_{v}.{ext}"
img.save(str(tile_dir / fname), format='WEBP', quality=80) img.save(str(tile_dir / fname), format='WEBP', quality=80)
return tile_dir return tile_dir
@ -124,6 +128,85 @@ def test_scan_tiles_multi_resolution(tmp_path):
assert resolutions == [0.2, 0.5] assert resolutions == [0.2, 0.5]
def test_res_suffix_str():
"""Le suffixe de résolution reflète le nommage du pipeline (miroir)."""
from lidar_pipeline.index import _res_suffix_str
assert _res_suffix_str(0.5) == ''
assert _res_suffix_str(0.2) == '_r0p2'
def test_collect_tile_metadata(tmp_path):
"""Les métadonnées lisent la méthode DTM et les dates/tailles des viz."""
import os
from datetime import datetime
from lidar_pipeline.index import _collect_tile_metadata
basename = "LHD_FXX_1000_6881_PTS_LAMB93_IGN69"
tile_dir = tmp_path / "visualisations" / basename
tile_dir.mkdir(parents=True)
viz_file = tile_dir / f"{basename}_hillshade_multi.webp"
viz_file.write_bytes(b"fake")
dtm_dir = tmp_path / "DTM"
dtm_dir.mkdir()
method_file = dtm_dir / f"{basename}_dtm_method.txt"
method_file.write_text("ign", encoding="utf-8")
# Dates déterministes : method.txt plus ancien que la viz
os.utime(method_file, (1600000000, 1600000000))
os.utime(viz_file, (1700000000, 1700000000))
fmt = lambda ts: datetime.fromtimestamp(ts).strftime('%Y-%m-%d %H:%M')
tile = {
'basename': basename, 'resolution': 0.5,
'dir_path': str(tile_dir),
'viz': {'hillshade_multi': {'filename': viz_file.name, 'ext': 'webp'}},
}
meta = _collect_tile_metadata(tile, dtm_dir)
assert meta['method'] == 'ign'
assert meta['generated'] == fmt(1600000000)
assert meta['viz']['hillshade_multi']['size'] == 4
assert meta['viz']['hillshade_multi']['date'] == fmt(1700000000)
def test_collect_tile_metadata_resolution_suffix(tmp_path):
"""Une tuile 0,2 m lit son sidecar _dtm_r0p2_method.txt dédié."""
from lidar_pipeline.index import _collect_tile_metadata
basename = "LHD_FXX_1000_6881_PTS_LAMB93_IGN69"
tile_dir = tmp_path / "visualisations" / (basename + "_r0p2")
tile_dir.mkdir(parents=True)
dtm_dir = tmp_path / "DTM"
dtm_dir.mkdir()
(dtm_dir / f"{basename}_dtm_r0p2_method.txt").write_text("smrf", encoding="utf-8")
tile = {'basename': basename, 'resolution': 0.2,
'dir_path': str(tile_dir),
'viz': {}}
meta = _collect_tile_metadata(tile, dtm_dir)
assert meta['method'] == 'smrf'
# La date vient du sidecar (écrit juste après la création du DTM)
assert meta['generated'] is not None
assert meta['viz'] == {}
def test_collect_tile_metadata_fallback_date(tmp_path):
"""Sans sidecar DTM, la date de génération remonte au plus ancien fichier viz."""
from lidar_pipeline.index import _collect_tile_metadata
basename = "LHD_FXX_1000_6881_PTS_LAMB93_IGN69"
tile_dir = tmp_path / "visualisations" / basename
tile_dir.mkdir(parents=True)
f = tile_dir / f"{basename}_svf.webp"
f.write_bytes(b"x")
tile = {'basename': basename, 'resolution': 0.5,
'dir_path': str(tile_dir),
'viz': {'svf': {'filename': f.name, 'ext': 'webp'}}}
meta = _collect_tile_metadata(tile, tmp_path / "DTM")
assert meta['method'] is None
assert meta['generated'] is not None
def test_build_index_generates_html(tmp_path): def test_build_index_generates_html(tmp_path):
"""build_index génère index.html et les vignettes.""" """build_index génère index.html et les vignettes."""
from lidar_pipeline.index import build_index from lidar_pipeline.index import build_index
@ -142,11 +225,18 @@ def test_build_index_generates_html(tmp_path):
content = html_path.read_text(encoding='utf-8') content = html_path.read_text(encoding='utf-8')
# Vérifie la présence des éléments clés # Vérifie la présence des éléments clés
assert "Carte continue LiDAR" in content assert "Carte LiDAR" in content
assert "LHD_FXX_1000_6881" in content assert "LHD_FXX_1000_6881" in content
assert "LHD_FXX_1001_6881" in content assert "LHD_FXX_1001_6881" in content
# Vérifie que le JSON intégré est valide # Vérifie que le JSON intégré est valide
assert "const DATA" in content assert "const TILES" in content
# Vérifie les assets de l'interface (CSS/JS séparés)
assets = output_dir / "assets"
assert (assets / "app.css").read_text(encoding='utf-8').startswith('/*')
app_js = (assets / "app.js").read_text(encoding='utf-8')
assert "Couches" in app_js or "layers" in app_js
assert 'assets/app.css' in content
assert 'assets/app.js' in content
# Vérifie les vignettes générées # Vérifie les vignettes générées
thumb_dir = output_dir / "index_thumbs" thumb_dir = output_dir / "index_thumbs"
assert thumb_dir.is_dir() assert thumb_dir.is_dir()
@ -154,6 +244,68 @@ def test_build_index_generates_html(tmp_path):
assert len(thumbs) >= 2 # au moins hillshade pour chaque tuile assert len(thumbs) >= 2 # au moins hillshade pour chaque tuile
def test_build_index_regenerates_stale_thumbnails(tmp_path):
"""Une tuile recalculée (source plus récente) régénère sa vignette."""
import os
import time
import numpy as np
from PIL import Image as PILImage
from lidar_pipeline.index import build_index
output_dir = tmp_path / "output"
vis_dir = output_dir / "visualisations"
vis_dir.mkdir(parents=True)
tile_dir = _make_fake_viz_dir(vis_dir, "a", 1000, 6881, ('hillshade_multi',))
assert build_index(output_dir) is not None
thumb_path = output_dir / "index_thumbs" / "LHD_FXX_1000_6881_PTS_LAMB93_IGN69_hillshade_multi.jpg"
assert thumb_path.exists()
m1 = thumb_path.stat().st_mtime
# Recalcul de la tuile : source réécrite avec une mtime plus récente
src = tile_dir / "LHD_FXX_1000_6881_PTS_LAMB93_IGN69_hillshade_multi.webp"
arr = np.random.randint(0, 255, (50, 50, 3), dtype=np.uint8)
PILImage.fromarray(arr).save(str(src), format='WEBP', quality=80)
os.utime(src, (m1 + 5, m1 + 5))
assert build_index(output_dir) is not None
m2 = thumb_path.stat().st_mtime
assert m2 > m1 # vignette régénérée
# Source non modifiée depuis → pas de régénération inutile
os.utime(src, (time.time() - 10, time.time() - 10))
assert build_index(output_dir) is not None
assert thumb_path.stat().st_mtime == m2
def test_build_subtiles_regenerates_stale_crops(tmp_path):
"""Une dalle 0,2 m recalculée régénère ses sous-tuiles (par visualisation)."""
import os
from lidar_pipeline.index import build_index
output_dir = tmp_path / "output"
vis_dir = output_dir / "visualisations"
vis_dir.mkdir(parents=True)
tile_dir = _make_fake_viz_dir(vis_dir, "a", 1000, 6881,
('hillshade_multi', 'aspect'), res_suffix='_r0p2')
assert build_index(output_dir) is not None
sub_dir = output_dir / "index_subtiles"
hill_avif = sub_dir / "LHD_FXX_1000_6881_PTS_LAMB93_IGN69_r0p2_hillshade_multi_0_0.avif"
aspect_avif = sub_dir / "LHD_FXX_1000_6881_PTS_LAMB93_IGN69_r0p2_aspect_0_0.avif"
assert hill_avif.exists() and aspect_avif.exists()
m_hill_1 = hill_avif.stat().st_mtime
m_aspect_1 = aspect_avif.stat().st_mtime
# Recalcul : seule la source hillshade est plus récente
src = tile_dir / "LHD_FXX_1000_6881_PTS_LAMB93_IGN69_hillshade_multi.webp"
os.utime(src, (m_hill_1 + 5, m_hill_1 + 5))
assert build_index(output_dir) is not None
assert hill_avif.stat().st_mtime > m_hill_1 # sous-tuiles hillshade régénérées
assert aspect_avif.stat().st_mtime == m_aspect_1 # aspect intact
def test_build_index_empty_returns_none(tmp_path): def test_build_index_empty_returns_none(tmp_path):
"""Aucune tuile → build_index retourne None sans crash.""" """Aucune tuile → build_index retourne None sans crash."""
from lidar_pipeline.index import build_index from lidar_pipeline.index import build_index
@ -175,26 +327,66 @@ def test_build_index_embeds_valid_json(tmp_path):
build_index(output_dir) build_index(output_dir)
content = (output_dir / "index.html").read_text(encoding='utf-8') content = (output_dir / "index.html").read_text(encoding='utf-8')
# Extrait le JSON entre "const DATA = " et ";" # Extrait le JSON entre "const TILES = " et la fin de déclaration
start = content.index("const DATA = ") + len("const DATA = ") start = content.index("const TILES = ") + len("const TILES = ")
# Trouve le ; de fin de déclaration # Trouve le ; de fin de déclaration
depth = 0 depth = 0
end = start end = start
for i, ch in enumerate(content[start:], start): for i, ch in enumerate(content[start:], start):
if ch == '{': if ch in ('{', '['):
depth += 1 depth += 1
elif ch == '}': elif ch in ('}', ']'):
depth -= 1 depth -= 1
if depth == 0: if depth == 0:
end = i + 1 end = i + 1
break break
data = json.loads(content[start:end]) data = json.loads(content[start:end])
assert 'tiles' in data assert len(data) > 0
assert 'bbox' in data assert data[0]['col'] == 1000
assert 'vizList' in data assert data[0]['row'] == 6881
assert len(data['tiles']) == 1
assert data['tiles'][0]['col'] == 1000
assert data['tiles'][0]['row'] == 6881 def test_attach_gps_bounds():
"""attach_gps_bounds ajoute des bounds GPS ordonnées (France métropolitaine)."""
from lidar_pipeline.index import attach_gps_bounds
tiles = [{'col': 1000, 'row': 6881}, {'col': 1042, 'row': 6900}]
attach_gps_bounds(tiles)
for t in tiles:
assert 'bounds' in t
(lat_s, lon_w), (lat_n, lon_e) = t['bounds']
assert lat_n > lat_s
assert lon_e > lon_w
# France métropolitaine
assert 41 < lat_s < 51
assert -5 < lon_w < 10
def test_attach_gps_bounds_row_is_north_edge():
"""Le numéro de ligne du fichier = bord NORD (convention LiDAR HD IGN).
Vérifié sur les bounds des DTM : X ∈ [col, col+1] km, Y ∈ [row-1, row] km.
La régression historique plaçait Y ∈ [row, row+1] (1 km trop au nord).
"""
from rasterio.warp import transform as warp_transform
from lidar_pipeline.index import attach_gps_bounds
col, row = 1054, 6882
tiles = [{'col': col, 'row': row}]
attach_gps_bounds(tiles)
corners = tiles[0]['corners']
# Référence exacte de la vraie cellule : SW, SE, NE, NW
xs = [col * 1000, (col + 1) * 1000, (col + 1) * 1000, col * 1000]
ys = [(row - 1) * 1000, (row - 1) * 1000, row * 1000, row * 1000]
lons, lats = warp_transform('EPSG:2154', 'EPSG:4326', xs, ys)
for k in range(4):
assert abs(corners[k][0] - lats[k]) < 1e-9
assert abs(corners[k][1] - lons[k]) < 1e-9
# L'ancienne convention (row = bord sud) serait décalée d'environ 1 km
lat_n = max(c[0] for c in corners)
assert abs(lat_n - max(lats)) < 1e-9 # bord nord = Y = row×1000
def test_pick_display_viz_prefers_hillshade(): def test_pick_display_viz_prefers_hillshade():
@ -203,3 +395,30 @@ def test_pick_display_viz_prefers_hillshade():
assert _pick_display_viz(['svf', 'hillshade_multi', 'slope']) == 'hillshade_multi' assert _pick_display_viz(['svf', 'hillshade_multi', 'slope']) == 'hillshade_multi'
assert _pick_display_viz(['svf', 'slope']) == 'svf' assert _pick_display_viz(['svf', 'slope']) == 'svf'
assert _pick_display_viz(['topo']) == 'topo' assert _pick_display_viz(['topo']) == 'topo'
def test_subdivision_k():
"""0,5 m/px (2000 px) reste entier ; 0,2 m/px (5000 px) est découpé en 2×2."""
from lidar_pipeline.index import _subdivision_k
assert _subdivision_k(0.5) == 1
assert _subdivision_k(0.2) == 2
assert _subdivision_k(1.0) == 1
def test_subtile_corners_grid():
"""Les sous-tuiles reconstruisent exactement la grille de la dalle."""
from lidar_pipeline.index import _subtile_corners
corners = [[10.0, 2.0], [10.0, 3.0], [11.0, 3.0], [11.0, 2.0]] # SW SE NE NW
k = 2
sw_quad = _subtile_corners(corners, 0, 0, k) # quadrant sud-ouest
ne_quad = _subtile_corners(corners, 1, 1, k) # quadrant nord-est
# Le quadrant SW partage le coin SW de la dalle
assert sw_quad[0] == corners[0]
# Le quadrant NE partage le coin NE de la dalle
assert ne_quad[2] == corners[2]
# Le quadrant SW a son coin NE au centre de la dalle
assert sw_quad[2] == [10.5, 2.5]
# Adjacence : bord est du SW = bord ouest du SE (0,0)-(1,0)
se_quad = _subtile_corners(corners, 1, 0, k)
assert sw_quad[1] == se_quad[0]
assert sw_quad[2] == se_quad[3]

View File

@ -71,3 +71,98 @@ class TestLidarArchaeoPipeline:
assert "test.laz" in names assert "test.laz" in names
assert "other.las" in names assert "other.las" in names
assert "readme.txt" not in names assert "readme.txt" not in names
class TestDtmMethodSidecar:
"""Méthode de classification enregistrée à côté du DTM (invalidation du cache)."""
def test_missing_sidecar_matches(self, tmp_path):
from lidar_pipeline.pipeline import LidarArchaeoPipeline
input_dir = tmp_path / "input"
input_dir.mkdir()
pipeline = LidarArchaeoPipeline(str(input_dir), str(tmp_path / "output"), ground_method='csf')
# Aucun sidecar écrit → cache conservé (considéré compatible).
assert pipeline._dtm_method_matches("tileA", "") is True
def test_matching_method(self, tmp_path):
from lidar_pipeline.pipeline import LidarArchaeoPipeline
input_dir = tmp_path / "input"
input_dir.mkdir()
pipeline = LidarArchaeoPipeline(str(input_dir), str(tmp_path / "output"), ground_method='csf')
pipeline._write_dtm_method("tileA", "")
assert pipeline._dtm_method_matches("tileA", "") is True
def test_different_method_invalidates_cache(self, tmp_path):
from lidar_pipeline.pipeline import LidarArchaeoPipeline
input_dir = tmp_path / "input"
input_dir.mkdir()
out = str(tmp_path / "output")
LidarArchaeoPipeline(str(input_dir), out, ground_method='ign')._write_dtm_method("tileA", "")
csf = LidarArchaeoPipeline(str(input_dir), out, ground_method='csf')
assert csf._dtm_method_matches("tileA", "") is False
def test_write_dtm_method(self, tmp_path):
from lidar_pipeline.pipeline import LidarArchaeoPipeline
input_dir = tmp_path / "input"
input_dir.mkdir()
pipeline = LidarArchaeoPipeline(str(input_dir), str(tmp_path / "output"), ground_method='smrf')
pipeline._write_dtm_method("tileA", "_r0p2")
sidecar = tmp_path / "output" / "DTM" / "tileA_dtm_r0p2_method.txt"
assert sidecar.exists()
assert sidecar.read_text(encoding="utf-8").strip() == "smrf"
assert pipeline._dtm_method_name("tileA", "_r0p2") == "smrf"
# Le sidecar est un fichier .txt : il ne gêne pas la recherche des DTM .tif.
dtm = tmp_path / "output" / "DTM" / "tileA_dtm_r0p2.tif"
dtm.touch()
assert [p.name for p in (tmp_path / "output" / "DTM").glob("*.tif")] == ["tileA_dtm_r0p2.tif"]
def test_force_images_regenerates_existing(self, tmp_path):
from lidar_pipeline.pipeline import LidarArchaeoPipeline
input_dir = tmp_path / "input"
input_dir.mkdir()
pipeline = LidarArchaeoPipeline(str(input_dir), str(tmp_path / "output"), output_format='avif')
calls = []
def fake_ortho(dem_file, basename, vis_dir, resolution):
calls.append(basename)
return vis_dir / f"{basename}_ortho.avif"
pipeline.viz_steps = [('ortho', fake_ortho)]
vis_dir = tmp_path / "output" / "visualisations" / "tileA"
vis_dir.mkdir(parents=True)
(vis_dir / "tileA_ortho.avif").touch()
dtm = tmp_path / "dtm.tif"
# Image existante, pas de force → ignorée (pas de régénération).
pipeline.generate_all_visualizations(dtm, "tileA", resolution=0.5, vis_dir=vis_dir, force_images=False)
assert calls == []
# Image existante, force_images=True → régénérée.
calls.clear()
pipeline.generate_all_visualizations(dtm, "tileA", resolution=0.5, vis_dir=vis_dir, force_images=True)
assert calls == ["tileA"]
class TestEffectiveGroundMethod:
def test_ign_label_encodes_classes(self):
"""Les classes IGN sont encodées dans l'étiquette de cache (reclassification)."""
from lidar_pipeline.pipeline import LidarArchaeoPipeline
import tempfile
with tempfile.TemporaryDirectory() as tmpdir:
p = LidarArchaeoPipeline(tmpdir, tmpdir, ground_method='ign',
ign_classes="sol,unclassified")
assert p._effective_ground_method() == "ign_1_2"
def test_ign_default_label(self):
from lidar_pipeline.pipeline import LidarArchaeoPipeline
import tempfile
with tempfile.TemporaryDirectory() as tmpdir:
p = LidarArchaeoPipeline(tmpdir, tmpdir, ground_method='ign')
assert p._effective_ground_method() == "ign"
def test_other_methods_unchanged(self):
from lidar_pipeline.pipeline import LidarArchaeoPipeline
import tempfile
with tempfile.TemporaryDirectory() as tmpdir:
p = LidarArchaeoPipeline(tmpdir, tmpdir, ground_method='smrf',
ign_classes="sol,unclassified")
assert p._effective_ground_method() == "smrf"

View File

@ -96,3 +96,57 @@ class TestApplyColormap:
result = tif_to_png(tif_file, tmp_path, 5.0) result = tif_to_png(tif_file, tmp_path, 5.0)
assert result is not None assert result is not None
assert result.exists() assert result.exists()
class TestTifToCrop:
"""Conversion TIF → dalle cartographique (tif_to_crop)."""
@staticmethod
def _write_named_tif(tmp_path, name, arr):
transform = from_bounds(660000, 6700000, 661000, 6701000, arr.shape[1], arr.shape[0])
tif_file = tmp_path / name
with rasterio.open(
tif_file, 'w', driver='GTiff', height=arr.shape[0], width=arr.shape[1],
count=1, dtype='float32', crs='EPSG:2154', transform=transform,
nodata=float('nan'), compress='lzw'
) as dst:
dst.write(arr.astype('float32'), 1)
return tif_file
def test_nodata_renders_black(self, tmp_path):
"""Le nodata restant est rendu en noir (comportement historique).
Les trous du MNT sont comblés en amont (interpolation dans
create_dtm_fast, tous modes) ; ce qui reste en nodata doit rester
visible en noir sur la dalle plutôt qu'inventé au rendu.
"""
from PIL import Image as PILImage
from lidar_pipeline.rendering import tif_to_crop
data = np.random.default_rng(7).normal(50, 10, (40, 40)).astype(np.float32)
data[15:25, 15:25] = np.nan
tif_file = self._write_named_tif(tmp_path, "LHD_test_slope.tif", data)
# WebP lossless : l'encodeur AVIF de l'image « saigne » légèrement les
# bords du noir même en lossless — on teste la logique nodata→noir,
# pas les artefacts du codec.
out = tif_to_crop(tif_file, tmp_path, 5.0, keep_tif=True,
quality=100, output_format='webp')
assert out is not None and out.exists()
rgb = np.asarray(PILImage.open(str(out)).convert('RGB'))
hole = rgb[15:25, 15:25, :]
assert np.all(hole == 0), "le nodata doit être rendu en noir"
def test_without_nodata(self, tmp_path):
"""Un TIF sans nodata est converti sans crash, taille préservée."""
from PIL import Image as PILImage
from lidar_pipeline.rendering import tif_to_crop
data = np.random.default_rng(7).normal(50, 10, (40, 40)).astype(np.float32)
tif_file = self._write_named_tif(tmp_path, "LHD_test_slope.tif", data)
out = tif_to_crop(tif_file, tmp_path, 5.0)
assert out is not None and out.exists()
img = PILImage.open(str(out))
assert img.size == (40, 40)

View File

@ -157,3 +157,61 @@ class TestRayTrace:
) )
assert pos.shape == (4, 2, rows, cols) assert pos.shape == (4, 2, rows, cols)
assert neg.shape == (4, 2, rows, cols) assert neg.shape == (4, 2, rows, cols)
class TestNodataPreserved:
"""Nodata préservé dans les rendus (comportement historique).
Les trous du MNT sont comblés en amont (create_dtm_fast, tous modes) ;
si un nodata subsiste malgré tout, hillshade/slope/aspect le restituent
(rendu noir en carte) au lieu d'inventer des valeurs interpolées.
"""
@staticmethod
def _dem_with_hole(synthetic_dem, tmp_path):
import rasterio
with rasterio.open(synthetic_dem) as src:
arr = src.read(1).copy()
profile = src.profile.copy()
arr[80:120, 80:120] = np.nan
dem_hole = tmp_path / "dem_hole.tif"
profile.update(dtype='float32', nodata=float('nan'))
with rasterio.open(dem_hole, 'w', **profile) as dst:
dst.write(arr.astype('float32'), 1)
return dem_hole
def test_aspect_solo_preserves_nodata(self, synthetic_dem, tmp_path):
from lidar_pipeline.visualizations import generate_aspect
dem_hole = self._dem_with_hole(synthetic_dem, tmp_path)
out = generate_aspect(dem_hole, "solo", tmp_path, 5.0)
assert out is not None and out.exists()
import rasterio
with rasterio.open(out) as src:
data = src.read(1)
assert np.isnan(data[80:120, 80:120]).all(), "le trou doit rester en nodata"
# Le gradient au bord du trou propage NaN sur un anneau de 1 px :
# on vérifie une zone éloignée du trou
assert not np.isnan(data[0:40, 0:40]).any(), "NaN loin du trou"
def test_aspect_shared_preserves_nodata(self, synthetic_dem, tmp_path):
from lidar_pipeline.visualizations import SharedDEM, generate_aspect
dem_hole = self._dem_with_hole(synthetic_dem, tmp_path)
shared = SharedDEM(dem_hole, 5.0)
out = generate_aspect(dem_hole, "partage", tmp_path, 5.0, shared=shared)
assert out is not None and out.exists()
import rasterio
with rasterio.open(out) as src:
data = src.read(1)
assert np.isnan(data[80:120, 80:120]).all(), "le trou doit rester en nodata"
assert not np.isnan(data[0:40, 0:40]).any(), "NaN loin du trou"
def test_slope_and_hillshade_preserve_nodata(self, synthetic_dem, tmp_path):
from lidar_pipeline.visualizations import generate_slope, generate_hillshade
dem_hole = self._dem_with_hole(synthetic_dem, tmp_path)
import rasterio
for gen, name in ((generate_slope, "p"), (generate_hillshade, "h")):
out = gen(dem_hole, name, tmp_path, 5.0)
assert out is not None and out.exists()
with rasterio.open(out) as src:
data = src.read(1)
assert np.isnan(data[80:120, 80:120]).any(), f"{out.name} : trou disparu"

View File

@ -0,0 +1,81 @@
"""Tests du serveur web de génération de zones (webapp)."""
def test_bbox_to_cells_single_km_cell():
"""Une bbox couvrant ~1 km² retourne la cellule L93 correspondante."""
from lidar_pipeline.webapp import bbox_to_cells
# Cellule 1054,6882 : X∈[1054000,1055000], Y∈[6881000,6882000] (L93)
from rasterio.warp import transform as warp_transform
lons, lats = warp_transform('EPSG:2154', 'EPSG:4326',
[1054100, 1054900], [6881100, 6881900])
cells = bbox_to_cells(min(lons), min(lats), max(lons), max(lats))
assert (1054, 6882) in cells
# La sélection reste locale : pas de cellule lointaine
for (c, r) in cells:
assert abs(c - 1054) <= 1 and abs(r - 6882) <= 1
def test_bbox_to_cells_empty_for_tiny_bbox():
"""Une bbox quasi ponctuelle ne sélectionne rien (rétrécie sous 1 m)."""
from lidar_pipeline.webapp import bbox_to_cells
assert bbox_to_cells(7.850000, 48.930000, 7.850001, 48.930001) == []
def test_processed_cells(tmp_path):
"""processed_cells lit les dossiers de visualisations."""
from lidar_pipeline.webapp import processed_cells
vis = tmp_path / "visualisations" / "LHD_FXX_1054_6882_PTS_LAMB93_IGN69_r0p2"
vis.mkdir(parents=True)
(vis / "LHD_FXX_1054_6882_PTS_LAMB93_IGN69_r0p2_aspect.avif").write_bytes(b"x")
assert processed_cells(tmp_path) == {(1054, 6882)}
def test_missing_cells_filters_processed(tmp_path):
"""Les cellules déjà traitées sont exclues, les autres gardent leurs coins."""
from lidar_pipeline.webapp import missing_cells_with_corners
vis = tmp_path / "visualisations" / "LHD_FXX_1054_6882_PTS_LAMB93_IGN69"
vis.mkdir(parents=True)
(vis / "LHD_FXX_1054_6882_PTS_LAMB93_IGN69_aspect.avif").write_bytes(b"x")
todo = missing_cells_with_corners([(1054, 6882), (1055, 6882)], tmp_path)
assert len(todo) == 1
assert todo[0]['col'] == 1055 and todo[0]['row'] == 6882
assert len(todo[0]['corners']) == 4 # SW, SE, NE, NW
def test_missing_cells_include_done(tmp_path):
"""include_done=True conserve les cellules déjà traitées (régénération)."""
from lidar_pipeline.webapp import missing_cells_with_corners
vis = tmp_path / "visualisations" / "LHD_FXX_1054_6882_PTS_LAMB93_IGN69"
vis.mkdir(parents=True)
(vis / "LHD_FXX_1054_6882_PTS_LAMB93_IGN69_aspect.avif").write_bytes(b"x")
todo = missing_cells_with_corners([(1054, 6882), (1055, 6882)], tmp_path,
include_done=True)
assert {(t['col'], t['row']) for t in todo} == {(1054, 6882), (1055, 6882)}
def test_build_command_regenerate():
"""regenerate=True ajoute --force --force-classification à la commande."""
from lidar_pipeline.webapp import _build_command
cmd = " ".join(_build_command([(1054, 6882)], regenerate=True))
assert "--force" in cmd
assert "--force-classification" in cmd
cmd = " ".join(_build_command([(1054, 6882)]))
assert "--force" not in cmd
assert "--force-classification" not in cmd
def test_build_command_ground_classification():
"""La commande utilise la méthode de classification demandée (défaut : ign)."""
from lidar_pipeline.webapp import _build_command, GROUND_CLASS_METHODS
# Défaut : ign (pré-classification)
cmd = _build_command([(1054, 6882)])
i = cmd.index("--ground-classification")
assert cmd[i + 1] == "ign"
# Chaque méthode valide est transmise telle quelle, avec ou sans régénération
for method in GROUND_CLASS_METHODS:
for regenerate in (False, True):
cmd = _build_command([(1054, 6882)], regenerate=regenerate, ground_class=method)
i = cmd.index("--ground-classification")
assert cmd[i + 1] == method
assert ("--force" in cmd) == regenerate
assert ("--force-classification" in cmd) == regenerate

254
lidar_pipeline/webapp.py Normal file
View File

@ -0,0 +1,254 @@
"""Serveur web de la carte LiDAR : sert l'index et expose l'API de génération.
Lancé via `./run.sh --serve [PORT]` (input/ monté en écriture pour permettre
le téléchargement IGN). Endpoints :
GET / → carte interactive (output/index.html)
GET /api/status → état de la génération en cours (ou dernière terminée)
POST /api/preview → cellules 1 km intersectant une bbox WGS84 (option
regenerate=true pour inclure celles déjà générées)
POST /api/generate → télécharge (géoplateforme IGN) puis traite des cellules
(option regenerate=true ajoute --force --force-classification ;
option ground_class choisit la méthode de classification du sol)
Fichiers statiques : /assets (interface), /index_thumbs, /index_subtiles,
/visualisations, /DTM.
Un seul job à la fois : la génération lance `python -m lidar_pipeline` en
sous-processus avec --fetch-tiles + --file, journalisé dans .generation.log.
"""
import math
import os
import subprocess
import sys
import threading
import time
from pathlib import Path
from fastapi import FastAPI, HTTPException
from fastapi.responses import FileResponse, JSONResponse
from pydantic import BaseModel, Field
INPUT_DIR = Path(os.environ.get("LIDAR_INPUT_DIR", "/data/input"))
OUTPUT_DIR = Path(os.environ.get("LIDAR_OUTPUT_DIR", "/data/output"))
PORT = int(os.environ.get("LIDAR_PORT", "8973"))
LOG_FILE = OUTPUT_DIR / ".generation.log"
MAX_CELLS = 400 # garde-fou : ~400 km² max par demande
app = FastAPI(title="Carte LiDAR — génération de zones")
# assets/ (CSS/JS de l'interface) est créé dès le démarrage pour que le monteur
# statique soit actif même avant la première génération de l'index.
_assets_dir = OUTPUT_DIR / "assets"
_assets_dir.mkdir(parents=True, exist_ok=True)
for name in ("index_thumbs", "index_subtiles", "visualisations", "DTM"):
_dir = OUTPUT_DIR / name
if _dir.exists():
from fastapi.staticfiles import StaticFiles
app.mount(f"/{name}", StaticFiles(directory=str(_dir)), name=name)
@app.get("/assets/{file_path:path}")
def assets(file_path: str):
"""Sert les fichiers de l'interface sans cache (régénérés à chaque rebuild)."""
base = _assets_dir.resolve()
p = (_assets_dir / file_path).resolve()
if base not in p.parents or not p.is_file():
raise HTTPException(404, f"asset introuvable : {file_path}")
return FileResponse(str(p), headers={"Cache-Control": "no-cache, must-revalidate"})
# Méthodes de classification du sol acceptées (mêmes valeurs que --ground-classification).
GROUND_CLASS_METHODS = ("auto", "ign", "smrf", "csf")
class PreviewRequest(BaseModel):
bbox: list = Field(..., description="[ouest, sud, est, nord] en WGS84")
regenerate: bool = Field(False, description="Inclure les tuiles déjà générées")
class GenerateRequest(BaseModel):
tiles: list = Field(..., description="liste [col, row] (entiers km L93)")
regenerate: bool = Field(False, description="Régénérer les tuiles déjà générées")
ground_class: str = Field("ign",
description="Méthode de classification du sol : "
"auto, ign, smrf, csf")
ign_classes: str = Field("sol",
description="Classes LAS pour le MNT IGN : liste noms ou "
"codes séparés par virgules — sol(2), "
"unclassified(1), eau(9), virtuel(66), "
"pont(17), sursol(64). Mode pur, "
"aucune retouche. (défaut: sol)")
bare_earth: bool = Field(False,
description="Sol nu : DTM au retour le plus bas de "
"chaque cellule (requalifie le point le plus "
"bas en terrain)")
# --- État du job de génération -------------------------------------------
_job = {"proc": None, "started": None, "returncode": None, "cmd": None}
_job_lock = threading.Lock()
def bbox_to_cells(w, s, e, n):
"""Cellules L93 de 1 km (col, row) intersectant une bbox WGS84.
Une cellule (col, row) couvre X ∈ [col, col+1] km, Y ∈ [row-1, row] km.
"""
from rasterio.warp import transform as warp_transform
lons, lats = warp_transform("EPSG:4326", "EPSG:2154", [w, e, w, e], [s, s, n, n])
# warp_transform renvoie (xs, ys) dans la CRS cible
min_x, max_x = min(lons) + 0.5, max(lons) - 0.5 # rétrécit d'1 m : bords exclus
min_y, max_y = min(lats) + 0.5, max(lats) - 0.5
if max_x <= min_x or max_y <= min_y:
return []
cols = range(int(math.floor(min_x / 1000)), int(math.floor(max_x / 1000)) + 1)
rows = range(int(math.floor(min_y / 1000)) + 1, int(math.floor(max_y / 1000)) + 2)
return [(c, r) for r in rows for c in cols]
def processed_cells(output_dir):
"""Ensemble des cellules (col, row) ayant déjà des visualisations."""
from .index import scan_tiles
tiles = scan_tiles(Path(output_dir) / "visualisations")
return {(t["col"], t["row"]) for t in tiles}
def missing_cells_with_corners(cells, output_dir, include_done=False):
"""Filtre les cellules déjà traitées et calcule leurs coins WGS84.
Retourne [{col, row, corners: [[lat, lon] × 4 SW,SE,NE,NW}].
"""
from .index import attach_gps_bounds
done = processed_cells(output_dir)
todo = [{"col": c, "row": r} for (c, r) in cells
if include_done or (c, r) not in done]
if todo:
attach_gps_bounds(todo)
return todo
@app.get("/")
def root():
index = OUTPUT_DIR / "index.html"
if not index.exists():
return JSONResponse({"erreur": "index.html introuvable — lancez d'abord le pipeline"},
status_code=404)
# index.html est régénéré à chaque passe du pipeline : on interdit le cache
# navigateur pour ne pas servir une version périmée (ex. menu de génération).
return FileResponse(
str(index), media_type="text/html",
headers={
"Cache-Control": "no-cache, no-store, must-revalidate",
"Pragma": "no-cache",
"Expires": "0",
})
@app.get("/api/status")
def status():
with _job_lock:
proc = _job["proc"]
running = proc is not None and proc.poll() is None
return {
"running": running,
"started": _job["started"],
"returncode": _job["returncode"],
"cmd": _job["cmd"],
"log": _tail_log(40),
}
def _tail_log(n_lines):
try:
lines = LOG_FILE.read_text(encoding="utf-8", errors="replace").splitlines()
return lines[-n_lines:]
except Exception:
return []
@app.post("/api/preview")
def preview(req: PreviewRequest):
if len(req.bbox) != 4:
raise HTTPException(400, "bbox attendue : [ouest, sud, est, nord]")
w, s, e, n = (float(v) for v in req.bbox)
cells = bbox_to_cells(w, s, e, n)
capped = len(cells) > MAX_CELLS
todo = missing_cells_with_corners(cells[:MAX_CELLS], OUTPUT_DIR,
include_done=req.regenerate)
return {"count": len(todo), "capped": capped, "cells": todo}
def _build_command(tiles, regenerate=False, ground_class="ign", bare_earth=False, ign_classes="sol"):
"""Commande de génération : téléchargement IGN + traitement des fichiers.
La classification du sol est choisie via `ground_class` (défaut : "ign",
pré-classification IGN ; le pipeline bascule sur SMRF si un fichier ne la
contient pas). Avec ign_classes, on choisit les classes LAS extraites
pour le MNT (mode pur, ex. "sol,unclassified" pour combler les trous
sans retouche). Avec regenerate=True, force la reclassification et la
régénération des visualisations des tuiles déjà présentes. Avec
bare_earth=True, le DTM est ramené au retour le plus bas de chaque
cellule (sol nu).
"""
from .fetch_ign import tile_filename
# -u : sortie non bufferisée — le journal .generation.log doit être
# lu en temps réel par /api/status (progression affichée dans l'UI).
cmd = [sys.executable, "-u", "-m", "lidar_pipeline", str(INPUT_DIR),
"-o", str(OUTPUT_DIR), "-r", "0.5,0.2", "--only", "aspect",
"--ground-classification", ground_class,
"--ign-classes", ign_classes]
if bare_earth:
cmd += ["--bare-earth"]
if regenerate:
cmd += ["--force", "--force-classification"]
if os.environ.get("LIDAR_GPU", "") == "1":
cmd += ["-g", "all", "-w", os.environ.get("LIDAR_WORKERS", "2")]
cmd += ["--fetch-tiles"]
cmd += [f"{c},{r}" for (c, r) in tiles]
cmd += ["--file"]
cmd += [tile_filename(c, r) for (c, r) in tiles]
return cmd
@app.post("/api/generate")
def generate(req: GenerateRequest):
tiles = []
for pair in req.tiles:
if not (isinstance(pair, list) and len(pair) == 2):
raise HTTPException(400, f"tuile invalide : {pair!r} (attendu [col, row])")
tiles.append((int(pair[0]), int(pair[1])))
if not tiles:
raise HTTPException(400, "aucune tuile fournie")
if len(tiles) > MAX_CELLS:
raise HTTPException(400, f"trop de tuiles ({len(tiles)}) — max {MAX_CELLS}")
if req.ground_class not in GROUND_CLASS_METHODS:
raise HTTPException(
400, f"méthode de classification invalide : {req.ground_class!r} "
f"(attendu : {', '.join(GROUND_CLASS_METHODS)})")
with _job_lock:
proc = _job["proc"]
if proc is not None and proc.poll() is None:
raise HTTPException(409, "une génération est déjà en cours")
cmd = _build_command(tiles, regenerate=req.regenerate, ground_class=req.ground_class, bare_earth=req.bare_earth, ign_classes=req.ign_classes)
LOG_FILE.parent.mkdir(parents=True, exist_ok=True)
log_fh = open(LOG_FILE, "w", encoding="utf-8")
_job.update({"proc": None, "started": time.time(), "returncode": None, "cmd": cmd})
p = subprocess.Popen(cmd, stdout=log_fh, stderr=subprocess.STDOUT,
cwd="/app" if Path("/app").exists() else None)
def _watch():
rc = p.wait()
log_fh.close()
with _job_lock:
_job["returncode"] = rc
threading.Thread(target=_watch, daemon=True).start()
_job["proc"] = p
return {"demarré": True, "tuiles": len(tiles), "commande": " ".join(cmd)}
if __name__ == "__main__":
import uvicorn
uvicorn.run(app, host="0.0.0.0", port=PORT)

79
run.sh
View File

@ -11,7 +11,9 @@
# -f / --force Régénérer tous les fichiers # -f / --force Régénérer tous les fichiers
# --keep-tif Conserver les fichiers TIFF # --keep-tif Conserver les fichiers TIFF
# --force-classification # --force-classification
# --ground-classification {auto,smrf,csf} # --ground-classification {auto,ign,smrf,csf}
# --ign-classes CLASSES
# Classes LAS pour le MNT IGN (ex: sol,unclassified ; défaut: sol)
# --quality N Qualité image 1-100 (défaut: 98) # --quality N Qualité image 1-100 (défaut: 98)
# --lossless Compression lossless # --lossless Compression lossless
# --format FMT Format de sortie : avif (défaut) ou webp # --format FMT Format de sortie : avif (défaut) ou webp
@ -44,8 +46,11 @@ if [ $# -eq 0 ]; then
echo " --force-classification" echo " --force-classification"
echo " Reclassifier le sol même si le fichier .las existe" echo " Reclassifier le sol même si le fichier .las existe"
echo " --keep-tif Conserver les TIFF pour régénérer les AVIF" echo " --keep-tif Conserver les TIFF pour régénérer les AVIF"
echo " --ground-classification {auto,smrf,csf}" echo " --ground-classification {auto,ign,smrf,csf}"
echo " Méthode de classification du sol (défaut: auto)" echo " Méthode de classification du sol (défaut: auto = pré-classification IGN)"
echo " --ign-classes CLASSES"
echo " Classes LAS pour le MNT IGN : sol,unclassified,eau,virtuel,pont,sursol"
echo " ou codes (ex: 2,1). Mode pur, aucune retouche. (défaut: sol)"
echo " --quality N Qualité image 1-100 (défaut: 98, 100=lossless)" echo " --quality N Qualité image 1-100 (défaut: 98, 100=lossless)"
echo " --lossless Compression lossless (équivalent à --quality 100)" echo " --lossless Compression lossless (équivalent à --quality 100)"
echo " --format FMT Format de sortie : avif (défaut) ou webp" echo " --format FMT Format de sortie : avif (défaut) ou webp"
@ -83,6 +88,7 @@ VERBOSE_FLAG=""
FORCE_FLAG="" FORCE_FLAG=""
FILE_ARGS="" FILE_ARGS=""
GROUND_METHOD="" GROUND_METHOD=""
IGN_CLASSES_FLAG=""
FORCE_CLASSIFY_FLAG="" FORCE_CLASSIFY_FLAG=""
KEEP_TIF_FLAG="" KEEP_TIF_FLAG=""
QUALITY="" QUALITY=""
@ -92,6 +98,9 @@ SKIP_FLAG=""
TEST_FLAG=0 TEST_FLAG=0
REBUILD_INDEX_FLAG="" REBUILD_INDEX_FLAG=""
NO_INDEX_FLAG="" NO_INDEX_FLAG=""
FETCH_TILES_ARGS=""
SERVE_FLAG=0
SERVE_PORT=8973
# Parse arguments manually (more robust than getopts for mixed short/long options) # Parse arguments manually (more robust than getopts for mixed short/long options)
while [ $# -gt 0 ]; do while [ $# -gt 0 ]; do
@ -117,6 +126,8 @@ while [ $# -gt 0 ]; do
--keep-tif) KEEP_TIF_FLAG="--keep-tif"; shift ;; --keep-tif) KEEP_TIF_FLAG="--keep-tif"; shift ;;
--ground-classification) GROUND_METHOD="$2"; shift 2 ;; --ground-classification) GROUND_METHOD="$2"; shift 2 ;;
--ground-classification=*) GROUND_METHOD="${1#--ground-classification=}"; shift ;; --ground-classification=*) GROUND_METHOD="${1#--ground-classification=}"; shift ;;
--ign-classes) IGN_CLASSES_FLAG="--ign-classes $2"; shift 2 ;;
--ign-classes=*) IGN_CLASSES_FLAG="--ign-classes=${1#--ign-classes=}"; shift ;;
--quality) QUALITY="--quality $2"; shift 2 ;; --quality) QUALITY="--quality $2"; shift 2 ;;
--lossless) QUALITY="--lossless"; shift ;; --lossless) QUALITY="--lossless"; shift ;;
--format) FORMAT_FLAG="--format $2"; shift 2 ;; --format) FORMAT_FLAG="--format $2"; shift 2 ;;
@ -125,7 +136,12 @@ while [ $# -gt 0 ]; do
--file) shift; while [ $# -gt 0 ] && [[ ! "$1" =~ ^- ]]; do FILE_ARGS="$FILE_ARGS $1"; shift; done ;; --file) shift; while [ $# -gt 0 ] && [[ ! "$1" =~ ^- ]]; do FILE_ARGS="$FILE_ARGS $1"; shift; done ;;
--rebuild-index) REBUILD_INDEX_FLAG="--rebuild-index"; shift ;; --rebuild-index) REBUILD_INDEX_FLAG="--rebuild-index"; shift ;;
--no-index) NO_INDEX_FLAG="--no-index"; shift ;; --no-index) NO_INDEX_FLAG="--no-index"; shift ;;
--test) TEST_FLAG=1 ;; --fetch-tiles) shift; while [ $# -gt 0 ] && [[ ! "$1" =~ ^- ]]; do FETCH_TILES_ARGS="$FETCH_TILES_ARGS $1"; shift; done ;;
--test) TEST_FLAG=1; shift ;;
--serve)
SERVE_FLAG=1; shift
if [ $# -gt 0 ] && [[ ! "$1" =~ ^- ]]; then SERVE_PORT="$1"; shift; fi
;;
-h|--help|-help) -h|--help|-help)
echo "Pipeline LiDAR Archéologique" echo "Pipeline LiDAR Archéologique"
echo "" echo ""
@ -141,14 +157,20 @@ while [ $# -gt 0 ]; do
echo " --force-classification" echo " --force-classification"
echo " Reclassifier le sol même si le fichier .las existe" echo " Reclassifier le sol même si le fichier .las existe"
echo " --keep-tif Conserver les TIFF pour régénérer les AVIF" echo " --keep-tif Conserver les TIFF pour régénérer les AVIF"
echo " --ground-classification {auto,smrf,csf}" echo " --ground-classification {auto,ign,smrf,csf}"
echo " Méthode de classification du sol (défaut: auto)" echo " Méthode de classification du sol (défaut: auto = pré-classification IGN)"
echo " --ign-classes CLASSES"
echo " Classes LAS pour le MNT IGN : sol,unclassified,eau,virtuel,pont,sursol"
echo " ou codes (ex: 2,1). Mode pur, aucune retouche. (défaut: sol)"
echo " --quality N Qualité image 1-100 (défaut: 98, 100=lossless)" echo " --quality N Qualité image 1-100 (défaut: 98, 100=lossless)"
echo " --lossless Compression lossless (équivalent à --quality 100)" echo " --lossless Compression lossless (équivalent à --quality 100)"
echo " --format FMT Format de sortie : avif (défaut) ou webp" echo " --format FMT Format de sortie : avif (défaut) ou webp"
echo " --only VIZ... Générer uniquement ces visualisations" echo " --only VIZ... Générer uniquement ces visualisations"
echo " --skip VIZ... Exclure ces visualisations" echo " --skip VIZ... Exclure ces visualisations"
echo " --file NOM... Traiter un ou plusieurs fichiers LAZ" echo " --file NOM... Traiter un ou plusieurs fichiers LAZ"
echo " --fetch-tiles COL,ROW..."
echo " Télécharger ces dalles depuis l'IGN puis les traiter"
echo " --serve [PORT] Servir la carte + API de génération de zones (défaut: 8973)"
echo " --test Exécuter les tests unitaires" echo " --test Exécuter les tests unitaires"
echo " -h Afficher cette aide" echo " -h Afficher cette aide"
echo "" echo ""
@ -181,6 +203,13 @@ done
# Trim FILE_ARGS whitespace # Trim FILE_ARGS whitespace
FILE_ARGS=$(echo "$FILE_ARGS" | xargs) FILE_ARGS=$(echo "$FILE_ARGS" | xargs)
FETCH_TILES_ARGS=$(echo "$FETCH_TILES_ARGS" | xargs)
# input/ en lecture seule sauf si --fetch-tiles (téléchargement IGN dans input/)
INPUT_RO=":ro"
if [ -n "$FETCH_TILES_ARGS" ]; then
INPUT_RO=""
fi
# Check for --test flag first # Check for --test flag first
if [ "$TEST_FLAG" -eq 1 ]; then if [ "$TEST_FLAG" -eq 1 ]; then
@ -207,6 +236,31 @@ fi
# Créer les répertoires s'ils n'existent pas # Créer les répertoires s'ils n'existent pas
mkdir -p "$INPUT_DIR" "$OUTPUT_DIR" mkdir -p "$INPUT_DIR" "$OUTPUT_DIR"
# Mode --serve : carte + API de génération (input/ en écriture)
if [ "$SERVE_FLAG" -eq 1 ]; then
echo "============================================"
echo " Carte LiDAR + API de génération"
echo "============================================"
echo " http://127.0.0.1:${SERVE_PORT}/"
echo " GPU : $([ -n "$GPU_FLAG" ] && echo 'OUI' || echo 'non')"
echo " Bouton « + Zone » : dessiner un rectangle → téléchargement IGN + rendu"
echo "============================================"
GPU_ENV_FLAG=""
[ -n "$GPU_FLAG" ] && GPU_ENV_FLAG="-e LIDAR_GPU=1"
if [ -n "$GPU_ARG" ] && [ "$GPU_ARG" != "all" ]; then
GPU_ENV_FLAG="-e LIDAR_GPU=1 -e CUDA_VISIBLE_DEVICES=${GPU_ARG} --gpus '""device=${GPU_ARG}""'"
fi
exec docker run --rm --init $GPU_FLAG $CUDA_ENV_FLAG \
--user 1000:1000 \
-p "${SERVE_PORT}:8973" \
-v "${INPUT_DIR}:/data/input" \
-v "${OUTPUT_DIR}:/data/output" \
-e LIDAR_INPUT_DIR=/data/input -e LIDAR_OUTPUT_DIR=/data/output \
$GPU_ENV_FLAG \
"$IMAGE_NAME" \
python3 -m uvicorn lidar_pipeline.webapp:app --host 0.0.0.0 --port 8973
fi
# Lancer le pipeline # Lancer le pipeline
echo "============================================" echo "============================================"
echo " Pipeline LiDAR Archéologique" echo " Pipeline LiDAR Archéologique"
@ -220,7 +274,7 @@ echo " Force classif.: $([ -n "$FORCE_CLASSIFY_FLAG" ] && echo 'OUI' || echo 'n
echo " Keep TIFF : $([ -n "$KEEP_TIF_FLAG" ] && echo 'OUI' || echo 'non')" echo " Keep TIFF : $([ -n "$KEEP_TIF_FLAG" ] && echo 'OUI' || echo 'non')"
echo " Qualité image : $([ -n "$QUALITY" ] && echo "$QUALITY" || echo '98')" echo " Qualité image : $([ -n "$QUALITY" ] && echo "$QUALITY" || echo '98')"
echo " Format : $([ -n "$FORMAT_FLAG" ] && echo "${FORMAT_FLAG#--format }" || echo 'avif')" echo " Format : $([ -n "$FORMAT_FLAG" ] && echo "${FORMAT_FLAG#--format }" || echo 'avif')"
echo " Classification sol : $([ -n "$GROUND_METHOD" ] && echo "$GROUND_METHOD" || echo 'auto')" echo " Classification sol : $([ -n "$GROUND_METHOD" ] && echo "$GROUND_METHOD" || echo 'auto')$([ -n "$IGN_CLASSES_FLAG" ] && echo " ${IGN_CLASSES_FLAG#--ign-classes }")"
if [ -n "$ONLY_FLAG" ]; then if [ -n "$ONLY_FLAG" ]; then
echo " Visualisations: uniquement${ONLY_FLAG#--only}" echo " Visualisations: uniquement${ONLY_FLAG#--only}"
elif [ -n "$SKIP_FLAG" ]; then elif [ -n "$SKIP_FLAG" ]; then
@ -229,6 +283,9 @@ fi
if [ -n "$FILE_ARGS" ]; then if [ -n "$FILE_ARGS" ]; then
echo " Fichiers :${FILE_ARGS}" echo " Fichiers :${FILE_ARGS}"
fi fi
if [ -n "$FETCH_TILES_ARGS" ]; then
echo " Dalles IGN :${FETCH_TILES_ARGS} (input/ monté en écriture)"
fi
echo "============================================" echo "============================================"
CMD_ARGS="-o /data/output -r $RESOLUTION -w $WORKERS $VERBOSE_FLAG $FORCE_FLAG $FORCE_CLASSIFY_FLAG $KEEP_TIF_FLAG $QUALITY $FORMAT_FLAG" CMD_ARGS="-o /data/output -r $RESOLUTION -w $WORKERS $VERBOSE_FLAG $FORCE_FLAG $FORCE_CLASSIFY_FLAG $KEEP_TIF_FLAG $QUALITY $FORMAT_FLAG"
@ -241,6 +298,9 @@ fi
if [ -n "$GROUND_METHOD" ]; then if [ -n "$GROUND_METHOD" ]; then
CMD_ARGS="$CMD_ARGS --ground-classification $GROUND_METHOD" CMD_ARGS="$CMD_ARGS --ground-classification $GROUND_METHOD"
fi fi
if [ -n "$IGN_CLASSES_FLAG" ]; then
CMD_ARGS="$CMD_ARGS $IGN_CLASSES_FLAG"
fi
if [ -n "$ONLY_FLAG" ]; then if [ -n "$ONLY_FLAG" ]; then
CMD_ARGS="$CMD_ARGS $ONLY_FLAG" CMD_ARGS="$CMD_ARGS $ONLY_FLAG"
fi fi
@ -256,6 +316,9 @@ fi
if [ -n "$NO_INDEX_FLAG" ]; then if [ -n "$NO_INDEX_FLAG" ]; then
CMD_ARGS="$CMD_ARGS $NO_INDEX_FLAG" CMD_ARGS="$CMD_ARGS $NO_INDEX_FLAG"
fi fi
if [ -n "$FETCH_TILES_ARGS" ]; then
CMD_ARGS="$CMD_ARGS --fetch-tiles $FETCH_TILES_ARGS"
fi
# Build CUDA_VISIBLE_DEVICES env var from GPU_ARG # Build CUDA_VISIBLE_DEVICES env var from GPU_ARG
CUDA_ENV_FLAG="" CUDA_ENV_FLAG=""
@ -271,7 +334,7 @@ fi
docker run --rm --init $GPU_FLAG $CUDA_ENV_FLAG \ docker run --rm --init $GPU_FLAG $CUDA_ENV_FLAG \
--user 1000:1000 \ --user 1000:1000 \
-v "${INPUT_DIR}:/data/input:ro" \ -v "${INPUT_DIR}:/data/input${INPUT_RO}" \
-v "${OUTPUT_DIR}:/data/output" \ -v "${OUTPUT_DIR}:/data/output" \
"$IMAGE_NAME" \ "$IMAGE_NAME" \
python3 -m lidar_pipeline /data/input \ python3 -m lidar_pipeline /data/input \