diff --git a/AGENTS.md b/AGENTS.md
index 7537301..3ca4bf3 100644
--- a/AGENTS.md
+++ b/AGENTS.md
@@ -8,6 +8,9 @@
- lint: not configured
- format: not configured
- 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`
## Conventions
@@ -21,6 +24,7 @@
- **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.
- **`_`-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
diff --git a/docker-compose.yml b/docker-compose.yml
index 65cc317..2832084 100644
--- a/docker-compose.yml
+++ b/docker-compose.yml
@@ -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:
- lidar:
+ # Carte interactive + API de génération (mode ./run.sh --serve)
+ serve:
build: .
- container_name: lidar-archeo
+ image: lidar-lidar
+ container_name: lidar-serve
+ init: true
user: "1000:1000"
+ gpus: all
+ ports:
+ - "8973:8973"
volumes:
- # Mount your LAZ files directory here
- - ./input:/data/input:ro
- # Output directory
+ # input/ en écriture : l'API y télécharge les dalles IGN manquantes
+ - ./input:/data/input
- ./output:/data/output
- # Optional: Mount a large data directory
- # - /path/to/your/laz/files:/data/input:ro
environment:
- TZ=Europe/Paris
- # Processing parameters
- - RESOLUTION=0.5
- - WHITEBOX_THREADS=4
- # Resource limits (adjust based on your system)
- deploy:
- resources:
- limits:
- cpus: '4'
- memory: 8G
- reservations:
- cpus: '2'
- memory: 4G
- # Override default command
- command: ["process_lidar.py", "/data/input", "-o", "/data/output", "-r", "0.5"]
+ - LIDAR_INPUT_DIR=/data/input
+ - LIDAR_OUTPUT_DIR=/data/output
+ # Les générations lancées depuis la carte utilisent le GPU
+ - LIDAR_GPU=1
+ - LIDAR_WORKERS=2
+ command: python3 -m uvicorn lidar_pipeline.webapp:app --host 0.0.0.0 --port 8973
+ restart: unless-stopped
- # 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:
build: .
+ image: lidar-lidar
container_name: lidar-jupyter
+ init: true
ports:
- "8888:8888"
volumes:
diff --git a/docs/GROUND_CLASSIFICATION.md b/docs/GROUND_CLASSIFICATION.md
new file mode 100644
index 0000000..aafd4fe
--- /dev/null
+++ b/docs/GROUND_CLASSIFICATION.md
@@ -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
diff --git a/lidar_pipeline/cli.py b/lidar_pipeline/cli.py
index e79206f..39b0f96 100644
--- a/lidar_pipeline/cli.py
+++ b/lidar_pipeline/cli.py
@@ -97,7 +97,10 @@ def main():
)
parser.add_argument(
"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(
"-o", "--output",
@@ -134,7 +137,16 @@ def main():
parser.add_argument(
"--force-classification",
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(
"--keep-tif",
@@ -143,9 +155,22 @@ def main():
)
parser.add_argument(
"--ground-classification",
- choices=["auto", "smrf", "csf"],
+ choices=["auto", "ign", "smrf", "csf"],
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(
"--quality",
@@ -185,6 +210,14 @@ def main():
default=None,
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(
"-v", "--verbose",
action="store_true",
@@ -253,6 +286,21 @@ def main():
logger.warning("Aucune tuile traitée trouvée — carte globale non générée")
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
# Parse --only and --skip: accept comma-separated values
only_viz = None
@@ -268,8 +316,10 @@ def main():
workers=args.workers,
force=args.force,
ground_method=args.ground_classification,
+ ign_classes=args.ign_classes,
force_classify=args.force_classification,
keep_tif=args.keep_tif,
+ bare_earth=args.bare_earth,
quality=quality,
only_viz=only_viz,
skip_viz=skip_viz,
@@ -315,30 +365,8 @@ def main():
logger.info(f"Traitement de {len(unique_files)} fichier(s) sélectionné(s)")
for laz_file in unique_files:
logger.info(f" → {laz_file.name}")
- for laz_file in unique_files:
- pipeline.process_file(laz_file)
-
- # 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}")
+ # Réutilise process_all : workers parallèles, résumé, index, nettoyage
+ pipeline.process_all(files=unique_files)
else:
pipeline.process_all()
except Exception as e:
diff --git a/lidar_pipeline/dtm.py b/lidar_pipeline/dtm.py
index 5eb8c78..e0f4a74 100644
--- a/lidar_pipeline/dtm.py
+++ b/lidar_pipeline/dtm.py
@@ -1,7 +1,10 @@
"""DTM generation from classified LiDAR point clouds.
-Handles ground classification via PDAL (SMRF or CSF) and DTM rasterisation
-using scipy binned_statistic_2d. Zones without LiDAR data remain as NaN.
+Handles ground classification via PDAL (IGN supplier pre-classification,
+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
@@ -16,8 +19,67 @@ from scipy.stats import binned_statistic_2d
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.
All methods include a ReturnNumber/NumberOfReturns >= 1 filter to handle
@@ -33,7 +95,10 @@ def _create_ground_pipeline(input_laz, output_las, method):
Args:
input_laz: Path to input LAZ/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:
JSON string of the PDAL pipeline.
@@ -44,6 +109,38 @@ def _create_ground_pipeline(input_laz, output_las, method):
"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 = {
"type": "filters.assign",
@@ -68,12 +165,6 @@ def _create_ground_pipeline(input_laz, output_las, method):
"multiplier": 3.0
}
- # Classification filter (ground points only)
- ground_filter = {
- "type": "filters.range",
- "limits": "Classification[2:2]"
- }
-
# Method-specific ground classification filter
if method == 'smrf':
ground_step = {
@@ -85,9 +176,12 @@ def _create_ground_pipeline(input_laz, output_las, method):
"scalar": 1.25
}
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 = {
"type": "filters.csf",
- "resolution": 0.5,
+ "resolution": 1.0,
"rigidness": 3,
"smooth": True,
"threshold": 0.5
@@ -119,6 +213,11 @@ def create_smrf_pipeline(input_laz, output_las):
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):
"""Create a PDAL pipeline JSON for CSF ground classification."""
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.
Returns:
- String: 'smrf' or 'csf'
+ String: 'ign', 'smrf' or 'csf'
"""
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")
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)
# Height variance (always available)
@@ -317,14 +432,17 @@ def detect_ground_method(laz_file):
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.
Args:
laz_file: Path to input LAZ/LAS file.
temp_dir: Directory for temporary files (pipeline.json, ground.las).
- method: Ground classification method ('auto', 'smrf', or 'csf').
+ method: Ground classification method ('auto', 'ign', 'smrf' or 'csf').
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:
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:
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
from .pipeline import _file_basename
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:
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}")
output_las.unlink()
- pipeline_json = _create_ground_pipeline(laz_file, output_las, method)
- pipeline_file = temp_dir / f"pipeline_{method}.json"
+ pipeline_json = _create_ground_pipeline(laz_file, output_las, method, ign_codes=ign_codes)
+ pipeline_file = temp_dir / f"pipeline_{method_label}.json"
with open(pipeline_file, 'w') as f:
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:
logger.error(f" ✗ Fichier ground vide (taille < 100 octets)")
output_las.unlink(missing_ok=True)
- # Fallback: if CSF produced no ground points, retry with SMRF
- if method == 'csf':
- return _fallback_to_smrf(laz_file, temp_dir, laz_base, force)
+ # Fallback: si la méthode ne produit aucun point sol, réessayer avec SMRF
+ if method in ('csf', 'ign'):
+ return _fallback_to_smrf(laz_file, temp_dir, laz_base, force, source=method_label)
return None
logger.info(f" ✓ Classification sol {method.upper()} terminée")
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)
logger.warning(f" ✗ Erreur classification PDAL ({method.upper()}): {error_msg}")
- # Fallback: if CSF failed, retry with SMRF
- if method == 'csf':
- return _fallback_to_smrf(laz_file, temp_dir, laz_base, force)
+ # Fallback: si CSF ou la pré-classification échouent, réessayer avec SMRF
+ if method in ('csf', 'ign'):
+ 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
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
-def _fallback_to_smrf(laz_file, temp_dir, laz_base, force=False):
- """Retry ground classification with SMRF when CSF fails.
+def _fallback_to_smrf(laz_file, temp_dir, laz_base, force=False, source='csf'):
+ """Retry ground classification with SMRF when CSF/IGN fails.
CSF (Cloth Simulation Filter) can fail on certain terrain types where
- SMRF (Simple Morphological Filter) succeeds. This fallback ensures
- processing continues even when auto-detection selects CSF incorrectly.
+ SMRF (Simple Morphological Filter) succeeds, and a file without usable
+ pre-classification produces an empty ground extract. This fallback ensures
+ processing continues even when the selected method fails.
Args:
laz_file: Path to input LAZ/LAS file.
temp_dir: Directory for temporary files.
laz_base: Base name for the file.
force: If True, reclassify even if output exists.
+ source: Method that failed ('csf' or 'ign').
Returns:
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
- csf_output = temp_dir / f"{laz_base}_ground_csf.las"
- if csf_output.exists():
- csf_output.unlink(missing_ok=True)
+ # Clean up failed output if it exists
+ failed_output = temp_dir / f"{laz_base}_ground_{source}.las"
+ if failed_output.exists():
+ failed_output.unlink(missing_ok=True)
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
-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.
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.
force: If True, regenerate even if DTM already exists.
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:
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 = 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))
if nan_count > 0:
total = dtm.size
@@ -558,8 +816,6 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, output
if filled_count > 0:
dtm = np.where(small_gap_mask, dtm_filled, dtm)
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
output_tif = dtm_dir / f"{basename}_dtm{output_suffix}.tif"
diff --git a/lidar_pipeline/fetch_ign.py b/lidar_pipeline/fetch_ign.py
new file mode 100644
index 0000000..a6d3443
--- /dev/null
+++ b/lidar_pipeline/fetch_ign.py
@@ -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
diff --git a/lidar_pipeline/index.py b/lidar_pipeline/index.py
index eb14b96..e29f5c8 100644
--- a/lidar_pipeline/index.py
+++ b/lidar_pipeline/index.py
@@ -1,13 +1,19 @@
"""Carte continue interactive des tuiles LiDAR traitées.
-Génère une page HTML unique (output/index.html) présentant toutes les tuiles
-1×1 km traitées, affichées en carte continue sans bordure, avec zoom/pan fluide
-et sélecteur de visualisation. Au clic sur une tuile, un modal affiche l'image
-avec légende et permet de télécharger en PDF.
+Génère l'interface web de consultation (servie par webapp.py via --serve) :
+ - output/index.html : coquille HTML minimaliste embarquant les données
+ - output/assets/app.css : styles (panneaux flottants type SIG)
+ - output/assets/app.js : logique carte (Leaflet), couches, infos tuiles
-Sortie:
- - output/index.html : page interactive auto-suffisante
- - output/index_thumbs/*.jpg : vignettes JPEG (~256px) par tuile/visualisation
+Fonctionnalités :
+ - Carte continue zoom/pan (Leaflet), tuiles rotées selon la projection
+ Lambert 93, vignettes ↔ images pleine résolution selon le zoom.
+ - Panneau de couches superposables avec opacité individuelle et
+ réordonnancement par glisser-déposer (persisté en localStorage).
+ - Panneau d'infos par tuile : résolution, méthode de classification du sol
+ (lue depuis output/DTM/*_dtm_method.txt), dates et tailles des fichiers.
+ - Génération de nouvelles zones via l'API web (./run.sh --serve) :
+ dessin d'un rectangle → téléchargement IGN + traitement.
Intégration:
- Appelé automatiquement à la fin de process_all() dans pipeline.py
@@ -17,12 +23,14 @@ Intégration:
import json
import logging
import re
+import time
+from datetime import datetime
from pathlib import Path
logger = logging.getLogger("lidar")
-# Noms d'affichage (français) pour le sélecteur de visualisation.
+# Noms d'affichage (français) pour le panneau de couches.
# Clé = mot-clé dans le nom de fichier de sortie (post-basename).
VIZ_LABELS = {
'hillshade_multi': 'Hillshade multidirectionnel',
@@ -43,31 +51,22 @@ VIZ_LABELS = {
'topo': 'Carte topographique IGN',
}
-# Informations de colormap pour génération de légende côté navigateur.
-# Utilisées par le JS pour dessiner la légende sur le canvas avant export PDF.
-VIZ_COLORMAPS = {
- 'hillshade_multi': {'cmap': 'gray', 'diverging': False},
- 'slope': {'cmap': 'inferno', 'diverging': False},
- 'aspect': {'cmap': 'twilight', 'diverging': False},
- 'mslrm': {'cmap': 'seismic', 'diverging': True},
- 'sailore': {'cmap': 'seismic', 'diverging': True},
- 'positive_openness': {'cmap': 'YlOrBr', 'diverging': False},
- 'negative_openness': {'cmap': 'PuBu', 'diverging': False},
- 'svf': {'cmap': 'hot_r', 'diverging': False},
- 'aniso_open': {'cmap': 'seismic', 'diverging': True},
- 'roughness': {'cmap': 'plasma', 'diverging': False},
- 'wavelet': {'cmap': 'cividis', 'diverging': False},
- 'flow_acc': {'cmap': 'YlGn', 'diverging': False},
- 'solar': {'cmap': 'gray', 'diverging': False},
- 'anomaly': {'cmap': 'YlOrRd', 'diverging': False},
- 'ortho': {'cmap': 'rgb', 'diverging': False},
- 'topo': {'cmap': 'rgb', 'diverging': False},
-}
-
-# Visualisation par défaut pour la vignette (si disponible).
+# Couche activée par défaut à l'ouverture de la carte.
DEFAULT_VIZ = 'hillshade_multi'
-# Ordre préféré pour le choix de la vignette de repli.
+# Visualisations découpées en sous-tuiles (quadrants 500 m) pour alléger la
+# carte — les autres viz restent disponibles en repli dalle entière.
+# Liste vide = découper toutes les visualisations disponibles.
+_CARTO_SUBTILED_VIZ = ('aspect', 'hillshade_multi')
+
+# Encodage AVIF des sous-tuiles (benchmark sur dalles d'aspect réelles :
+# q55 ≈ −45 % vs WebP q82, visuellement propre sur rampes de couleur).
+# speed 0 (lent/meilleur) → 10 (rapide) ; 6+8 quasi identique en taille,
+# on garde 8 (2× plus rapide) pour 2500×2500 px.
+_SUBTILE_AVIF_QUALITY = 55
+_SUBTILE_AVIF_SPEED = 8
+
+# Ordre préféré des couches (bas → haut de pile) et choix de la vignette de repli.
_VIZ_FALLBACK_ORDER = [
'hillshade_multi', 'svf', 'slope', 'mslrm', 'positive_openness',
'negative_openness', 'aspect', 'sailore', 'aniso_open', 'roughness',
@@ -114,6 +113,13 @@ def _strip_res_suffix(dirname):
return dirname, 0.5
+def _res_suffix_str(resolution):
+ """Suffixe de nommage d'une résolution (miroir de pipeline._res_suffix)."""
+ if resolution == 0.5:
+ return ''
+ return '_r' + f'{resolution}'.replace('.', 'p')
+
+
def scan_tiles(vis_dir):
"""Scanne le dossier des visualisations pour inventorier les tuiles traitées.
@@ -198,6 +204,154 @@ def compute_bbox(tiles):
}
+def compute_zones(tiles, proximity_threshold=15):
+ """Regroupe les tuiles en zones géographiques par clustering de proximité.
+
+ Les tuiles à moins de proximity_threshold km les unes des autres
+ sont regroupées dans la même zone. Les zones sont triées par taille
+ décroissante (plus grande zone en premier).
+
+ Args:
+ tiles: Liste de dictionnaires de tuiles.
+ proximity_threshold: Distance maximale en km pour regrouper deux tuiles.
+
+ Returns:
+ Liste de dictionnaires de zones:
+ {label, tiles, bbox}
+ """
+ if not tiles:
+ return []
+
+ # Union-Find pour le clustering
+ parent = list(range(len(tiles)))
+
+ def find(x):
+ while parent[x] != x:
+ parent[x] = parent[parent[x]]
+ x = parent[x]
+ return x
+
+ def union(x, y):
+ px, py = find(x), find(y)
+ if px != py:
+ parent[px] = py
+
+ # Regrouper les tuiles proches
+ for i in range(len(tiles)):
+ for j in range(i + 1, len(tiles)):
+ dc = abs(tiles[i]['col'] - tiles[j]['col'])
+ dr = abs(tiles[i]['row'] - tiles[j]['row'])
+ if max(dc, dr) <= proximity_threshold:
+ union(i, j)
+
+ # Construire les zones
+ zone_members = {}
+ for i in range(len(tiles)):
+ root = find(i)
+ if root not in zone_members:
+ zone_members[root] = []
+ zone_members[root].append(tiles[i])
+
+ zones = []
+ for i, zone_tiles in enumerate(zone_members.values(), 1):
+ zone_bbox = compute_bbox(zone_tiles)
+ zones.append({
+ 'label': f'Zone {i} ({len(zone_tiles)} tuiles)',
+ 'tiles': zone_tiles,
+ 'bbox': zone_bbox,
+ })
+
+ # Trier par taille décroissante
+ zones.sort(key=lambda z: len(z['tiles']), reverse=True)
+ return zones
+
+
+def _approx_l93_to_wgs84(x_m, y_m):
+ """Approximation affine Lambert 93 → WGS84 (fallback sans rasterio).
+
+ Origine exacte: (700000, 6600000) L93 ↔ (3.0°E, 46.5°N).
+ Précision de l'ordre du km — utilisée seulement si rasterio/PROJ
+ est indisponible (jamais le cas dans le Docker).
+ """
+ import math
+ lat = 46.5 + (y_m - 6600000.0) / 111320.0
+ lon = 3.0 + (x_m - 700000.0) / (111320.0 * math.cos(math.radians(47.0)))
+ return lon, lat
+
+
+def attach_gps_bounds(tiles):
+ """Attache à chaque tuile ses coins GPS pour l'affichage Leaflet.
+
+ Chaque tuile 1×1 km est définie par son coin nord-ouest en km L93
+ (col, row) → X ∈ [col, col+1] km, Y ∈ [row-1, row] km.
+ (Vérifié sur les bounds des DTM : X_min = col×1000, Y_max = row×1000.)
+
+ Utilise rasterio.warp (conversion PROJ exacte) si disponible,
+ sinon l'approximation affine _approx_l93_to_wgs84.
+
+ Ajoute à chaque tuile:
+ corners: [[lat, lon] × 4] dans l'ordre SW, SE, NE, NW
+ bounds : [[lat_sud, lon_ouest], [lat_nord, lon_est]]
+ """
+ try:
+ from rasterio.warp import transform as warp_transform
+ xs = []
+ ys = []
+ for t in tiles:
+ # SW, SE, NE, NW — Y du bord sud = (row-1)×1000, bord nord = row×1000
+ xs.extend([t['col'] * 1000, (t['col'] + 1) * 1000,
+ (t['col'] + 1) * 1000, t['col'] * 1000])
+ ys.extend([(t['row'] - 1) * 1000, (t['row'] - 1) * 1000,
+ t['row'] * 1000, t['row'] * 1000])
+ lons, lats = warp_transform('EPSG:2154', 'EPSG:4326', xs, ys)
+ ok = True
+ except Exception as e:
+ logger.debug(f"Coins GPS approximatifs (rasterio indisponible: {e})")
+ lons = None
+ lats = None
+ ok = False
+
+ for i, t in enumerate(tiles):
+ if ok:
+ corners = [[lats[4 * i], lons[4 * i]],
+ [lats[4 * i + 1], lons[4 * i + 1]],
+ [lats[4 * i + 2], lons[4 * i + 2]],
+ [lats[4 * i + 3], lons[4 * i + 3]]]
+ else:
+ corners = []
+ for cx in (t['col'] * 1000, (t['col'] + 1) * 1000):
+ for cy in ((t['row'] - 1) * 1000, t['row'] * 1000):
+ lon, lat = _approx_l93_to_wgs84(cx, cy)
+ corners.append([lat, lon])
+ # Reordonner SW, SE, NE, NW (la boucle donne SW, NW, SE, NE)
+ corners = [corners[0], corners[2], corners[3], corners[1]]
+ t['corners'] = corners
+ t['bounds'] = [[min(c[0] for c in corners), min(c[1] for c in corners)],
+ [max(c[0] for c in corners), max(c[1] for c in corners)]]
+ return ok
+
+
+def _mtime(path):
+ """Mtime d'un fichier, ou None si inaccessible."""
+ try:
+ return Path(path).stat().st_mtime
+ except OSError:
+ return None
+
+
+def _cached_file_fresh(path, src_mtime):
+ """True si un fichier cache existe et est plus récent que sa source.
+
+ Sert à invalider vignettes et sous-tuiles quand une tuile est recalculée :
+ l'image source (AVIF/WebP) étant réécrite, sa mtime devient plus récente
+ que celle du cache, qui doit alors être régénéré.
+ """
+ cached_mtime = _mtime(path)
+ if cached_mtime is None:
+ return False
+ return src_mtime is None or cached_mtime >= src_mtime
+
+
def generate_thumbnail(src_path, thumb_path, max_size=256):
"""Génère une vignette JPEG depuis une image AVIF/WebP existante.
@@ -246,11 +400,197 @@ def _pick_display_viz(viz_keys):
return sorted(viz_keys)[0]
-def build_index(output_dir, output_format='avif'):
- """Génère la carte continue HTML des tuiles traitées.
+def _subdivision_k(resolution, tile_m=1000, target_px=2500):
+ """Facteur de découpage k (grille k×k) pour alléger le rendu carte.
- Scanne output_dir/visualisations/, génère les vignettes JPEG, puis écrit
- output_dir/index.html (auto-suffisant) + output_dir/index_thumbs/.
+ Au 0,2 m/px une dalle de 1 km fait 5000×5000 px (~100 Mo décodés dans le
+ navigateur) : on la découpe en quadrants de 500 m (k=2, 2500×2500 px).
+ À 0,5 m/px (2000 px) la dalle reste entière (k=1).
+ """
+ px = max(1, int(round(tile_m / resolution)))
+ return max(1, int(round(px / target_px)))
+
+
+def _subtile_corners(corners, i, j, k):
+ """Coins WGS84 [SW, SE, NE, NW] de la sous-tuile (i, j) d'un découpage k×k.
+
+ i : indice vers l'est (0..k-1), j : indice vers le nord (0..k-1).
+ Interpolation bilinéaire des coins de la dalle — le quadrilatère projeté
+ est quasi un parallélogramme à cette échelle (erreur écran < 1 px).
+ """
+ sw, se, ne, nw = corners
+
+ def lerp(p, q, u):
+ return [p[0] + (q[0] - p[0]) * u, p[1] + (q[1] - p[1]) * u]
+
+ def at(u, v):
+ return lerp(lerp(sw, se, u), lerp(nw, ne, u), v)
+
+ u0, u1 = i / k, (i + 1) / k
+ v0, v1 = j / k, (j + 1) / k
+ return [at(u0, v0), at(u1, v0), at(u1, v1), at(u0, v1)]
+
+
+def _build_subtiles(tile, offered_viz_keys, output_dir, sub_dir_name):
+ """Découpe une dalle en sous-tuiles (crops AVIF) pour la carte interactive.
+
+ Ne découpe que les visualisations proposées dans le panneau de couches. Retourne
+ la liste des entrées display (une par sous-tuile), ou None si le découpage
+ n'est pas nécessaire/possible (la dalle entière sera alors affichée).
+ """
+ k = _subdivision_k(tile['resolution'])
+ if k <= 1:
+ return None
+ try:
+ from PIL import Image as PILImage
+ except ImportError:
+ return None
+
+ sub_dir = output_dir / sub_dir_name
+ sub_dir.mkdir(parents=True, exist_ok=True)
+
+ entries = {}
+ for j in range(k):
+ for i in range(k):
+ corners = _subtile_corners(tile['corners'], i, j, k)
+ entries[(i, j)] = {
+ 'col': tile['col'], 'row': tile['row'],
+ 'name': tile['name'], 'dir_name': tile['dir_name'],
+ 'resolution': tile['resolution'],
+ 'bounds': [[min(c[0] for c in corners), min(c[1] for c in corners)],
+ [max(c[0] for c in corners), max(c[1] for c in corners)]],
+ 'corners': corners,
+ 'display_viz': tile['display_viz'],
+ 'viz': {},
+ 'meta': tile.get('meta'),
+ 'sub_i': i, 'sub_j': j, 'sub_k': k,
+ 'size_km': round(1.0 / k, 3),
+ }
+
+ try:
+ resample = getattr(PILImage, 'LANCZOS', 1)
+ for viz_key in offered_viz_keys:
+ info = tile['viz'].get(viz_key)
+ if not info:
+ continue
+ stems = {key: f"{tile['dir_name']}_{viz_key}_{key[0]}_{key[1]}"
+ for key in entries}
+ # Régénère si au moins un fichier manque ou est périmé (dalle
+ # source recalculée depuis — comme les vignettes)
+ src = output_dir / info['full']
+ src_mtime = _mtime(src)
+ all_exist = all(
+ _cached_file_fresh(output_dir / sub_dir_name / (stem + '.avif'), src_mtime)
+ and _cached_file_fresh(output_dir / sub_dir_name / (stem + '_thumb.webp'), src_mtime)
+ for stem in stems.values())
+ if not all_exist:
+ logger.info(f" Sous-tuiles recalculées : {tile['dir_name']}/{viz_key} "
+ f"({len(stems)} découpages)")
+ try:
+ img = PILImage.open(str(src))
+ img.load()
+ except Exception as e:
+ logger.debug(f"Sous-tuilage {viz_key} impossible ({src.name}): {e}")
+ return None
+ if img.mode not in ('RGB', 'L'):
+ img = img.convert('RGB')
+ W, H = img.size
+ for (i, j), stem in stems.items():
+ # Image : ligne 0 = nord → la sous-tuile j (nord) est en haut
+ left, right = round(W * i / k), round(W * (i + 1) / k)
+ top = round(H * (1 - (j + 1) / k))
+ bottom = round(H * (1 - j / k))
+ quad = img.crop((left, top, right, bottom))
+ try:
+ quad.save(str(output_dir / sub_dir_name / (stem + '.avif')),
+ format='AVIF', quality=_SUBTILE_AVIF_QUALITY,
+ speed=_SUBTILE_AVIF_SPEED)
+ except Exception as e:
+ logger.warning(f"Encodage AVIF impossible ({stem}), "
+ f"sous-tuilage abandonné : {e}")
+ return None
+ # Supprime l'ancien crop .webp d'une génération précédente
+ (output_dir / sub_dir_name / (stem + '.webp')).unlink(missing_ok=True)
+ scale = min(1.0, 256 / max(quad.size))
+ if scale < 1.0:
+ quad = quad.resize((max(1, int(quad.size[0] * scale)),
+ max(1, int(quad.size[1] * scale))), resample)
+ quad.save(str(output_dir / sub_dir_name / (stem + '_thumb.webp')),
+ format='WEBP', quality=80)
+ for (i, j), stem in stems.items():
+ entries[(i, j)]['viz'][viz_key] = {
+ 'thumb': f"{sub_dir_name}/{stem}_thumb.webp",
+ 'full': f"{sub_dir_name}/{stem}.avif",
+ }
+ except Exception as e:
+ logger.warning(f"Sous-tuilage abandonné pour {tile['dir_name']}: {e}")
+ return None
+
+ usable = [e for e in entries.values() if e['viz']]
+ if not usable:
+ return None
+ # Repli dalle entière pour les visualisations non découpées : elles restent
+ # superposables en couche (images plus lourdes, mais fonctionnelles).
+ sub_keys = set(offered_viz_keys)
+ for entry in usable:
+ for viz_key, info in tile['viz'].items():
+ if viz_key not in sub_keys:
+ entry['viz'][viz_key] = dict(info)
+ return usable
+
+
+def _collect_tile_metadata(tile, dtm_dir):
+ """Rassemble les métadonnées de génération d'une tuile.
+
+ Lit la méthode de classification du sol depuis le sidecar DTM
+ (output/DTM/{basename}_dtm{suffix}_method.txt, écrit par pipeline.py),
+ et les dates/tailles des fichiers de visualisation.
+
+ Returns:
+ {method: str|None, generated: str|None,
+ viz: {viz_key: {date: str, size: int}}}
+ """
+ meta = {'method': None, 'generated': None, 'viz': {}}
+ suffix = _res_suffix_str(tile['resolution'])
+ method_file = Path(dtm_dir) / f"{tile['basename']}_dtm{suffix}_method.txt"
+ if not method_file.exists() and suffix:
+ # La classification du sol est partagée entre résolutions : repli sur
+ # le sidecar de la résolution primaire si le spécifique manque.
+ method_file = Path(dtm_dir) / f"{tile['basename']}_dtm_method.txt"
+ try:
+ if method_file.exists():
+ method = method_file.read_text(encoding='utf-8').strip()
+ if method:
+ meta['method'] = method
+ # Le sidecar est écrit juste après la création du DTM :
+ # sa date ≈ date de génération de la dalle.
+ meta['generated'] = datetime.fromtimestamp(
+ method_file.stat().st_mtime).strftime('%Y-%m-%d %H:%M')
+ except OSError as e:
+ logger.debug(f"Métadonnées illisibles {method_file.name}: {e}")
+
+ for viz_key, info in tile['viz'].items():
+ try:
+ st = (Path(tile['dir_path']) / info['filename']).stat()
+ meta['viz'][viz_key] = {
+ 'date': datetime.fromtimestamp(st.st_mtime).strftime('%Y-%m-%d %H:%M'),
+ 'size': st.st_size,
+ }
+ except OSError:
+ continue
+
+ if meta['generated'] is None and meta['viz']:
+ dates = [v['date'] for v in meta['viz'].values()]
+ meta['generated'] = min(dates)
+ return meta
+
+
+def build_index(output_dir, output_format='avif'):
+ """Génère la carte interactive HTML des tuiles traitées, organisée par zones.
+
+ Scanne output_dir/visualisations/, collecte les métadonnées de génération,
+ génère les vignettes JPEG, écrit les assets (CSS/JS) puis
+ output_dir/index.html.
Args:
output_dir: dossier de sortie racine (contient visualisations/).
@@ -261,112 +601,186 @@ def build_index(output_dir, output_format='avif'):
"""
output_dir = Path(output_dir)
vis_dir = output_dir / 'visualisations'
+ dtm_dir = output_dir / 'DTM'
+ t_start = time.time()
tiles = scan_tiles(vis_dir)
if not tiles:
logger.info("Aucune tuile traitée trouvée — index global non généré")
return None
- bbox = compute_bbox(tiles)
- assert bbox is not None # garanti par le test tiles non vide ci-dessus
+ # Bounds GPS par tuile (géoréférencement exact pour la carte Leaflet)
+ attach_gps_bounds(tiles)
+
+ # Une seule tuile par position (col, row) : on garde la résolution la plus
+ # fine disponible. Sinon les versions 0,5 m et 0,2 m d'une même dalle se
+ # superposent exactement sur la carte et celle ajoutée en dernier dans le
+ # DOM (la moins résolue, tri croissant) masque l'autre.
+ best_by_pos = {}
+ for t in tiles:
+ key = (t['col'], t['row'])
+ if key not in best_by_pos or t['resolution'] < best_by_pos[key]['resolution']:
+ best_by_pos[key] = t
+ tiles = sorted(best_by_pos.values(),
+ key=lambda t: (t['resolution'], -t['row'], t['col']))
+
+ # Détecter les zones géographiques (pour info uniquement)
+ zones = compute_zones(tiles)
+ logger.info(f" {len(zones)} zone(s) détectée(s)")
+
thumb_dir = output_dir / 'index_thumbs'
thumb_dir.mkdir(parents=True, exist_ok=True)
- # Collecte toutes les visualisations disponibles (pour le sélecteur).
+ # Collecte toutes les visualisations disponibles (pour le panneau de couches).
all_viz_keys = set()
for t in tiles:
all_viz_keys.update(t['viz'].keys())
+ # Restreint le découpage en sous-tuiles aux visualisations choisies
+ if _CARTO_SUBTILED_VIZ:
+ sub_viz = [v for v in _CARTO_SUBTILED_VIZ if v in all_viz_keys]
+ if sub_viz:
+ logger.info(f" Sous-tuilage limité à : {', '.join(sub_viz)}")
+ else:
+ sub_viz = list(all_viz_keys)
+
# Génère les vignettes et construit les données pour le HTML.
- tile_records = []
+ zone_records = []
thumbs_generated = 0
thumbs_failed = 0
- for t in tiles:
- viz_thumbs = {}
- for viz_key, info in t['viz'].items():
- src = Path(t['dir_path']) / info['filename']
- thumb_name = f"{t['dir_name']}_{viz_key}.jpg"
- thumb_path = thumb_dir / thumb_name
- # Régénère seulement si manquante
- if not thumb_path.exists():
- if generate_thumbnail(src, thumb_path):
- thumbs_generated += 1
+ tile_idx = 0
+ n_tiles = len(tiles)
+ logger.info(f" Vignettes : {n_tiles} tuile(s) × {len(all_viz_keys)} visualisation(s)")
+ for zone in zones:
+ zone_tile_records = []
+ for t in zone['tiles']:
+ tile_idx += 1
+ viz_thumbs = {}
+ regen = 0
+ for viz_key, info in t['viz'].items():
+ src = Path(t['dir_path']) / info['filename']
+ thumb_name = f"{t['dir_name']}_{viz_key}.jpg"
+ thumb_path = thumb_dir / thumb_name
+ # Régénère si manquante ou périmée (tuile recalculée depuis)
+ if not _cached_file_fresh(thumb_path, _mtime(src)):
+ if generate_thumbnail(src, thumb_path):
+ thumbs_generated += 1
+ regen += 1
+ else:
+ thumbs_failed += 1
+ continue
else:
- thumbs_failed += 1
- continue
- else:
- thumbs_generated += 1
- viz_thumbs[viz_key] = {
- 'thumb': f"index_thumbs/{thumb_name}",
- 'full': f"visualisations/{t['dir_name']}/{info['filename']}",
- }
+ thumbs_generated += 1
+ viz_thumbs[viz_key] = {
+ 'thumb': f"index_thumbs/{thumb_name}",
+ 'full': f"visualisations/{t['dir_name']}/{info['filename']}",
+ }
- if not viz_thumbs:
- continue
+ if regen:
+ logger.info(f" [{tile_idx}/{n_tiles}] {t['dir_name']} — "
+ f"{regen} vignette(s) régénérée(s)")
- display_viz = _pick_display_viz(viz_thumbs.keys())
- tile_records.append({
- 'col': t['col'],
- 'row': t['row'],
- 'name': t['basename'],
- 'dir_name': t['dir_name'],
- 'resolution': t['resolution'],
- 'display_viz': display_viz,
- 'viz': viz_thumbs,
- })
+ if not viz_thumbs:
+ continue
- if not tile_records:
+ display_viz = _pick_display_viz(viz_thumbs.keys())
+ tile_meta = _collect_tile_metadata(t, dtm_dir)
+ zone_tile_records.append({
+ 'col': t['col'],
+ 'row': t['row'],
+ 'name': t['basename'],
+ 'dir_name': t['dir_name'],
+ 'resolution': t['resolution'],
+ 'bounds': t.get('bounds'),
+ 'corners': t.get('corners'),
+ 'display_viz': display_viz,
+ 'viz': viz_thumbs,
+ 'meta': tile_meta,
+ })
+
+ if zone_tile_records:
+ zone_records.append({
+ 'label': zone['label'],
+ 'tiles': zone_tile_records,
+ 'bbox': zone['bbox'],
+ })
+
+ if not zone_records:
logger.warning("Aucune vignette générée — index global abandonné")
return None
- # HTML avec données intégrées.
- html = _render_html(tile_records, bbox, all_viz_keys, output_format)
+ # Calculer le bbox global (pour la stat bar)
+ global_bbox = compute_bbox(tiles)
+
+ # HTML avec données intégrées : liste plate des quads affichables.
+ # Les dalles 0,2 m (5000×5000 px) sont découpées en sous-tuiles 500 m
+ # (quadrants 2500×2500 px) pour alléger mémoire et chargements navigateur.
+ sub_dir_name = 'index_subtiles'
+ display_tiles = []
+ n_dalles = 0
+ n_sous = 0
+ for zr in zone_records:
+ for t in zr['tiles']:
+ n_dalles += 1
+ subs = _build_subtiles(t, sub_viz, output_dir, sub_dir_name)
+ if subs:
+ display_tiles.extend(subs)
+ n_sous += len(subs)
+ else:
+ display_tiles.append(t)
+ if n_sous:
+ logger.info(f" {n_sous} sous-tuiles générée(s) pour {n_dalles} dalle(s)")
+
+ _write_assets(output_dir)
+ html = _render_html(display_tiles, global_bbox, all_viz_keys, output_format)
html_path = output_dir / 'index.html'
html_path.write_text(html, encoding='utf-8')
- logger.info(f"Index global généré : {html_path}")
- logger.info(f" {len(tile_records)} tuile(s) • {thumbs_generated} vignette(s) générée(s)"
+ logger.info(f"Index global généré : {html_path} ({time.time() - t_start:.1f}s)")
+ logger.info(f" {len(tiles)} tuile(s) • {thumbs_generated} vignette(s) générée(s)"
+ (f" • {thumbs_failed} échec(s)" if thumbs_failed else ""))
- logger.info(f" Grille : {bbox['min_col']}-{bbox['max_col']} km E × "
- f"{bbox['min_row']}-{bbox['max_row']} km N")
+ logger.info(f" Grille : {global_bbox['min_col']}-{global_bbox['max_col']} km E × "
+ f"{global_bbox['min_row']}-{global_bbox['max_row']} km N")
return html_path
-def _render_html(tile_records, bbox, all_viz_keys, output_format):
- """Construit le HTML complet avec CSS et JS (carte continue + modal PDF)."""
- # Ordre des viz dans le sélecteur (selon ordre préféré puis alpha).
+def _write_assets(output_dir):
+ """Écrit les fichiers statiques de l'interface (CSS + JS) dans output/assets/."""
+ assets_dir = Path(output_dir) / 'assets'
+ assets_dir.mkdir(parents=True, exist_ok=True)
+ (assets_dir / 'app.css').write_text(_APP_CSS, encoding='utf-8')
+ (assets_dir / 'app.js').write_text(_APP_JS, encoding='utf-8')
+
+
+def _render_html(tiles, global_bbox, all_viz_keys, output_format):
+ """Construit la coquille HTML (données intégrées, styles/JS dans assets/)."""
+ # Ordre des couches dans le panneau (selon ordre préféré puis alpha).
ordered_viz = [v for v in _VIZ_FALLBACK_ORDER if v in all_viz_keys]
for v in sorted(all_viz_keys):
if v not in ordered_viz:
ordered_viz.append(v)
- data_json = json.dumps({
- 'tiles': tile_records,
- 'bbox': bbox,
- 'vizList': ordered_viz,
- }, ensure_ascii=False)
+ tiles_json = json.dumps(tiles, ensure_ascii=False)
viz_meta_json = json.dumps({
- k: {'label': VIZ_LABELS.get(k, k), 'cmap': VIZ_COLORMAPS.get(k, {}).get('cmap', 'terrain'),
- 'diverging': VIZ_COLORMAPS.get(k, {}).get('diverging', False)}
+ k: {'label': VIZ_LABELS.get(k, k)}
for k in ordered_viz
}, ensure_ascii=False)
- options_html = '\n'.join(
- f' '
- for v in ordered_viz
- )
+ # Couche activée par défaut : DEFAULT_VIZ si proposée, sinon la 1re
+ default_viz = DEFAULT_VIZ if DEFAULT_VIZ in ordered_viz else (ordered_viz[0] if ordered_viz else DEFAULT_VIZ)
+ stats_json = json.dumps({'default_viz': default_viz, 'n_tiles': len(tiles)},
+ ensure_ascii=False)
- n_tiles = len(tile_records)
- grid_w = bbox['max_col'] - bbox['min_col'] + 1
- grid_h = bbox['max_row'] - bbox['min_row'] + 1
+ n_tiles = len(tiles)
+ grid_w = global_bbox['max_col'] - global_bbox['min_col'] + 1
+ grid_h = global_bbox['max_row'] - global_bbox['min_row'] + 1
return _HTML_TEMPLATE.format(
- data_json=data_json,
+ tiles_json=tiles_json,
viz_meta_json=viz_meta_json,
- options_html=options_html,
+ stats_json=stats_json,
n_tiles=n_tiles,
grid_w=grid_w,
grid_h=grid_h,
@@ -379,546 +793,1077 @@ _HTML_TEMPLATE = """
-Carte continue LiDAR
-
+Carte LiDAR
+
+
-
-
Carte continue LiDAR
- {n_tiles} tuile(s) • {grid_w}×{grid_h} km
-
-
-
-