diff --git a/AGENTS.md b/AGENTS.md index 039b14d..b980c1b 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -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. diff --git a/lidar_pipeline/dtm.py b/lidar_pipeline/dtm.py index c125906..275c852 100644 --- a/lidar_pipeline/dtm.py +++ b/lidar_pipeline/dtm.py @@ -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 name in names: - matches += list(directory.glob(f"{name}_*.las")) - matches += list(directory.glob(f"{name}_*.laz")) + for search_dir in search_dirs: + for name in names: + matches += list(search_dir.glob(f"{name}_*.las")) + matches += list(search_dir.glob(f"{name}_*.laz")) if matches: neighbors.append(sorted(matches)[0]) else: diff --git a/lidar_pipeline/pipeline.py b/lidar_pipeline/pipeline.py index 7321110..7b8f624 100644 --- a/lidar_pipeline/pipeline.py +++ b/lidar_pipeline/pipeline.py @@ -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)") diff --git a/lidar_pipeline/tests/test_dtm.py b/lidar_pipeline/tests/test_dtm.py index 1f13b8c..9becae8 100644 --- a/lidar_pipeline/tests/test_dtm.py +++ b/lidar_pipeline/tests/test_dtm.py @@ -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