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 = """
-