Isoler les dalles voisines de raccord dans input/edge_neighbors/

This commit is contained in:
Antoine Jacquin
2026-09-19 23:43:50 +02:00
parent 5ef0232c80
commit 3236d56a4b
4 changed files with 50 additions and 12 deletions

View File

@ -32,7 +32,7 @@
- **Calage vertical des faisceaux de vol** : chaque tuile mélange plusieurs passes (1-2 `PointSourceId` par passe) parfois biaisées verticalement de quelques cm (±2,5 cm mesurés sur 1000_6882). `create_dtm_fast` mesure l'offset robuste de chaque faisceau (points sol, maille 1 m, surface médiane itérée 3×) et retranche les offsets ≥ 0,5 cm (`STRIP_ALIGN_THRESHOLD` dans `dtm.py`) avant rastérisation. Offsets calculés **par tuile** (ils dérivent le long d'une ligne de vol : jamais de table globale), mémoïsés par LAS sol, consignés dans `DTM/*_dtm*_stripalign.json` (sidecar de cache : absent, ou version/seuil/paramètres différents ⇒ régénération du DTM). Désactivable : `--no-strip-align`.
- **Gigue intra-faisceau (2ᵉ passe du calage)** : les lignes de balayage successives d'une MÊME passe peuvent être décalées verticalement de façon aléatoire (vibration capteur / bruit haute fréquence de trajectoire) — un offset constant par faisceau n'y suffit pas. `_strip_jitter_offsets` découpe chaque faisceau en fenêtres de temps GPS (`STRIP_JITTER_BIN` = 0,1 s, origine de temps propre à chaque faisceau), mesure l'offset robuste de chaque fenêtre contre la surface médiane des AUTRES faisceaux (maille 1 m partagée, ≥ `STRIP_JITTER_MIN_CELLS` = 40 cellules), lisse la série (médiane glissante `STRIP_JITTER_SMOOTH` = 5 fenêtres), la borne à ± `STRIP_JITTER_MAX` (10 cm) puis l'interpole au temps GPS de chaque point (`_apply_strip_jitter`) ; en recouvrement à deux faisceaux, chacun reçoit une série (chacun absorbe sa part). Requiert la dimension `gps_time` (silencieusement ignorée sinon). Sidecar version 2 (séries dans `jitter`), couverte par `--no-strip-align`.
- **Openness sous-échantillonnée** : `generate_openness` calcule le lancé de rayons (l'étape la plus coûteuse : 532 s/tuile à 0,2 m sur CPU) sur une grille décimée par blocs (`OPENNESS_DOWNSAMPLE = 2` : max par bloc en positive, min en négative — préserve les reliefs qui bornent l'horizon) puis rééchantillonne en bilinéaire. Coût ÷ facteur³ : 532 s → 40 s (×13). Signal archéologique préservé (corr. 0,93 après lissage) ; la texture de bruit sub-métrique disparaît. `--openness-downsample 1` = pleine résolution. SVF et openness anisotrope ne sont PAS concernés.
- **Raccord des bords entre tuiles** : les rendus à grand noyau (openness/SVF : rayons 100 m ; LRM : 15 m) tronquent leur fenêtre au bord de dalle — bandes d'artefacts à chaque changement de tuile. `--edge-buffer N` (défaut 0 = off ; case « Raccord des bords » de la webapp, `EDGE_BUFFER_METERS` = 100 m dans `webapp.py`) fait rastériser le DTM sur la **dalle nominale 1 km alignée sur la grille** plus une bande de N m remplie avec les points sol des 8 LAZ voisines (`_neighbor_ground_points` dans `dtm.py` : PDAL en flux, découpe + filtre de classes IGN ; voisine absente de input/ = téléchargement automatique depuis le catalogue IGN avant le run (dédupliqué sur tout le lot, `_fetch_edge_neighbors` dans `pipeline.py`) ; introuvable ou échec = bande vide). Les visualisations calculent sur l'emprise étendue puis `rendering.py` (`_core_tile_window`, via `tif_to_crop`/`tif_to_png`) recadre les sorties sur la dalle 1 km exacte lue dans le nom LHD — les AVIF restent des carrés 1 km alignés dans la mosaïque. Tampon consigné dans le tag GeoTIFF `LIDAR_EDGE_BUFFER` du DTM : changer `--edge-buffer` invalide le cache DTM automatiquement (tag absent = 0). Bandes voisines non calées par faisceaux (contexte seul, recadrée hors image finale). Coût : ~7 s de lecture par voisine + ~44 % de pixels en plus à 100 m/0,2 m. Nom hors pattern LHD : option ignorée (bornes d'en-tête, pas de recadrage).
- **Raccord des bords entre tuiles** : les rendus à grand noyau (openness/SVF : rayons 100 m ; LRM : 15 m) tronquent leur fenêtre au bord de dalle — bandes d'artefacts à chaque changement de tuile. `--edge-buffer N` (défaut 0 = off ; case « Raccord des bords » de la webapp, `EDGE_BUFFER_METERS` = 100 m dans `webapp.py`) fait rastériser le DTM sur la **dalle nominale 1 km alignée sur la grille** plus une bande de N m remplie avec les points sol des 8 LAZ voisines (`_neighbor_ground_points` dans `dtm.py` : PDAL en flux, découpe + filtre de classes IGN ; voisine absente = téléchargement automatique depuis le catalogue IGN avant le run, **isolée dans `input/edge_neighbors/`** pour ne pas gonfler le corpus des passes globales (dédupliqué sur tout le lot, `_fetch_edge_neighbors` dans `pipeline.py`) ; introuvable ou échec = bande vide). Les visualisations calculent sur l'emprise étendue puis `rendering.py` (`_core_tile_window`, via `tif_to_crop`/`tif_to_png`) recadre les sorties sur la dalle 1 km exacte lue dans le nom LHD — les AVIF restent des carrés 1 km alignés dans la mosaïque. Tampon consigné dans le tag GeoTIFF `LIDAR_EDGE_BUFFER` du DTM : changer `--edge-buffer` invalide le cache DTM automatiquement (tag absent = 0). Bandes voisines non calées par faisceaux (contexte seul, recadrée hors image finale). Coût : ~7 s de lecture par voisine + ~44 % de pixels en plus à 100 m/0,2 m. Nom hors pattern LHD : option ignorée (bornes d'en-tête, pas de recadrage).
- **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.

View File

@ -1043,8 +1043,14 @@ def _interpolate_holes(dtm, downsample=8):
# images finales sont recadrées sur la dalle exacte (cf. rendering.py).
EDGE_BUFFER_TAG = "LIDAR_EDGE_BUFFER" # tag GeoTIFF : tampon utilisé (m)
# Sous-dossier de input/ où sont isolées les dalles voisines téléchargées
# pour le seul raccord des bords : elles ne font pas partie du corpus de
# tuiles à rendre (les scans globaux n'énumèrent que input/ à plat).
EDGE_NEIGHBORS_DIRNAME = "edge_neighbors"
# Décalages des 8 voisines d'une dalle (col, row) en km
_NEIGHBOR_OFFSETS = [(-1, -1), (0, -1), (1, -1), (-1, 0),
(1, 0), (-1, 1), (0, 1), (1, 1)]
@ -1055,12 +1061,19 @@ def _tile_coords(name):
def _neighbor_laz_files(source_laz):
"""Liste les 8 fichiers LAZ/LAS adjacents à `source_laz` dans son dossier."""
"""Liste les 8 fichiers LAZ/LAS adjacents à `source_laz`.
Recherche dans le dossier de la source, puis dans le sous-dossier
edge_neighbors/ : les dalles voisines téléchargées pour le seul raccord
des bords y sont isolées — sinon elles s'accumulent à plat dans input/ et
chaque passe globale annexe un anneau de tuiles à rendre en plus.
"""
coords = _tile_coords(source_laz)
if coords is None:
return []
col, row = coords
directory = Path(source_laz).parent
search_dirs = [directory, directory / EDGE_NEIGHBORS_DIRNAME]
neighbors = []
for dcol, drow in _NEIGHBOR_OFFSETS:
nc, nr = col + dcol, row + drow
@ -1068,9 +1081,10 @@ def _neighbor_laz_files(source_laz):
# on accepte aussi la variante sans remplissage pour les dalles exotiques.
names = {f"LHD_FXX_{nc:04d}_{nr:04d}", f"LHD_FXX_{nc}_{nr}"}
matches = []
for search_dir in search_dirs:
for name in names:
matches += list(directory.glob(f"{name}_*.las"))
matches += list(directory.glob(f"{name}_*.laz"))
matches += list(search_dir.glob(f"{name}_*.las"))
matches += list(search_dir.glob(f"{name}_*.laz"))
if matches:
neighbors.append(sorted(matches)[0])
else:

View File

@ -456,17 +456,20 @@ class LidarArchaeoPipeline:
"""Télécharge les dalles LAZ voisines manquantes (raccord des bords).
La bande de raccord lit les 8 LAZ adjacentes de chaque tuile ; celles
absentes de input/ sont téléchargées depuis le catalogue IGN avant de
lancer les workers, pour que la bande soit remplie jusqu'au bord de la
zone même pour des tuiles jamais rendues. La liste est dédupliquée sur
absentes sont téléchargées depuis le catalogue IGN avant de lancer
les workers, dans le sous-dossier input/edge_neighbors/ : elles ne
doivent PAS rejoindre input/ à plat, sinon les passes globales
(« tout input/ ») les comptent comme des tuiles à rendre et la zone
grandit d'un anneau à chaque relance. La liste est dédupliquée sur
tout le lot : la couronne d'un bloc contigu ne coûte qu'un passage.
Une dalle introuvable (zone non publiée) ou en échec laisse simplement
la bande vide — le rendu continue (best-effort).
"""
if self.edge_buffer <= 0:
return
from .dtm import _tile_coords, _NEIGHBOR_OFFSETS
from .dtm import _tile_coords, _NEIGHBOR_OFFSETS, EDGE_NEIGHBORS_DIRNAME
from .fetch_ign import fetch_tiles, tile_filename
edge_dir = self.input_dir / EDGE_NEIGHBORS_DIRNAME
wanted = set()
for laz_file in files:
coords = _tile_coords(Path(laz_file).name)
@ -477,16 +480,19 @@ class LidarArchaeoPipeline:
wanted.add((col + dcol, row + drow))
missing = sorted(
cr for cr in wanted
if not (self.input_dir / tile_filename(cr[0], cr[1])).exists())
if not (self.input_dir / tile_filename(cr[0], cr[1])).exists()
and not (edge_dir / tile_filename(cr[0], cr[1])).exists())
if not missing:
return
edge_dir.mkdir(parents=True, exist_ok=True)
logger.info(f"Raccord des bords : {len(missing)} dalle(s) voisine(s) "
f"absente(s) de input/ — téléchargement depuis le catalogue IGN")
f"absente(s) — téléchargement dans {EDGE_NEIGHBORS_DIRNAME}/ "
f"(catalogue IGN)")
t0 = time.time()
from concurrent.futures import ThreadPoolExecutor
with ThreadPoolExecutor(max_workers=4) as pool:
fetched = [path for batch in pool.map(
lambda spec: fetch_tiles(self.input_dir, [spec]), missing)
lambda spec: fetch_tiles(edge_dir, [spec]), missing)
for path in batch]
logger.info(f"Raccord des bords : {len(fetched)}/{len(missing)} dalle(s) "
f"voisine(s) téléchargée(s) ({time.time() - t0:.0f} s)")

