From 5f018fe2b65970ea8c8cfd0c59ef862027c9f976 Mon Sep 17 00:00:00 2001 From: Antoine Jacquin Date: Sat, 19 Sep 2026 20:44:44 +0200 Subject: [PATCH] Raccorder les bords de tuiles voisines et encadrer les rendus en cours MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Le MNT de chaque dalle s'étend d'une bande de 100 m remplie avec les points sol des 8 tuiles voisines (option « Raccord des bords ») : les rendus à grand noyau (openness, SVF, LRM) deviennent continus d'une dalle à l'autre, les images restent recadrées sur le kilomètre exact. Pendant une génération, la carte encadre les dalles du run : orange pulsant en cours de rendu, rouge en échec — les coins WGS84 sont portés par /api/status pour toutes les dalles non terminées. La file de génération trie les dalles du nord au sud et passe à 10 workers GPU. --- AGENTS.md | 5 +- docker-compose.local-2m.yml | 2 +- docker-compose.worker.yml | 2 +- docker-compose.yml | 2 +- docs/GROUND_CLASSIFICATION.md | 11 +- lidar_pipeline/cli.py | 27 +- lidar_pipeline/dtm.py | 569 ++++++++++++++++++++----- lidar_pipeline/index.py | 58 ++- lidar_pipeline/pipeline.py | 122 +++++- lidar_pipeline/rendering.py | 50 ++- lidar_pipeline/tests/test_dtm.py | 263 ++++++++---- lidar_pipeline/tests/test_pipeline.py | 65 +++ lidar_pipeline/tests/test_rendering.py | 43 +- lidar_pipeline/tests/test_webapp.py | 35 +- lidar_pipeline/webapp.py | 68 ++- 15 files changed, 1075 insertions(+), 247 deletions(-) diff --git a/AGENTS.md b/AGENTS.md index f7057a9..23d84ec 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -20,6 +20,7 @@ ## Conventions - **Generation is 0.2 m only** (policy): `/api/generate` (`GENERATE_RESOLUTIONS` in `webapp.py`), the compose `process` command and the CLI `-r` default all produce 0.2 m exclusively; 0.5 m stays available via explicit `-r 0.5`. Completeness detection (`complete_cells`) requires the viz at 0.2 m only. +- **Génération du nord au sud** : les tuiles sont traitées par ligne décroissante (row = nord en km), colonnes croissantes — `find_laz_files` (pipeline.py) pour les passes batch et `_resolve_request` (webapp.py) pour les runs lancés depuis la carte. Les workers prennent les fichiers dans l'ordre de soumission : la carte se remplit de haut en bas pendant un run (`--file` explicite au CLI = ordre utilisateur préservé). Parallélisme de génération : `LIDAR_WORKERS` (10 dans les compose ; fallback 10 dans `webapp.py`). - **Sub-tuilage intégral** : `_CARTO_SUBTILED_VIZ` (vide dans `index.py`) découpe TOUTES les couches en quadrants 500 m à 0,2 m ; ortho/topo sont encodées en AVIF q75 (`_SUBTILE_DETAIL_VIZ`) contre q55 pour les rampes de couleur. Une couche qui échoue à la découpe retombe en dalle entière (`_fallback_full_dalle`) sans pénaliser les autres. - **Bilingual naming**: all code identifiers are English; every user-facing string, log message, argparse help, and comment is French. - **Adding a visualization requires 4 edits**: (1) `generate_X()` in `visualizations.py`, (2) entry in `VIZ_STEPS` in `pipeline.py`, (3) entry in `COLORMAPS` in `rendering.py`, (4) entry in `VIZ_LEGENDS` in `index.py` (title/legend/description + sampled cmap gradient — single text source merged into `COLORMAPS` at import, also used by the export mosaic legend in `export.py`). Missing any one breaks the pipeline. @@ -28,8 +29,10 @@ - **Logger is always `logging.getLogger("lidar")`**, never `__name__`. All modules route through this single logger so worker processes can configure it. - **Filename special-cases** in `_expected_output_path()`: `pos_open` → `positive_openness`, `neg_open` → `negative_openness`, `hillshade` → `hillshade_multi`. - **Default output is AVIF**, not WebP. Use `--format webp` for WebP. Quality default is 98. -- **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 différents ⇒ régénération du DTM), appliqués aussi au plancher `--bare-earth`. Désactivable : `--no-strip-align`. +- **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). - **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/docker-compose.local-2m.yml b/docker-compose.local-2m.yml index d8f44a9..219bb39 100644 --- a/docker-compose.local-2m.yml +++ b/docker-compose.local-2m.yml @@ -35,7 +35,7 @@ services: - LIDAR_INPUT_DIR=/data/input - LIDAR_OUTPUT_DIR=/data/output - LIDAR_GPU=1 - - LIDAR_WORKERS=2 + - LIDAR_WORKERS=10 command: python3 -m uvicorn lidar_pipeline.webapp:app --host 0.0.0.0 --port 8973 restart: unless-stopped diff --git a/docker-compose.worker.yml b/docker-compose.worker.yml index b645b29..2520cf6 100644 --- a/docker-compose.worker.yml +++ b/docker-compose.worker.yml @@ -33,7 +33,7 @@ services: - LIDAR_OUTPUT_DIR=/data/output # Les générations lancées depuis une webapp distante utilisent le GPU - LIDAR_GPU=1 - - LIDAR_WORKERS=5 + - LIDAR_WORKERS=10 # Protéger l'API si le réseau n'est pas de confiance : même valeur que # LIDAR_REMOTE_TOKEN sur chaque webapp distante (sinon, laisser commenté) # - LIDAR_API_TOKEN=change-moi diff --git a/docker-compose.yml b/docker-compose.yml index 6d9b522..a677aa8 100644 --- a/docker-compose.yml +++ b/docker-compose.yml @@ -30,7 +30,7 @@ services: - LIDAR_OUTPUT_DIR=/data/output # Les générations lancées depuis la carte utilisent le GPU - LIDAR_GPU=1 - - LIDAR_WORKERS=2 + - LIDAR_WORKERS=10 # Machine de traitement pour une webapp distante (Raspberry Pi) : # décommenter et définir la même valeur en LIDAR_REMOTE_TOKEN là-bas # (protège /api/generate, /api/preview, /api/rebuild, /api/sync). diff --git a/docs/GROUND_CLASSIFICATION.md b/docs/GROUND_CLASSIFICATION.md index aafd4fe..0efeac4 100644 --- a/docs/GROUND_CLASSIFICATION.md +++ b/docs/GROUND_CLASSIFICATION.md @@ -113,12 +113,11 @@ interpolation. Le KPConv/RandLA-Net est plus précis en 3D pur mais plus lourd - 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. + - Le MNT n'est complété que pour les petits trous (< 1 m, `fillnodata`) : + les grands trous (forêt dense, relief raide où le sol est sous-classé) + restent en nodata (noir dans les rendus). Volontairement pas de plancher + au retour le plus bas : sous canopée dense ce retour est la végétation, + qui imprimerait les arbres dans le MNT. - 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 diff --git a/lidar_pipeline/cli.py b/lidar_pipeline/cli.py index 0c4f90e..960c68d 100644 --- a/lidar_pipeline/cli.py +++ b/lidar_pipeline/cli.py @@ -143,13 +143,6 @@ def main(): "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( "--openness-downsample", type=int, @@ -159,13 +152,27 @@ def main(): "grille décimée par blocs puis rééchantillonnage. Défaut : 2 (~8× plus " "rapide, rendu quasi identique) ; 1 = pleine résolution" ) + parser.add_argument( + "--edge-buffer", + type=float, + default=0.0, + metavar="METRES", + help="Raccord des bords : étendre le MNT d'une bande de N mètres remplie avec " + "les points sol des 8 tuiles LAZ voisines, pour que les rendus à grand " + "noyau (openness, SVF, LRM) soient continus d'une tuile à l'autre. " + "100 m couvre tous les rayons ; les images finales sont recadrées sur " + "la dalle 1 km exacte. 0 = désactivé (défaut). Un changement de valeur " + "régénère les DTM concernés." + ) parser.add_argument( "--no-strip-align", action="store_true", help="Désactiver le calage vertical des faisceaux de vol. Par défaut, les écarts " "verticaux ≥ 0,5 cm entre lignes de vol (PointSourceId) d'une tuile sont mesurés " - "sur les points sol et corrigés avant rastérisation du MNT (offsets consignés " - "dans DTM/*_stripalign.json)" + "sur les points sol et corrigés avant rastérisation du MNT, puis la gigue " + "intra-faisceau (décalages aléatoires des lignes de balayage successives d'une " + "même passe, type vibration capteur) est corrigée par fenêtres de temps GPS " + "(offsets et séries consignés dans DTM/*_stripalign.json)" ) parser.add_argument( "--keep-tif", @@ -360,9 +367,9 @@ def main(): ign_classes=args.ign_classes, force_classify=args.force_classification, keep_tif=args.keep_tif, - bare_earth=args.bare_earth, strip_align=not args.no_strip_align, openness_downsample=args.openness_downsample, + edge_buffer=args.edge_buffer, quality=quality, only_viz=only_viz, skip_viz=skip_viz, diff --git a/lidar_pipeline/dtm.py b/lidar_pipeline/dtm.py index 0da25a9..c125906 100644 --- a/lidar_pipeline/dtm.py +++ b/lidar_pipeline/dtm.py @@ -42,14 +42,74 @@ IGN_CLASS_NAMES = { # selon la campagne, rien n'est codé en dur) et on le retranche avant # rastérisation. Les offsets dérivent le long d'une ligne de vol (signe inversé # entre tuiles voisines mesuré) : le calcul est donc par tuile, jamais global. -STRIP_ALIGN_VERSION = 1 +# +# 2ᵉ passe — gigue intra-faisceau : le décalage vertical peut aussi varier au +# fil d'une MÊME passe (vibration capteur / bruit haute fréquence de la +# trajectoire) : les lignes de balayage successives d'un faisceau apparaissent +# alors décalées verticalement de façon aléatoire. Chaque faisceau est découpé +# en fenêtres de temps GPS, l'offset robuste de chaque fenêtre est mesuré contre +# la surface médiane des autres faisceaux (même maille 1 m), la série est lissée +# (médiane glissante) puis interpolée au temps GPS de chaque point. Requiert la +# dimension gps_time (ignorée silencieusement sinon). +STRIP_ALIGN_VERSION = 2 STRIP_ALIGN_THRESHOLD = 0.005 # m : écart mini pour corriger un faisceau (0,5 cm) STRIP_ALIGN_CELL = 1.0 # m : maille de comparaison des faisceaux STRIP_ALIGN_MIN_SHARED = 500 # cellules sol communes mini pour valider un offset +STRIP_JITTER_BIN = 0.1 # s : durée d'une fenêtre de temps GPS (gigue) +STRIP_JITTER_SMOOTH = 5 # fenêtres : largeur de la médiane glissante +STRIP_JITTER_MIN_CELLS = 40 # cellules sol partagées mini pour valider une fenêtre +STRIP_JITTER_MAX = 0.10 # m : amplitude maxi d'une correction de gigue (garde-fou) # Mémo des offsets par fichier : la classification est partagée entre # résolutions, le même LAS sol est rasterisé à 0,5 m puis 0,2 m. _STRIP_OFFSETS_CACHE = {} +_STRIP_JITTER_CACHE = {} + + +def _strip_surface_grid(x, y, z, inv, n_sources, cell=STRIP_ALIGN_CELL): + """Surfaces sol par faisceau sur maille 1 m (moyenne des points bas). + + Pour chaque faisceau et chaque maille : moyenne des points situés à moins + de 0,5 m du minimum du faisceau dans la maille (robuste à la végétation + résiduelle), mailles à ≥ 2 points seulement. + + Returns: + (allcells, grid, key) : cellules triées, grid[faisceau, cellule] = Z + moyen (NaN si absent) et clé de maille de CHAQUE point ; None si + aucune surface n'est peuplée. + """ + x0 = np.floor(np.min(x) / cell) * cell + y0 = np.floor(np.min(y) / cell) * cell + xi = ((x - x0) / cell).astype(np.int64) + yi = ((y - y0) / cell).astype(np.int64) + ny = int(yi.max()) + 1 + key = xi * ny + yi + + def _surface(k): + m = inv == k + kk, zz = key[m], z[m] + order = np.argsort(kk, kind='stable') + k_s, z_s = kk[order], zz[order] + starts = np.flatnonzero(np.r_[True, k_s[1:] != k_s[:-1]]) + mins = np.minimum.reduceat(z_s, starts) + ukey = k_s[starts] + thr = mins[np.searchsorted(ukey, kk)] + sel = np.flatnonzero(zz <= thr + 0.5) + k2, z2 = kk[sel], zz[sel] + cnt = np.bincount(k2) + sums = np.bincount(k2, weights=z2) + v = np.flatnonzero(cnt >= 2) + return v, sums[v] / cnt[v] + + surfaces = [_surface(k) for k in range(n_sources)] + populated = [c for c, _ in surfaces if len(c)] + if not populated: + return None + allcells = np.unique(np.concatenate(populated)) + grid = np.full((n_sources, len(allcells)), np.nan) + for k, (c, zs) in enumerate(surfaces): + grid[k, np.searchsorted(allcells, c)] = zs + return allcells, grid, key def _strip_vertical_offsets(x, y, z, psid, cell=STRIP_ALIGN_CELL, @@ -78,37 +138,10 @@ def _strip_vertical_offsets(x, y, z, psid, cell=STRIP_ALIGN_CELL, us, inv = np.unique(psid, return_inverse=True) if len(us) < 2: return {} - x0 = np.floor(np.min(x) / cell) * cell - y0 = np.floor(np.min(y) / cell) * cell - xi = ((x - x0) / cell).astype(np.int64) - yi = ((y - y0) / cell).astype(np.int64) - ny = int(yi.max()) + 1 - key = xi * ny + yi - - def _surface(k): - m = inv == k - kk, zz = key[m], z[m] - order = np.argsort(kk, kind='stable') - k_s, z_s = kk[order], zz[order] - starts = np.flatnonzero(np.r_[True, k_s[1:] != k_s[:-1]]) - mins = np.minimum.reduceat(z_s, starts) - ukey = k_s[starts] - thr = mins[np.searchsorted(ukey, kk)] - sel = np.flatnonzero(zz <= thr + 0.5) - k2, z2 = kk[sel], zz[sel] - cnt = np.bincount(k2) - sums = np.bincount(k2, weights=z2) - v = np.flatnonzero(cnt >= 2) - return v, sums[v] / cnt[v] - - surfaces = [_surface(k) for k in range(len(us))] - populated = [c for c, _ in surfaces if len(c)] - if not populated: + built = _strip_surface_grid(x, y, z, inv, len(us), cell) + if built is None: return {} - allcells = np.unique(np.concatenate(populated)) - grid = np.full((len(us), len(allcells)), np.nan) - for k, (c, zs) in enumerate(surfaces): - grid[k, np.searchsorted(allcells, c)] = zs + allcells, grid, _key = built comparable = np.sum(~np.isnan(grid), axis=0) >= 2 offsets = np.zeros(len(us)) @@ -123,6 +156,177 @@ def _strip_vertical_offsets(x, y, z, psid, cell=STRIP_ALIGN_CELL, for k, p in enumerate(us) if abs(offsets[k]) >= threshold} +def _rolling_median(values, window): + """Médiane glissante (fenêtre tronquée aux bords) sur un petit vecteur.""" + n = len(values) + if window <= 1 or n == 0: + return np.array(values, dtype=np.float64) + half = window // 2 + return np.array([np.median(values[max(0, i - half):i + half + 1]) + for i in range(n)]) + + +def _strip_jitter_offsets(x, y, z, psid, t, cell=STRIP_ALIGN_CELL, + bin_seconds=STRIP_JITTER_BIN, + smooth=STRIP_JITTER_SMOOTH, + min_cells=STRIP_JITTER_MIN_CELLS, + max_corr=STRIP_JITTER_MAX, + threshold=STRIP_ALIGN_THRESHOLD): + """Mesure la gigue verticale intra-faisceau par fenêtres de temps GPS. + + Le calage constant retire un offset par faisceau, mais le décalage + vertical peut aussi varier au fil d'une même passe (vibration capteur, + bruit haute fréquence de la trajectoire). Chaque faisceau est découpé en + fenêtres de temps GPS ; l'offset robuste de chaque fenêtre est mesuré + contre la surface médiane des AUTRES faisceaux (maille 1 m, même + sélection de points bas que le calage constant — en recouvrement à deux, + une référence incluant le faisceau testé ne révélerait que la moitié du + décalage), puis la série est lissée (médiane glissante) pour ne pas + suivre le bruit de mesure. z doit déjà être corrigé des offsets constants. + + Args: + x, y, z, psid: coordonnées, Z (constant-calé) et PointSourceId des + points sol. + t: temps GPS (s) de chaque point. + cell: taille de maille de comparaison (m). + bin_seconds: durée d'une fenêtre de temps (s). + smooth: largeur de la médiane glissante (fenêtres). + min_cells: cellules partagées minimales pour valider une fenêtre. + max_corr: amplitude maximale d'une correction (garde-fou, m). + threshold: amplitude de série sous laquelle un faisceau est + considéré comme stable (pas de correction). + + Returns: + dict {psid: (times, corrections)} des corrections à SOUSTRAIRE, + interpolables linéairement au temps GPS de chaque point ; vide si + rien de mesurable (faisceau unique, temps absent, pas de + recouvrement). + """ + us, inv = np.unique(psid, return_inverse=True) + if len(us) < 2: + return {} + t = np.asarray(t, dtype=np.float64) + if t.size != z.size or not np.isfinite(t).all(): + return {} + built = _strip_surface_grid(x, y, z, inv, len(us), cell) + if built is None: + return {} + allcells, grid, key = built + covered = np.sum(~np.isnan(grid), axis=0) >= 2 + + # Référence d'un faisceau = médiane des autres faisceaux. + ref = np.full(grid.shape, np.nan) + if len(us) == 2: + ref[0], ref[1] = grid[1], grid[0] + else: + import warnings + with warnings.catch_warnings(): + warnings.simplefilter("ignore", RuntimeWarning) # tranches tout-NaN + for k in range(len(us)): + ref[k] = np.nanmedian(np.delete(grid, k, axis=0), axis=0) + + # Résidu vertical de chaque point contre la référence de son faisceau. + idx = np.minimum(np.searchsorted(allcells, key), len(allcells) - 1) + ref_point = ref[inv, idx] + hit = (allcells[idx] == key) & covered[idx] & np.isfinite(ref_point) + residual = np.full(np.asarray(z).shape, np.nan) + residual[hit] = np.asarray(z)[hit] - ref_point[hit] + + # Origine de temps propre à chaque faisceau : des passes d'une même tuile + # peuvent être espacées de plusieurs heures, les fenêtres restent ainsi + # dense autour du vol réel (et un saut de semaine GPS ne crée pas de + # géantes plages vides). + origin = np.full(len(us), np.inf) + np.minimum.at(origin, inv, t) + tb = np.floor((t - origin[inv]) / bin_seconds).astype(np.int64) + n_bins = int(tb.max()) + 1 + group = inv * n_bins + tb + + # Médiane robuste par groupe (faisceau, fenêtre) : les résidus finis sont + # triés en tête de segment, les NaN (sans référence) sont ignorés. + order = np.lexsort((np.where(np.isfinite(residual), residual, np.inf), group)) + g_s, r_s = group[order], residual[order] + starts = np.flatnonzero(np.r_[True, g_s[1:] != g_s[:-1]]) + ends = np.r_[starts[1:], len(g_s)] + + gids, times_c, med = [], [], [] + for s, e in zip(starts, ends): + finite = np.isfinite(r_s[s:e]) + if int(finite.sum()) < min_cells: + continue + gids.append(int(g_s[s])) + med.append(float(np.median(r_s[s:e][finite]))) + if not gids: + return {} + + gids = np.asarray(gids, dtype=np.int64) + med = np.asarray(med, dtype=np.float64) + result = {} + for k in range(len(us)): + selk = gids // n_bins == k + if int(selk.sum()) < 3: # série trop courte : correction non fiable + continue + b = gids[selk] % n_bins + cs = np.clip(_rolling_median(med[selk], smooth), -max_corr, max_corr) + if float(np.max(np.abs(cs))) < threshold: + continue + result[int(us[k])] = (origin[k] + (b + 0.5) * bin_seconds, cs) + return result + + +def _apply_strip_jitter(psid, t, jitter): + """Corrections de gigue interpolées au temps GPS de chaque point. + + Args: + psid, t: PointSourceId et temps GPS (s) de chaque point. + jitter: dict {psid: (times, corrections)} issu de _strip_jitter_offsets. + + Returns: + array des corrections à SOUSTRAIRE (0 pour les points sans série). + """ + corr = np.zeros(len(t)) + if not jitter: + return corr + psid = np.asarray(psid) + t = np.asarray(t, dtype=np.float64) + for p, (times, cs) in jitter.items(): + m = psid == p + if m.any(): + corr[m] = np.interp(t[m], times, cs) + return corr + + +def _strip_jitter_for_file(las_file, las, offsets): + """Gigue temporelle d'un LAS sol (offsets constants mesurés), mémoïsée.""" + try: + p = Path(las_file) + cache_key = (str(p), p.stat().st_mtime_ns) + except OSError: + cache_key = (str(las_file), 0) + if cache_key in _STRIP_JITTER_CACHE: + return _STRIP_JITTER_CACHE[cache_key] + jitter = {} + try: + t = np.asarray(las.gps_time, dtype=np.float64) + psid = np.asarray(las.point_source_id) + except AttributeError: + t, psid = None, None # dimensions absentes : pas de gigue mesurable + if (t is not None and psid is not None + and len(t) == len(las.points) == len(psid)): + z = np.asarray(las.z, dtype=np.float64) + if offsets: + lut = np.zeros(65536) + for p_, off in offsets.items(): + lut[int(p_) & 0xFFFF] = off + z = z - lut[np.asarray(las.point_source_id, dtype=np.int64)] + jitter = _strip_jitter_offsets( + np.asarray(las.x, dtype=np.float64), + np.asarray(las.y, dtype=np.float64), + z, psid, t) + _STRIP_JITTER_CACHE[cache_key] = jitter + return jitter + + def _strip_offsets_for_file(las_file, las): """Offsets de calage d'un LAS sol, mémoïsés par (chemin, mtime).""" try: @@ -148,18 +352,31 @@ def _strip_offsets_for_file(las_file, las): return offsets -def _write_strip_align_sidecar(dtm_dir, basename, output_suffix, offsets): - """Consigne les offsets de calage appliqués (version et seuil inclus). +def _write_strip_align_sidecar(dtm_dir, basename, output_suffix, offsets, + jitter=None): + """Consigne les calages appliqués (version, seuil, offsets et gigue). Le sidecar sert de suivi de cache : un DTM sans sidecar, ou produit avec - une version/un seuil différents, est régénéré. Il est écrit même quand - aucun offset n'a été appliqué, pour ne pas re-mesurer une tuile déjà - connue comme bien alignée. + une version/un seuil/des paramètres de gigue différents, est régénéré. Il + est écrit même quand aucune correction n'a été appliquée, pour ne pas + re-mesurer une tuile déjà connue comme bien alignée. """ + jitter_payload = {} + for p, (times, cs) in (jitter or {}).items(): + cs = np.asarray(cs, dtype=np.float64) + jitter_payload[str(p)] = { + "bins": int(len(times)), + "rms_m": round(float(np.sqrt(np.mean(cs ** 2))), 4), + "max_m": round(float(np.max(np.abs(cs))), 4), + "series_m": [round(float(c), 4) for c in cs], + } payload = { "version": STRIP_ALIGN_VERSION, "threshold": STRIP_ALIGN_THRESHOLD, "offsets": offsets, + "jitter_bin": STRIP_JITTER_BIN, + "jitter_smooth": STRIP_JITTER_SMOOTH, + "jitter": jitter_payload, } try: sidecar = Path(dtm_dir) / f"{basename}_dtm{output_suffix}_stripalign.json" @@ -813,73 +1030,139 @@ def _interpolate_holes(dtm, downsample=8): return filled, int(holes.sum()) -def _min_return_grid(laz_file, width, height, bounds, chunk_size=2_000_000, - strip_offsets=None): - """Rasterize the per-cell minimum z (lowest return) of the full point cloud. +# ------------------------------------------------------------ +# Raccord des bords (edge buffer avec les tuiles adjacentes) +# ------------------------------------------------------------ +# Les visualisations à grand noyau (openness/SVF : rayons jusqu'à 100 m, +# LRM : 15 m) tronquent leur fenêtre au bord de la dalle : les rendus +# présentent alors une bande d'artefacts à chaque changement de tuile. Avec +# edge_buffer > 0, le MNT est rastérisé sur la dalle nominale 1 km ÉTENDUE +# d'une bande de `edge_buffer` mètres remplie avec les points sol des 8 +# tuiles LAZ voisines (lecture PDAL en flux : découpe + filtre de classes). +# Les visualisations voient ainsi le terrain réel au-delà du bord, et les +# images finales sont recadrées sur la dalle exacte (cf. rendering.py). +EDGE_BUFFER_TAG = "LIDAR_EDGE_BUFFER" # tag GeoTIFF : tampon utilisé (m) - 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. +# 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)] + + +def _tile_coords(name): + """Coordonnées (col, row) en km d'un nom de fichier LHD, ou None.""" + from .index import parse_basename_coords + return parse_basename_coords(Path(name).name) + + +def _neighbor_laz_files(source_laz): + """Liste les 8 fichiers LAZ/LAS adjacents à `source_laz` dans son dossier.""" + coords = _tile_coords(source_laz) + if coords is None: + return [] + col, row = coords + directory = Path(source_laz).parent + neighbors = [] + for dcol, drow in _NEIGHBOR_OFFSETS: + nc, nr = col + dcol, row + drow + # Noms LHD : col/row en km sur 4 chiffres complétés (ex. 0638_6628) ; + # 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")) + if matches: + neighbors.append(sorted(matches)[0]) + else: + logger.debug(f" Voisine absente : LHD_FXX_{nc:04d}_{nr:04d} (bande de bord vide)") + return neighbors + + +def _neighbor_ground_points(source_laz, bounds, classes): + """Points sol des tuiles voisines dans `bounds` (raccord des bords). + + Lecture PDAL en flux par voisine : découpe sur l'emprise étendue puis + filtre de classes (mêmes codes que le MNT). Best-effort : une voisine + illisible ou absente est ignorée — la bande correspondante reste vide. 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. - strip_offsets: Optionnel : dict {point_source_id: offset} issu du - calage des faisceaux, retranché aux Z du nuage complet pour - rester cohérent avec le MNT calé. + source_laz: LAZ de la tuile traitée (repère pour trouver les voisines). + bounds: (min_x, min_y, max_x, max_y) de l'emprise étendue. + classes: codes LAS à extraire (ex. [2] = sol). Returns: - (height, width) float32 array of per-cell min z (NaN where no point). + (xs, ys, zs) concaténés (tableaux vides si aucune voisine). """ - import laspy + import tempfile + + neighbors = _neighbor_laz_files(source_laz) + if not neighbors: + return (np.empty(0),) * 3 + 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]] + codes = sorted(set(int(c) for c in classes)) or [2] + limits = ",".join(f"Classification[{c}:{c}]" for c in codes) - lut = None - if strip_offsets: - lut = np.zeros(65536) - for p, off in strip_offsets.items(): - lut[int(p) & 0xFFFF] = off + xs, ys, zs = [], [], [] + found = 0 + for neighbor in neighbors: + tmp_path = None + try: + with tempfile.NamedTemporaryFile(suffix='.las', delete=False) as tmp: + tmp_path = tmp.name + pipeline = json.dumps({ + "pipeline": [ + {"type": "readers.las", "filename": str(neighbor)}, + {"type": "filters.crop", + "bounds": f"([{min_x},{max_x}],[{min_y},{max_y}])"}, + {"type": "filters.range", "limits": limits}, + {"type": "writers.las", "filename": tmp_path}, + ] + }) + result = subprocess.run( + ["pdal", "pipeline", "--stdin"], + input=pipeline, capture_output=True, text=True, timeout=300 + ) + if result.returncode != 0: + raise RuntimeError(result.stderr[:200]) + import laspy + las = laspy.read(tmp_path) + if len(las.points) > 0: + xs.append(np.asarray(las.x, dtype=np.float64)) + ys.append(np.asarray(las.y, dtype=np.float64)) + zs.append(np.asarray(las.z, dtype=np.float64)) + found += 1 + except Exception as e: + logger.debug(f" Voisine {Path(neighbor).name} ignorée : {e}") + finally: + if tmp_path: + try: + Path(tmp_path).unlink(missing_ok=True) + except Exception: + pass - 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) - if lut is not None: - try: - z = z - lut[np.asarray(points.point_source_id, dtype=np.int64)] - except AttributeError: - pass - 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) + if not found: + logger.warning(" Aucune voisine lisible — bande de bord vide") + return (np.empty(0),) * 3 + logger.info(f" Raccord bords : {found} voisine(s), " + f"{sum(len(a) for a in xs):,} pts sol") + return np.concatenate(xs), np.concatenate(ys), np.concatenate(zs) + +def read_dtm_edge_buffer(dtm_path): + """Tampon de raccord enregistré dans un DTM (m ; 0 si absent/illisible).""" 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 + with rasterio.open(dtm_path) as src: + tags = src.tags() + return float(tags.get(EDGE_BUFFER_TAG, 0.0) or 0.0) + except Exception: + return 0.0 def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, - output_suffix="", source_laz=None, bare_earth=False, - pure=False, strip_align=True): + output_suffix="", source_laz=None, + pure=False, strip_align=True, edge_buffer=0.0, + neighbor_classes=None): """Create DTM using fast binning method with gap filling. Args: @@ -889,19 +1172,26 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, 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. + source_laz: Optionnel : LAZ complet de la tuile traitée (repère pour + lire les 8 tuiles voisines du raccord de bords, cf. edge_buffer). 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). strip_align: Si True (défaut), mesure et corrige les écarts verticaux - entre faisceaux de vol (PointSourceId) avant rastérisation ; les - offsets ≥ STRIP_ALIGN_THRESHOLD (0,5 cm) sont consignés dans un - sidecar *_dtm*_stripalign.json. + entre faisceaux de vol (PointSourceId) avant rastérisation : offsets + constants ≥ STRIP_ALIGN_THRESHOLD (0,5 cm) puis gigue intra-faisceau + par fenêtres de temps GPS ; tout est consigné dans un sidecar + *_dtm*_stripalign.json. + edge_buffer: Raccord des bords en mètres (0 = désactivé). Le MNT couvre + alors la dalle nominale 1 km étendue de cette bande, remplie avec + les points sol des 8 tuiles LAZ voisines (source_laz requis) ; les + visualisations calculent sur l'emprise étendue puis les images sont + recadrées sur la dalle exacte (rendering.py). Le tampon est inscrit + dans le tag GeoTIFF LIDAR_EDGE_BUFFER pour l'invalidation du cache. + neighbor_classes: codes LAS extraits chez les voisines (défaut : les + classes IGN du MNT, ex. [2]). Les voisines sont lues dans leur + pré-classification fournisseur, même si la tuile centrale est + classée SMRF/CSF (bande de contexte, quelques cm d'écart au pire). Returns: Path to output DTM GeoTIFF, or None on failure. @@ -934,6 +1224,8 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, # Calage vertical des faisceaux de vol avant rastérisation : best-effort, # en cas d'échec de la mesure on continue non calé (jamais d'abort). strip_offsets = {} + strip_jitter = {} + gps_time = None if strip_align: try: strip_offsets = _strip_offsets_for_file(las_file, las) @@ -946,6 +1238,24 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, else: logger.debug(" Calage faisceaux : aucun écart >= " f"{STRIP_ALIGN_THRESHOLD * 100:.1f} cm, rien à corriger") + # 2ᵉ passe : gigue intra-faisceau (temps GPS requis, ignorée sinon). + try: + gps_time = np.asarray(las.gps_time, dtype=np.float64) + except AttributeError: + gps_time = None + if gps_time is not None and len(gps_time) == len(las.points): + try: + strip_jitter = _strip_jitter_for_file(las_file, las, strip_offsets) + except Exception as e: + logger.warning(f" Mesure de la gigue intra-faisceau impossible ({e}) — gigue non corrigée") + strip_jitter = {} + if strip_jitter: + for p_, (jt, jc) in sorted(strip_jitter.items()): + logger.info(f" Gigue PSID {p_} : ±{np.max(np.abs(jc)) * 100:.1f} cm " + f"(rms {np.sqrt(np.mean(np.asarray(jc) ** 2)) * 100:.1f} cm, " + f"{len(jt)} fenêtres de {STRIP_JITTER_BIN:g} s)") + else: + logger.debug(" Gigue intra-faisceau : rien à corriger") try: @@ -955,6 +1265,32 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, width = int(np.ceil((max_x - min_x) / resolution)) height = int(np.ceil((max_y - min_y) / resolution)) + # Raccord des bords : emprise = dalle nominale 1 km (alignée sur la + # grille multi-tuiles) + bande de edge_buffer mètres remplie par les + # points sol des voisines. Sinon : bornes de l'en-tête (historique). + used_edge_buffer = 0.0 + if edge_buffer > 0: + coords = _tile_coords(source_laz or las_file) + if coords is not None: + buffer_px = max(1, int(round(edge_buffer / resolution))) + buffer_m = buffer_px * resolution + col_km, row_km = coords + # Grille LHD : (col, row) = coin nord-ouest en km → + # X ∈ [col, col+1] km, Y ∈ [row-1, row] km (bord nord = row). + min_x = float(col_km) * 1000.0 + max_x = min_x + 1000.0 + max_y = float(row_km) * 1000.0 + min_y = max_y - 1000.0 + ext_bounds = (min_x - buffer_m, min_y - buffer_m, + max_x + buffer_m, max_y + buffer_m) + width = int(round(1000.0 / resolution)) + 2 * buffer_px + height = width + min_x, min_y, max_x, max_y = ext_bounds + used_edge_buffer = float(edge_buffer) + else: + logger.warning(" Raccord des bords impossible : nom de " + f"fichier non LHD ({basename}) — tuile seule") + logger.debug(f" Bounds: X[{min_x:.1f}, {max_x:.1f}] Y[{min_y:.1f}, {max_y:.1f}]") logger.debug(f" Grid: {width}x{height} pixels ({len(las.points):,} points)") logger.info(f" Rasterisation {width}x{height} ({len(las.points):,} points)...") @@ -967,6 +1303,20 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, for p, off in strip_offsets.items(): lut[int(p) & 0xFFFF] = off zs = zs - lut[np.asarray(las.point_source_id, dtype=np.int64)] + if strip_jitter and gps_time is not None: + zs = zs - _apply_strip_jitter(las.point_source_id, gps_time, strip_jitter) + + # Points sol des tuiles voisines dans la bande de raccord (best-effort, + # non calés par faisceau : bande de contexte, l'image finale est + # recadrée sur la dalle avant livraison). + if used_edge_buffer > 0 and source_laz is not None: + nx, ny, nz = _neighbor_ground_points( + source_laz, (min_x, min_y, max_x, max_y), + neighbor_classes if neighbor_classes is not None else [2]) + if len(nx): + xs = np.concatenate([xs, nx]) + ys = np.concatenate([ys, ny]) + zs = np.concatenate([zs, nz]) stat = binned_statistic_2d( xs, ys, zs, @@ -981,20 +1331,9 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, # 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), - strip_offsets=strip_offsets) - # Cellules sans sol mesuré (NaN) : le retour le plus bas devient - # la mesure — sinon la comparaison NaN est fausse et la cellule - # retombe sur l'interpolation en fin de passe. - has_min = ~np.isnan(min_grid) - lower = has_min & (min_grid < dtm) - lower |= has_min & np.isnan(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") + # Volontairement PAS de plancher au retour le plus bas : sous canopée + # dense ce retour est la végétation, qui imprimerait les arbres dans + # le MNT. # Fill small gaps (< 1 m from data) precisely — comme avant nan_count = np.count_nonzero(np.isnan(dtm)) @@ -1026,10 +1365,14 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, compress='lzw' ) as dst: dst.write(dtm.astype('float32'), 1) + if used_edge_buffer > 0: + # Tampon de raccord inscrit dans le fichier : changement de + # --edge-buffer ⇒ invalidation automatique du cache DTM. + dst.update_tags(**{EDGE_BUFFER_TAG: f"{used_edge_buffer:g}"}) if strip_align: _write_strip_align_sidecar(dtm_dir, basename, output_suffix, - strip_offsets) + strip_offsets, strip_jitter) logger.info(f" ✓ DTM créé: {output_tif.name}") return output_tif diff --git a/lidar_pipeline/index.py b/lidar_pipeline/index.py index ddd75a7..1d07fb6 100644 --- a/lidar_pipeline/index.py +++ b/lidar_pipeline/index.py @@ -1352,8 +1352,8 @@ _HTML_TEMPLATE = """ -