View File

@ -720,6 +720,24 @@ class TestEdgeBuffer:
src.touch()
assert _neighbor_laz_files(src) == []
def test_neighbor_found_in_edge_subdir(self, tmp_output_dir):
"""Une voisine isolée dans edge_neighbors/ est trouvée ; la priorité
reste à une dalle à plat dans input/."""
from lidar_pipeline.dtm import _neighbor_laz_files, EDGE_NEIGHBORS_DIRNAME
base = "LHD_FXX_0637_6627_PTS_LAMB93_IGN69"
(tmp_output_dir / f"{base}.copc.laz").touch()
edge = tmp_output_dir / EDGE_NEIGHBORS_DIRNAME
edge.mkdir()
# Voisine uniquement dans le sous-dossier de raccord
(edge / "LHD_FXX_0638_6628_PTS_LAMB93_IGN69.copc.laz").touch()
# Voisine présente aux deux endroits : la version input/ gagne
(edge / "LHD_FXX_0636_6626_PTS_LAMB93_IGN69.copc.laz").touch()
(tmp_output_dir / "LHD_FXX_0636_6626_PTS_LAMB93_IGN69.copc.laz").touch()
found = _neighbor_laz_files(tmp_output_dir / f"{base}.copc.laz")
by_name = {f.name: f for f in found}
assert by_name["LHD_FXX_0638_6628_PTS_LAMB93_IGN69.copc.laz"].parent == edge
assert by_name["LHD_FXX_0636_6626_PTS_LAMB93_IGN69.copc.laz"].parent == tmp_output_dir
def test_buffered_dtm_extends_into_neighbor(self, tmp_output_dir):
"""MNT 24x24 (dalle 20x20 + bande 100 m), bande EST remplie à z=20 par la voisine."""
from lidar_pipeline.dtm import create_dtm_fast, read_dtm_edge_buffer, EDGE_BUFFER_TAG