"""DTM generation from classified LiDAR point clouds. Handles ground classification via PDAL (IGN supplier pre-classification, SMRF or CSF) and DTM rasterisation using scipy binned_statistic_2d (GPU-first via gpu.bin_mean_2d, scipy as fallback). Gaps without LiDAR data (common in complex/rocky terrain) are filled with a terrain-aware interpolation so the DTM stays continuous. """ import json import logging import subprocess import time from pathlib import Path import numpy as np import rasterio from rasterio.transform import from_bounds from scipy.stats import binned_statistic_2d from .gpu import bin_mean_2d logger = logging.getLogger("lidar") # Classes LAS exploitables de la pré-classification LiDAR HD (noms → codes) IGN_CLASS_NAMES = { "sol": 2, "unclassified": 1, "non-classe": 1, "eau": 9, "virtuel": 66, "pont": 17, "sursol": 64, } # Calage vertical des faisceaux de vol (strip alignment). Une tuile LiDAR HD # est couverte par plusieurs passes d'acquisition, chacune portée par un ou # deux PointSourceId (faisceaux). Les passes sont bien alignées horizontalement # mais certains faisceaux portent un biais vertical de quelques cm (mesuré # jusqu'à ~5 cm sur LHD_FXX_1000_6882 : les deux faisceaux d'une passe écartés # de ±2,5 cm). À 0,2 m/px ces écarts créent des marches et du bruit aux # coutures des zones de recouvrement. On mesure l'offset robuste de chaque # faisceau sur les points sol de la tuile elle-même (les PointSourceId changent # 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. # # 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). # # 3ᵉ passe — décalage ligne à ligne : les fenêtres de 0,1 s (lissées sur # 0,5 s) regroupent ~15 lignes de balayage (~150 lignes/s) et ne travaillent # qu'en recouvrement. Or deux lignes SUCCESSIVES d'une même passe peuvent # différer de 1 à 2 cm (motif alterné, mesuré sur LHD_FXX_0999_6882) : stries # fines perpendiculaires au vol sur tout le MNT. Chaque faisceau est découpé # en lignes (sauts de scan_angle en dents de scie), le décalage robuste de # chaque ligne (décalage ET inclinaison le long de la ligne : le roulis # bascule les lignes) est mesuré contre la surface de SON faisceau (maille 0,5 m, # boîte 1,5 m, plan local ajusté à la position réelle du point, 3 itérations # pour que la ligne ne fausse pas sa propre référence), et seule la composante # ligne à ligne est retirée (série − lissage gaussien σ 3 lignes ; une médiane # glissante suivrait un motif alterné au lieu de l'effacer) : les variations # lentes restent aux passes précédentes. Fonctionne sans recouvrement ; requiert # gps_time et scan_angle. # # Ajustement conjoint (prioritaire) : la surface d'un faisceau absorbe toute # erreur plus large que sa boîte de référence ; là où plusieurs faisceaux se # recouvrent, chaque ligne (décalage + inclinaison) est donc recalée contre le # consensus des AUTRES faisceaux, à toutes les échelles (roulis lent d'une # passe entière, lignes isolées très décalées), par pas amortis et itérés. # La correction sur son propre faisceau ne sert plus qu'aux lignes sans # recouvrement. STRIP_ALIGN_VERSION = 3 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) STRIP_LINE_CELL = 0.5 # m : maille de la surface de référence d'un faisceau STRIP_LINE_BOX = 3 # mailles : lissage de la référence (boîte 1,5 m) STRIP_LINE_WINDOW = 3 # lignes : σ du lissage gaussien retiré (garde le ligne-à-ligne) STRIP_LINE_ITERS = 3 # itérations (atténuation par la ligne elle-même < 1 %) STRIP_LINE_MAX = 0.05 # m : correction maxi d'une ligne (garde-fou) STRIP_LINE_MIN_POINTS = 100 # points sol mini pour mesurer une ligne STRIP_LINE_GAP = 0.05 # s : trou de temps qui coupe une ligne (fin de passe) STRIP_LINE_MIN_RMS = 0.002 # m : faisceau laissé tel quel sous ce niveau STRIP_LINE_MODEL = "conjoint+decalage+inclinaison+profil-angle" # modèle (consigné, invalide le cache) _SCAN_ANGLE_UNIT = 0.006 # ° par unité de scan_angle (LAS 1.4) STRIP_ANGLE_BIN = 1.0 # ° : pas du profil de correction par faisceau selon l'angle STRIP_ANGLE_MIN_POINTS = 500 # points en recouvrement mini par classe d'angle (sinon valeur voisine) STRIP_LINE_SUBSAMPLE = 3 # estimation sur 1 point sur 3 (correction appliquée à tous) STRIP_JOINT_CELL = 1.0 # m : maille du consensus des autres faisceaux STRIP_JOINT_ITERS = 8 # itérations maxi de l'ajustement conjoint STRIP_JOINT_DAMPING = 0.5 # pas amorti : deux faisceaux se rapprochent sans se croiser STRIP_JOINT_TOL = 0.001 # m : arrêt quand le pas moyen passe sous 1 mm STRIP_JOINT_MIN_POINTS = 30 # points en recouvrement mini pour recaler une ligne STRIP_JOINT_MAX = 0.15 # m : correction maxi d'un point (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 = {} _STRIP_LINES_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, threshold=STRIP_ALIGN_THRESHOLD, min_shared=STRIP_ALIGN_MIN_SHARED): """Mesure les offsets verticaux relatifs entre faisceaux d'une tuile. Méthode : surface sol par faisceau (moyenne des points à moins de 0,5 m du minimum de chaque maille de 1 m, robuste à la végétation résiduelle), puis offset de chaque faisceau = médiane de l'écart à la surface médiane de référence (itéré 3 fois), sur les seules cellules couvertes par au moins deux faisceaux. Ne corrige rien d'autre que le vertical : le calage horizontal mesuré sur les données est excellent (≤ 1 cm). Args: x, y, z, psid: coordonnées et PointSourceId des points sol. cell: taille de maille de comparaison (m). threshold: seuil (m) en dessous duquel un offset est ignoré. min_shared: cellules communes minimales pour valider un faisceau. Returns: dict {psid: offset} des offsets à SOUSTRAIRE (z - offset), ne contenant que les |offset| >= threshold ; vide si rien à corriger (faisceau unique, pas de recouvrement, tuile déjà alignée). """ us, inv = np.unique(psid, return_inverse=True) if len(us) < 2: return {} built = _strip_surface_grid(x, y, z, inv, len(us), cell) if built is None: return {} allcells, grid, _key = built comparable = np.sum(~np.isnan(grid), axis=0) >= 2 offsets = np.zeros(len(us)) for _ in range(3): ref = np.nanmedian(grid + offsets[:, None], axis=0) for k in range(len(us)): m = ~np.isnan(grid[k]) & comparable if int(m.sum()) >= min_shared: offsets[k] = np.median(grid[k][m] - ref[m]) return {int(p): round(float(offsets[k]), 3) 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 _module_of(arr): """Module (numpy ou cupy) d'un tableau.""" mod = type(arr).__module__ if mod.startswith("cupy"): import cupy return cupy return np def _array_module(use_gpu): """cupy si le GPU est actif et demandé, sinon numpy.""" if use_gpu: from . import gpu as _gpu if _gpu.is_gpu_active() and _gpu._cp is not None: return _gpu._cp return np def _scan_line_ids(t, angle, gap=STRIP_LINE_GAP): """Identifiant de ligne de balayage de points d'UN faisceau triés par temps. Nouvelle ligne à chaque saut de scan_angle de plus de la moitié de l'amplitude dans le sens opposé au balayage (retour de la dent de scie) ou à chaque trou de temps > gap (fin de passe). Les trous laissés par la végétation retirée dans une ligne ne la coupent pas. """ m = _module_of(angle) angle = m.asarray(angle, dtype=m.float64) if len(angle) == 0: return m.zeros(0, dtype=m.int64) amp = float(m.percentile(angle, 99) - m.percentile(angle, 1)) d = m.diff(angle) moving = d[d != 0] sweep = float(m.sign(m.median(moving))) if len(moving) else 1.0 # Retour de la dent de scie : grand saut de sens OPPOSÉ au balayage. Un # trou de végétation fait aussi sauter l'angle, mais dans le sens du # balayage : il ne coupe pas la ligne. new = m.concatenate([m.ones(1, dtype=bool), (d * sweep < -0.5 * amp) | (m.diff(t) > gap)]) return m.cumsum(new) - 1 def _group_median(values, groups, n_groups, min_count): """Médiane des valeurs finies par groupe (vectorisée), NaN sous min_count.""" finite = np.isfinite(values) order = np.lexsort((np.where(finite, values, np.inf), groups)) g = groups[order] v = values[order] nf = np.bincount(groups[finite], minlength=n_groups) start = np.searchsorted(g, np.arange(n_groups)) med = np.full(n_groups, np.nan) ok = nf >= max(1, min_count) lo = start[ok] + (nf[ok] - 1) // 2 hi = start[ok] + nf[ok] // 2 med[ok] = 0.5 * (v[lo] + v[hi]) return med def _line_fit(r, w, line, u, n_lines, min_points): """Moindres carrés pondérés par ligne : r ≈ a + b·u (sommes par bincount). Returns: (a, b, ok) ; b vaut 0 là où l'étendue de u ne permet pas d'estimer une inclinaison (a seul, moyenne pondérée). """ m = _module_of(r) rw = m.where(w, r, 0.0) wf = w.astype(m.float64) s0 = m.bincount(line, weights=wf, minlength=n_lines) s1 = m.bincount(line, weights=wf * u, minlength=n_lines) s2 = m.bincount(line, weights=wf * u * u, minlength=n_lines) sr = m.bincount(line, weights=rw, minlength=n_lines) sur = m.bincount(line, weights=rw * u, minlength=n_lines) det = s0 * s2 - s1 * s1 ok = s0 >= min_points tilt = ok & (det > 1e-2 * m.maximum(s0, 1) ** 2) a = m.where(ok, sr / m.maximum(s0, 1), 0.0) b = m.zeros(n_lines) safe = m.where(tilt, det, 1.0) a = m.where(tilt, (s2 * sr - s1 * sur) / safe, a) b = m.where(tilt, (s0 * sur - s1 * sr) / safe, b) return a, b, ok def _robust_mask(r, floor=0.03, k=5.0): """Résidus finis sous k MAD (plancher floor) : écarte végétation basse, points mal classés et bords de trous sans trier chaque ligne.""" m = _module_of(r) f = m.isfinite(r) if not bool(f.any()): return f mad = 1.4826 * float(m.median(m.abs(r[f]))) return f & (m.abs(r) < max(floor, k * mad)) def _scan_line_corrections_beam(x, y, z, t, angle, cell=STRIP_LINE_CELL, box=STRIP_LINE_BOX, window=STRIP_LINE_WINDOW, iters=STRIP_LINE_ITERS, max_corr=STRIP_LINE_MAX, min_points=STRIP_LINE_MIN_POINTS, subsample=STRIP_LINE_SUBSAMPLE): """Corrections ligne à ligne d'un faisceau contre SA surface (à SOUSTRAIRE). Chaque ligne est modélisée par un décalage ET une inclinaison le long de la ligne (a + b·u, u = scan_angle normalisé) : une erreur de roulis bascule la ligne, un bout plus haut que l'autre. Seule la composante ligne à ligne est retirée (série − gaussienne σ window lignes). Estimation sur 1 point sur subsample, moindres carrés tronqués (5 MAD) vectorisés. Returns: (corr par point dans l'ordre d'entrée, décalages par ligne, inclinaisons par ligne au bord de fauchée). """ from scipy.ndimage import gaussian_filter1d, uniform_filter n = len(z) o = np.argsort(t, kind="stable") line_all = np.empty(n, dtype=np.int64) line_all[o] = _scan_line_ids(t[o], angle[o]) nl = int(line_all.max()) + 1 if n else 0 tot_a, tot_b = np.zeros(nl), np.zeros(nl) if nl < 10 * window: return np.zeros(n), tot_a, tot_b u_all = np.asarray(angle, dtype=np.float64) / max(np.percentile(np.abs(angle), 99), 1e-9) sel = np.arange(n) % max(1, int(subsample)) == 0 xs, ys, zs = x[sel], y[sel], np.asarray(z, dtype=np.float64)[sel] line, u = line_all[sel], u_all[sel] # Référence = plan local : moyennes (x, y, z) par boîte de mailles, pente # de la surface lissée ; évalué à la position réelle du point (la moyenne # d'une maille n'est pas en son centre : sur une pente, une interpolation # au centre crée un biais qui dépend de la position de la ligne). x0, y0 = xs.min(), ys.min() ix = np.floor((xs - x0) / cell).astype(np.int64) iy = np.floor((ys - y0) / cell).astype(np.int64) W, H = int(ix.max()) + 1, int(iy.max()) + 1 flat = iy * W + ix count_b = uniform_filter(np.bincount(flat, minlength=W * H).reshape(H, W).astype(np.float64), box) valid = count_b > 0 inv_c = np.where(valid, 1.0 / np.maximum(count_b, 1e-12), np.nan) mean_x = uniform_filter(np.bincount(flat, weights=xs - x0, minlength=W * H).reshape(H, W), box) * inv_c mean_y = uniform_filter(np.bincount(flat, weights=ys - y0, minlength=W * H).reshape(H, W), box) * inv_c dxp = (xs - x0) - mean_x.ravel()[flat] dyp = (ys - y0) - mean_y.ravel()[flat] min_pts = max(10, min_points // max(1, int(subsample))) z_work = zs.copy() for _ in range(iters): mean_z = uniform_filter(np.bincount(flat, weights=z_work, minlength=W * H).reshape(H, W), box) * inv_c gz_y, gz_x = np.gradient(np.where(valid, mean_z, np.nanmean(mean_z)), cell) r = z_work - (mean_z.ravel()[flat] + gz_x.ravel()[flat] * dxp + gz_y.ravel()[flat] * dyp) a, b, ok = _line_fit(r, _robust_mask(r), line, u, nl, min_pts) if ok.sum() < 10 * window: break idx = np.flatnonzero(ok) a = np.interp(np.arange(nl), idx, a[ok]) b = np.interp(np.arange(nl), idx, b[ok]) # Seule la composante ligne à ligne est retirée (une médiane glissante # suivrait un motif alterné au lieu de l'effacer) da = a - gaussian_filter1d(a, window, mode="nearest") db = b - gaussian_filter1d(b, window, mode="nearest") da[~ok] = 0.0 db[~ok] = 0.0 tot_a += da tot_b += db z_work = zs - np.clip(tot_a[line] + tot_b[line] * u, -max_corr, max_corr) corr = np.clip(tot_a[line_all] + tot_b[line_all] * u_all, -max_corr, max_corr) return corr, tot_a, tot_b def _joint_line_corrections(x, y, z, psid, t, angle, cell=STRIP_JOINT_CELL, iters=STRIP_JOINT_ITERS, damping=STRIP_JOINT_DAMPING, tol=STRIP_JOINT_TOL, min_points=STRIP_JOINT_MIN_POINTS, subsample=STRIP_LINE_SUBSAMPLE, max_corr=STRIP_JOINT_MAX, use_gpu=True): """Ajustement conjoint des lignes de tous les faisceaux contre le consensus des AUTRES faisceaux (décalage + inclinaison par ligne). À chaque itération, le résidu d'un point est mesuré contre la moyenne des autres faisceaux de sa maille (ramenée à sa position par la pente de la surface), chaque ligne est ajustée par moindres carrés tronqués, et une fraction damping du pas est appliquée à toutes les lignes à la fois ; l'altitude moyenne est recentrée (pas de dérive d'ensemble). Sur GPU (CuPy) si disponible, repli numpy sur toute erreur. Returns: (corr à SOUSTRAIRE par point, id de ligne global par point, lignes recalées (bool par ligne), nombre d'itérations) — numpy. """ m = _array_module(use_gpu) if m is not np: try: return _joint_line_corrections_impl(m, x, y, z, psid, t, angle, cell, iters, damping, tol, min_points, subsample, max_corr) except Exception as e: logger.warning(f" Ajustement conjoint GPU impossible ({e}) — repli CPU") return _joint_line_corrections_impl(np, x, y, z, psid, t, angle, cell, iters, damping, tol, min_points, subsample, max_corr) def _joint_line_corrections_impl(m, x, y, z, psid, t, angle, cell, iters, damping, tol, min_points, subsample, max_corr): def host(a): return a.get() if m is not np else a n = len(z) psid_h = np.asarray(psid) beams, bidx_h = np.unique(psid_h, return_inverse=True) # Classe d'angle par faisceau : profil d'étalonnage commun à toutes les # lignes d'un faisceau (écart non linéaire en travers de la fauchée). abin_h = np.round(np.asarray(angle, dtype=np.float64) * _SCAN_ANGLE_UNIT / STRIP_ANGLE_BIN).astype(np.int64) amin = int(abin_h.min()) if n else 0 n_ab = int(abin_h.max()) - amin + 1 if n else 1 abin_h = bidx_h * n_ab + (abin_h - amin) x, y = m.asarray(x, dtype=m.float64), m.asarray(y, dtype=m.float64) z, t = m.asarray(z, dtype=m.float64), m.asarray(t, dtype=m.float64) angle = m.asarray(angle, dtype=m.float64) bidx = m.asarray(bidx_h) gl = m.zeros(n, dtype=m.int64) u = m.zeros(n) base = 0 for i in range(len(beams)): idx = m.flatnonzero(bidx == i) o = m.argsort(t[idx]) gl[idx[o]] = _scan_line_ids(t[idx][o], angle[idx][o]) + base base = int(gl[idx].max()) + 1 u[idx] = angle[idx] / max(float(m.percentile(m.abs(angle[idx]), 99)), 1e-9) n_lines = base touched = m.zeros(n_lines, dtype=bool) if len(beams) < 2 or n == 0: return np.zeros(n), host(gl), host(touched), 0 abin = m.asarray(abin_h) n_bins = len(beams) * n_ab sel = m.arange(0, n, max(1, int(subsample))) xs, ys, zs = x[sel], y[sel], z[sel] gs, us, bs, cs = gl[sel], u[sel], bidx[sel], abin[sel] x0, y0 = float(xs.min()), float(ys.min()) ix = m.floor((xs - x0) / cell).astype(m.int64) iy = m.floor((ys - y0) / cell).astype(m.int64) W, H = int(ix.max()) + 1, int(iy.max()) + 1 nk = W * H key = iy * W + ix bkey = bs * nk + key dxc = (xs - x0) - (ix + 0.5) * cell dyc = (ys - y0) - (iy + 0.5) * cell n_all = m.bincount(key, minlength=nk) n_own = m.bincount(bkey, minlength=len(beams) * nk) n_other = n_all[key] - n_own[bkey] overlap = n_other >= 2 min_pts = max(10, min_points // max(1, int(subsample))) A = m.zeros(n_lines) B = m.zeros(n_lines) C = m.zeros(n_bins) min_bin = max(30, STRIP_ANGLE_MIN_POINTS // max(1, int(subsample))) it = 0 for it in range(1, iters + 1): zc = zs - m.clip(A[gs] + B[gs] * us + C[cs], -max_corr, max_corr) s_all = m.bincount(key, weights=zc, minlength=nk) s_own = m.bincount(bkey, weights=zc, minlength=len(beams) * nk) mean = m.where(n_all > 0, s_all / m.maximum(n_all, 1), m.nan).reshape(H, W) gy, gx = m.gradient(m.where(m.isfinite(mean), mean, m.nanmean(mean)), cell) ref = ((s_all[key] - s_own[bkey]) / m.maximum(n_other, 1) + gx.ravel()[key] * dxc + gy.ravel()[key] * dyc) r = m.where(overlap, zc - ref, m.nan) w = m.zeros(len(r), dtype=bool) for i in range(len(beams)): mb = bs == i w[mb] = _robust_mask(r[mb]) a, b, ok = _line_fit(r, w, gs, us, n_lines, min_pts) touched |= ok # Profil par faisceau et classe d'angle : moyenne tronquée du résidu # restant après le pas de ligne ; sa moyenne (portée par les lignes) # est retirée faisceau par faisceau. Sa pente est conservée : quand la # fauchée ne traverse la dalle qu'en partie, l'inclinaison des lignes # est indéterminée et seul le profil peut la porter. r2 = m.where(w, r - (a[gs] + b[gs] * us), 0.0) wf = w.astype(m.float64) cnt = m.bincount(cs, weights=wf, minlength=n_bins) prof = m.where(cnt >= min_bin, m.bincount(cs, weights=r2, minlength=n_bins) / m.maximum(cnt, 1), 0.0) su = m.bincount(cs, weights=wf * us, minlength=n_bins) ub = m.where(cnt > 0, su / m.maximum(cnt, 1), 0.0) for i in range(len(beams)): sl = slice(i * n_ab, (i + 1) * n_ab) k = cnt[sl] >= min_bin if int(k.sum()) >= 3: wk = cnt[sl][k] uk, pk = ub[sl][k], prof[sl][k] um, pm = float((wk * uk).sum() / wk.sum()), float((wk * pk).sum() / wk.sum()) fitted = prof[sl] - pm # Classe trop pauvre (bord de fauchée, classe incomplète) : # valeur de la classe valide la plus proche, pas zéro — c'est # précisément au bord que l'écart est le plus fort. idx = m.arange(n_ab, dtype=m.float64) prof[sl] = m.interp(idx, idx[k], fitted[k]) else: prof[sl] = 0.0 A += damping * a B += damping * b C += damping * prof A -= float(m.mean(A[gs] + B[gs] * us + C[cs])) # pas de dérive d'ensemble step_lines = float(m.mean(m.abs(a[ok]))) * damping if bool(ok.any()) else 0.0 step_prof = float(m.max(m.abs(prof))) * damping if max(step_lines, step_prof) < tol: break corr = m.clip(A[gl] + B[gl] * u + C[abin], -max_corr, max_corr) return host(corr), host(gl), host(touched), it def _rolling_median_fast(values, window): """Médiane glissante centrée (bords répétés), vectorisée.""" from numpy.lib.stride_tricks import sliding_window_view half = window // 2 return np.median(sliding_window_view(np.pad(values, half, mode="edge"), window), axis=1) def _scan_line_corrections(x, y, z, psid, t, angle, min_rms=STRIP_LINE_MIN_RMS): """Corrections ligne à ligne de tous les faisceaux d'une tuile. 1. Ajustement conjoint contre les autres faisceaux (lignes en recouvrement). 2. Lignes sans recouvrement : composante ligne à ligne contre la surface de leur propre faisceau. z doit être déjà calé (offsets constants et gigue temporelle). Returns: (corrections à SOUSTRAIRE par point, {psid: (lignes, rms, max)} des faisceaux corrigés). """ psid = np.asarray(psid) z = np.asarray(z, dtype=np.float64) corr, gl, touched, _ = _joint_line_corrections(x, y, z, psid, t, angle) stats = {} for p in np.unique(psid): m = psid == p if int(m.sum()) < 50 * STRIP_LINE_MIN_POINTS: continue alone = ~touched[gl[m]] if alone.mean() > 0.05: c_self, per_line, _ = _scan_line_corrections_beam( x[m], y[m], z[m] - corr[m], t[m], angle[m]) c_beam = corr[m] + np.where(alone, c_self, 0.0) else: c_beam = corr[m] rms = float(np.sqrt(np.mean(c_beam ** 2))) if m.any() else 0.0 if rms < min_rms: corr[m] = 0.0 continue corr[m] = c_beam n_lines = int(len(np.unique(gl[m]))) stats[int(p)] = (n_lines, rms, float(np.max(np.abs(c_beam)))) return corr, stats def _scan_angle(las): """scan_angle (LAS 1.4) ou scan_angle_rank (LAS ≤ 1.3), None si absent.""" for name in ("scan_angle", "scan_angle_rank"): try: return np.asarray(getattr(las, name), dtype=np.float64) except AttributeError: continue return None def _scan_lines_for_file(las_file, las, z_aligned, t): """Corrections ligne à ligne d'un LAS sol, mémoïsées par (chemin, mtime).""" 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_LINES_CACHE: return _STRIP_LINES_CACHE[cache_key] angle = _scan_angle(las) result = (np.zeros(len(z_aligned)), {}) if angle is not None and len(angle) == len(z_aligned): result = _scan_line_corrections( np.asarray(las.x, dtype=np.float64), np.asarray(las.y, dtype=np.float64), z_aligned, np.asarray(las.point_source_id), t, angle) _STRIP_LINES_CACHE[cache_key] = result 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: 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_OFFSETS_CACHE: return _STRIP_OFFSETS_CACHE[cache_key] try: psid = np.asarray(las.point_source_id) except AttributeError: psid = None # dimension absente (producteur tiers) : pas de calage if psid is not None and len(psid) == len(las.points): offsets = _strip_vertical_offsets( np.asarray(las.x, dtype=np.float64), np.asarray(las.y, dtype=np.float64), np.asarray(las.z, dtype=np.float64), psid) else: offsets = {} _STRIP_OFFSETS_CACHE[cache_key] = offsets return offsets def _write_strip_align_sidecar(dtm_dir, basename, output_suffix, offsets, jitter=None, lines=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/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, "line_window": STRIP_LINE_WINDOW, "line_cell": STRIP_LINE_CELL, "line_model": STRIP_LINE_MODEL, "lines": {str(p): {"lines": n, "rms_m": round(r, 4), "max_m": round(m, 4)} for p, (n, r, m) in (lines or {}).items()}, } try: sidecar = Path(dtm_dir) / f"{basename}_dtm{output_suffix}_stripalign.json" sidecar.write_text(json.dumps(payload, ensure_ascii=False), encoding="utf-8") except Exception as e: logger.warning(f" Écriture sidecar calage faisceaux impossible: {e}") def _strip_lidar_ext(path): """Extract base name from a LAZ/LAS file (mirrors pipeline._file_basename).""" name = Path(path).name for ext in ('.copc.laz', '.copc.las', '.laz', '.las'): if name.lower().endswith(ext): return name[:-len(ext)] return Path(path).stem def parse_ign_classes(spec): """Convertit une liste de classes IGN (noms ou codes) en codes LAS triés. Args: spec: Chaîne séparée par virgules, ex. "sol,unclassified" ou "2,1". Returns: Liste triée de codes LAS uniques. Raises: ValueError: Si un élément n'est ni un nom connu ni un code LAS 0-255, ou si la liste est vide. """ codes = set() for token in str(spec).split(","): token = token.strip().lower() if not token: continue if token in IGN_CLASS_NAMES: codes.add(IGN_CLASS_NAMES[token]) else: try: code = int(token) except ValueError: raise ValueError( f"Classe IGN inconnue: '{token}' " f"(noms: {', '.join(sorted(IGN_CLASS_NAMES))} ou code LAS 0-255)") if not 0 <= code <= 255: raise ValueError(f"Code LAS hors bornes (0-255): {code}") codes.add(code) if not codes: raise ValueError("Aucune classe IGN fournie") return sorted(codes) def ign_method_label(codes): """Étiquette de méthode encodant les classes IGN (ex. 'ign_1_2'). 'ign' seul = sol uniquement (code 2), rétrocompatible avec les fichiers de classification existants. Toute autre combinaison est encodée dans le nom pour invalider le cache et déclencher la reclassification. """ codes = sorted(codes) if codes == [2]: return "ign" return "ign_" + "_".join(str(c) for c in codes) def _create_ground_pipeline(input_laz, output_las, method, ign_codes=None): """Create a PDAL pipeline JSON for ground classification. All methods include a ReturnNumber/NumberOfReturns >= 1 filter to handle LiDAR HD files that may contain points with invalid return numbers. Pre-processing steps (PDAL recommended workflow): 1. Reset Classification to 0 2. ELM (Extended Local Minimum) — mark low outliers as noise (Classification=7) 3. Statistical outlier removal 4. Ground classification (SMRF or CSF) 5. Extract ground points (Classification=2) Args: input_laz: Path to input LAZ/LAS file. output_las: Path to output classified LAS file. method: Ground classification method ('ign', 'smrf' or 'csf'). ign_codes: LAS class codes to extract with the 'ign' method (default: [2] = sol). Multiple ranges on Classification are logically ORed by filters.range (documented PDAL semantics). Returns: JSON string of the PDAL pipeline. """ # Common ReturnNumber filter for LiDAR HD compatibility return_filter = { "type": "filters.range", "limits": "ReturnNumber[1:],NumberOfReturns[1:]" } # Classification filter (ground points only) ground_filter = { "type": "filters.range", "limits": "Classification[2:2]" } # LiDAR HD IGN : le fichier est pré-classifié par le fournisseur. # On réutilise la classification telle quelle (mode pur) : les classes # extraites sont paramétrables — par défaut le sol seul (2), mais on peut # ajouter p.ex. unclassified (1) pour combler les trous sans retouche. # Les plages multiples sur Classification sont combinées en OU logique # par filters.range (sémantique PDAL documentée). if method == 'ign': codes = sorted(ign_codes) if ign_codes else [2] ign_filter = { "type": "filters.range", "limits": ",".join(f"Classification[{c}:{c}]" for c in codes) } pipeline = { "pipeline": [ str(input_laz), return_filter, ign_filter, { "type": "writers.las", "filename": str(output_las), "extra_dims": "all" } ] } return json.dumps(pipeline) # Reset Classification to 0 before preprocessing reset_classification = { "type": "filters.assign", "assignment": "Classification[:]=0" } # ELM (Extended Local Minimum) — mark low outliers as noise # Parameters tuned for rocky limestone terrain with low vegetation: # - cell=5.0m: fine resolution to capture rocky relief # - threshold=2.0m: high threshold to avoid marking rock outcrops as noise elm_filter = { "type": "filters.elm", "cell": 5.0, "threshold": 2.0 } # Statistical outlier removal outlier_filter = { "type": "filters.outlier", "method": "statistical", "mean_k": 8, "multiplier": 3.0 } # Method-specific ground classification filter if method == 'smrf': ground_step = { "type": "filters.smrf", "ignore": "Classification[7:7]", "slope": 1.0, "window": 16.0, "threshold": 0.5, "scalar": 1.25 } elif method == 'csf': # resolution 1.0 m : un cloth à 0.5 m (= 4 M particules pour 1 km²) # rend la classification ~4× plus lente sans gain visible sur le MNT # (la résolution finale du MNT vient de la rasterisation, pas du cloth). ground_step = { "type": "filters.csf", "resolution": 1.0, "rigidness": 3, "smooth": True, "threshold": 0.5 } else: raise ValueError(f"Méthode de classification inconnue: {method}") pipeline = { "pipeline": [ str(input_laz), return_filter, reset_classification, elm_filter, outlier_filter, ground_step, ground_filter, { "type": "writers.las", "filename": str(output_las), "extra_dims": "all" } ] } return json.dumps(pipeline) def create_smrf_pipeline(input_laz, output_las): """Create a PDAL pipeline JSON for SMRF ground classification.""" return _create_ground_pipeline(input_laz, output_las, 'smrf') def create_ign_pipeline(input_laz, output_las): """Create a PDAL pipeline JSON using the IGN supplier pre-classification.""" return _create_ground_pipeline(input_laz, output_las, 'ign') def create_csf_pipeline(input_laz, output_las): """Create a PDAL pipeline JSON for CSF ground classification.""" return _create_ground_pipeline(input_laz, output_las, 'csf') def validate_laz(laz_file): """Integrity check for a LAZ/LAS file. Verifies that both the header AND point data are readable. Some corrupted COPC files have valid headers but inaccessible point data (LazrsError: failed to fill whole buffer). Such files must be re-downloaded. Returns: True if file is readable and contains accessible points, False otherwise. """ import laspy try: with laspy.open(str(laz_file)) as f: header = f.header point_count = header.point_count if point_count == 0: logger.error(f" ✗ Fichier vide (0 points): {laz_file.name}") logger.error(f" → Re-télécharger depuis https://ign.fr/lidar-hd") return False # Verify point data is actually accessible (not just header metadata) try: for _ in f.chunk_iterator(1000): break # Read just one chunk to confirm data integrity except Exception as chunk_err: logger.error(f" ✗ Données inaccessibles (fichier corrompu?): {laz_file.name}") logger.error(f" Erreur: {chunk_err}") logger.error(f" → Re-télécharger depuis https://ign.fr/lidar-hd") return False return True except Exception: pass # Fallback: try PDAL (handles COPC v1.1 that laspy can't read) try: result = subprocess.run( ["pdal", "info", str(laz_file), "--summary"], capture_output=True, text=True, timeout=30 ) if result.returncode == 0: # Check point count from PDAL info output import json as _json try: info = _json.loads(result.stdout) count = info.get('summary', {}).get('num_points', 0) if count == 0: logger.error(f" ✗ Fichier vide (0 points PDAL): {laz_file.name}") logger.error(f" → Re-télécharger depuis https://ign.fr/lidar-hd") return False except Exception: pass # Can't parse — assume valid # Verify PDAL can actually read point data (not just header) try: test_result = subprocess.run( ["pdal", "info", str(laz_file), "--point", "1"], capture_output=True, text=True, timeout=60 ) if test_result.returncode != 0: logger.error(f" ✗ Données inaccessibles (PDAL): {laz_file.name}") logger.error(f" → Re-télécharger depuis https://ign.fr/lidar-hd") return False except (subprocess.TimeoutExpired, FileNotFoundError): pass # Timeout — assume valid, will fail later if corrupted return True logger.error(f" ✗ Fichier illisible: {laz_file.name}") logger.error(f" PDAL: {result.stderr.strip()[:200]}") except (subprocess.TimeoutExpired, FileNotFoundError): logger.error(f" ✗ Impossible de vérifier le fichier: {laz_file.name}") logger.error(f" → Re-télécharger depuis https://ign.fr/lidar-hd") return False def _read_with_pdal(laz_file): """Read a LAZ/LAS file via PDAL when laspy fails (e.g. COPC v1.1). Returns a laspy.LasData object, or None on failure. """ import subprocess import tempfile import os tmp_path = None try: # Convert COPC to LAS via PDAL, then read with laspy with tempfile.NamedTemporaryFile(suffix='.las', delete=False) as tmp: tmp_path = tmp.name pipeline = json.dumps({ "pipeline": [ str(laz_file), {"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: logger.warning(f" PDAL conversion échouée: {result.stderr[:200]}") return None import laspy las = laspy.read(tmp_path) if len(las.points) == 0: logger.warning(f" PDAL: conversion réussie mais 0 points") return None return las except Exception as e: logger.warning(f" PDAL fallback échoué: {e}") return None finally: # Toujours supprimer le LAS temporaire, même en cas d'exception # (timeout subprocess, laspy.read échoué…) : sinon fuite de ~Go. if tmp_path: try: os.unlink(tmp_path) except Exception: pass def detect_ground_method(laz_file): """Detect the best ground classification method based on point cloud statistics. Auto-selects between SMRF and CSF: - SMRF: fast, robust for most natural terrain (PDAL recommended default) - CSF: cloth simulation, better for complex/urban terrain Falls back to SMRF if the file cannot be read or attributes are missing. Args: laz_file: Path to input LAZ/LAS file. Returns: String: 'ign', 'smrf' or 'csf' """ import laspy # Try laspy first, then PDAL for COPC files las = None try: las = laspy.read(str(laz_file)) except Exception as e: logger.warning(f" laspy: {e}") logger.info(f" → Lecture via PDAL pour auto-détection...") las = _read_with_pdal(laz_file) if las is None: logger.info(f" → Méthode: SMRF (défaut — lecture impossible)") return 'smrf' _LAST_READ.clear() _LAST_READ[str(laz_file)] = las # réutilisé par l'extraction IGN total_points = len(las.points) if total_points == 0: logger.warning(f" Nuage vide (0 points) — méthode par défaut: SMRF") return 'smrf' # LiDAR HD IGN : les données livrées sont pré-classifiées par le fournisseur # (classe 2 = sol). C'est la base la plus rapide (~10 s) et de référence. # Le MNT est ensuite complété par le retour le plus bas par cellule + # interpolation (voir create_dtm_fast), ce qui « rattrape » les trous de la # pré-classification (forêt dense / relief). On la préfère donc dès qu'une # part raisonnable des points est classée sol, plutôt que de refiltrer. try: cls = np.asarray(las.classification, dtype=np.int32) ground_ratio = float(np.mean(cls == 2)) except Exception: ground_ratio = 0.0 if ground_ratio >= 0.2: logger.info(f" → Méthode: IGN (pré-classification fournisseur — " f"{ground_ratio * 100:.1f}% de points classe 2)") return 'ign' z = np.array(las.z) # Height variance (always available) z_std = float(np.std(z)) z_range = float(np.max(z) - np.min(z)) # Try to get NumberOfReturns (may not exist in all point formats) single_return_ratio = 0.0 try: num_returns = np.array(las.NumberOfReturns) single_return_count = int(np.sum(num_returns == 1)) single_return_ratio = single_return_count / total_points if total_points > 0 else 0 except AttributeError: logger.debug(" NumberOfReturns non disponible — utilisation de la variance Z uniquement") logger.info(f" Analyse du nuage: {total_points:,} points, " f"ratio_retours_uniques={single_return_ratio:.2f}, " f"écart_Z={z_std:.1f}m, amplitude_Z={z_range:.1f}m") # Decision logic: # - High single-return ratio (>0.6) → urban (buildings, roads) → CSF (cloth simulation) # - High elevation variance (>30m) → complex/mountainous terrain → CSF # - Default → SMRF (fast, robust for most natural terrain) if single_return_ratio > 0.6: method = 'csf' reason = f"ratio retours uniques={single_return_ratio:.2f} > 0.6 → milieu urbain" elif z_std > 30: method = 'csf' reason = f"écart_Z={z_std:.1f}m > 30m → terrain complexe" else: method = 'smrf' reason = f"terrain naturel standard" logger.info(f" → Méthode: {method.upper()} ({reason})") return method _LAST_READ = {} def _extract_ign_ground(laz_file, output_las, codes): """Extraction directe (laspy) des classes IGN choisies vers un LAS. Même résultat que le pipeline PDAL (filtres ReturnNumber ≥ 1, NumberOfReturns ≥ 1, classes) mais sans relecture ni conversion : ~4 s au lieu de ~9 s par dalle, et la lecture de la détection automatique est réutilisée. Format de points et échelles du fichier source conservés. Returns: True si le fichier a été écrit avec au moins un point. """ import laspy las = _LAST_READ.pop(str(laz_file), None) if las is None: las = laspy.read(str(laz_file)) keep = ((np.asarray(las.return_number) >= 1) & (np.asarray(las.number_of_returns) >= 1) & np.isin(np.asarray(las.classification), np.asarray(sorted(codes)))) if not keep.any(): return False header = laspy.LasHeader(point_format=las.header.point_format, version=las.header.version) header.scales = las.header.scales header.offsets = las.header.offsets out = laspy.LasData(header) out.points = las.points[keep] out.write(str(output_las)) return True def classify_ground(laz_file, temp_dir, method='auto', force=False, ign_classes="sol"): """Classify ground points using PDAL ground classification filter. Args: laz_file: Path to input LAZ/LAS file. temp_dir: Directory for temporary files (pipeline.json, ground.las). method: Ground classification method ('auto', 'ign', 'smrf' or 'csf'). force: If True, reclassify even if output file already exists. ign_classes: Classes LAS extraites par la méthode IGN (noms ou codes séparés par virgules, ex. "sol,unclassified"). Ignoré pour les autres méthodes. Returns: Path to classified ground LAS file, or None on failure. """ import laspy # Auto-detect method if requested if method == 'auto': method = detect_ground_method(laz_file) logger.info(f" Classification sol: {method.upper()} (auto)") else: logger.info(f" Classification sol: {method.upper()} (forcé)") # Les classes IGN sont encodées dans le nom de fichier (ex. ign_1_2) # pour qu'un changement de classes invalide le cache et déclenche la # reclassification. ign_codes = parse_ign_classes(ign_classes) if method == 'ign' else None method_label = ign_method_label(ign_codes) if ign_codes else method laz_base = _strip_lidar_ext(laz_file) output_las = temp_dir / f"{laz_base}_ground_{method_label}.las" if output_las.exists() and not force: logger.info(f" Classification {method.upper()} déjà effectuée — fichier existant réutilisé") return output_las if force and output_las.exists(): logger.info(f" Reclassification forcée — suppression de {output_las.name}") output_las.unlink() # Pré-classification IGN : extraction directe, PDAL en secours. if method == 'ign': try: if _extract_ign_ground(laz_file, output_las, ign_codes or [2]): logger.info(f" ✓ Classification sol IGN terminée (extraction directe)") return output_las logger.warning(" Aucun point des classes IGN demandées — repli SMRF") output_las.unlink(missing_ok=True) return _fallback_to_smrf(laz_file, temp_dir, laz_base, force, source=method_label) except Exception as e: logger.warning(f" Extraction IGN directe impossible ({e}) — pipeline PDAL") output_las.unlink(missing_ok=True) _LAST_READ.clear() pipeline_json = _create_ground_pipeline(laz_file, output_las, method, ign_codes=ign_codes) pipeline_file = temp_dir / f"pipeline_{method_label}.json" with open(pipeline_file, 'w') as f: f.write(pipeline_json) try: subprocess.run( ["pdal", "pipeline", str(pipeline_file)], capture_output=True, check=True ) # Verify that ground file has points if output_las.exists() and output_las.stat().st_size < 100: logger.error(f" ✗ Fichier ground vide (taille < 100 octets)") output_las.unlink(missing_ok=True) # Fallback: si la méthode ne produit aucun point sol, réessayer avec SMRF if method in ('csf', 'ign'): return _fallback_to_smrf(laz_file, temp_dir, laz_base, force, source=method_label) return None logger.info(f" ✓ Classification sol {method.upper()} terminée") return output_las except subprocess.CalledProcessError as e: error_msg = e.stderr.decode() if e.stderr else str(e) logger.warning(f" ✗ Erreur classification PDAL ({method.upper()}): {error_msg}") # Fallback: si CSF ou la pré-classification échouent, réessayer avec SMRF if method in ('csf', 'ign'): return _fallback_to_smrf(laz_file, temp_dir, laz_base, force, source=method_label) # Try repairing file with laspy if PDAL fails on EVLR/VLR if 'VLR' in error_msg or 'Invalid' in error_msg: logger.info(f" → Tentative de réparation du fichier avec laspy...") repaired_las = temp_dir / f"{laz_base}_repaired.las" if _repair_laz_with_laspy(laz_file, repaired_las): # Retry PDAL pipeline with repaired file pipeline_json = _create_ground_pipeline(repaired_las, output_las, method) with open(pipeline_file, 'w') as f: f.write(pipeline_json) try: subprocess.run( ["pdal", "pipeline", str(pipeline_file)], capture_output=True, check=True ) logger.info(f" ✓ Classification sol {method.upper()} terminée (fichier réparé)") return output_las except subprocess.CalledProcessError as e2: error_msg2 = e2.stderr.decode() if e2.stderr else str(e2) logger.error(f" ✗ Échec classification même après réparation: {error_msg2}") else: logger.error(f" ✗ Impossible de réparer le fichier") return None def _fallback_to_smrf(laz_file, temp_dir, laz_base, force=False, source='csf'): """Retry ground classification with SMRF when CSF/IGN fails. CSF (Cloth Simulation Filter) can fail on certain terrain types where SMRF (Simple Morphological Filter) succeeds, and a file without usable pre-classification produces an empty ground extract. This fallback ensures processing continues even when the selected method fails. Args: laz_file: Path to input LAZ/LAS file. temp_dir: Directory for temporary files. laz_base: Base name for the file. force: If True, reclassify even if output exists. source: Method that failed ('csf' or 'ign'). Returns: Path to classified ground LAS file, or None on failure. """ logger.info(f" → Basculement {source.upper()} → SMRF (fallback)") # Clean up failed output if it exists failed_output = temp_dir / f"{laz_base}_ground_{source}.las" if failed_output.exists(): failed_output.unlink(missing_ok=True) output_las = temp_dir / f"{laz_base}_ground_smrf.las" if output_las.exists() and not force: logger.info(f" Classification SMRF déjà existante — fichier réutilisé") return output_las pipeline_json = _create_ground_pipeline(laz_file, output_las, 'smrf') pipeline_file = temp_dir / "pipeline_smrf.json" with open(pipeline_file, 'w') as f: f.write(pipeline_json) try: subprocess.run( ["pdal", "pipeline", str(pipeline_file)], capture_output=True, check=True ) if output_las.exists() and output_las.stat().st_size < 100: logger.error(f" ✗ Fichier ground SMRF vide (taille < 100 octets)") output_las.unlink(missing_ok=True) return None logger.info(f" ✓ Classification sol SMRF terminée (fallback depuis CSF)") return output_las except subprocess.CalledProcessError as e: error_msg = e.stderr.decode() if e.stderr else str(e) logger.error(f" ✗ Échec classification SMRF (fallback): {error_msg}") return None def _repair_laz_with_laspy(input_laz, output_las): """Try to repair a corrupt LAZ file by re-reading with laspy and saving as LAS. Works around PDAL errors like 'Invalid Extended VLR size' by stripping problematic VLR/EVLR metadata during re-save. Args: input_laz: Path to corrupt LAZ/LAS file. output_las: Path for repaired LAS output. Returns: True if repair succeeded, False otherwise. """ import laspy try: las = laspy.read(str(input_laz)) las.write(str(output_las)) logger.info(f" ✓ Fichier réparé via laspy ({len(las.points):,} points)") return True except Exception as e: logger.warning(f" ✗ Réparation laspy échouée: {e}") return False # ------------------------------------------------------------ # Comblement des vides entre points (petits trous) # ------------------------------------------------------------ # À 0,2 m, ~80 % des pixels ne reçoivent aucun point : le MNT est comblé entre # les points mesurés. Un comblement « à distance fixe » (tout pixel à moins de # 1 m d'un point) faisait grossir chaque point isolé en pastille plate de 2 m # et extrapolait une bande de 1 m dans chaque trou (plateaux recopiés, fausses # pentes : pastilles roses et liserés colorés dans le relief orienté). Le # comblement est désormais borné à l'ENVELOPPE des points : fermeture # morphologique dont le rayon suit l'espacement local des points (court en # zone dense, long sous couvert clairsemé), puis suppression des îlots isolés. GAP_FILL_TAG = "LIDAR_GAP_FILL" # tag GeoTIFF : version du comblement GAP_FILL_VERSION = 3 # 1 (absent) = distance fixe 1 m ; 3 = + densité GAP_RADIUS_K = 1.5 # rayon = K × espacement local des points GAP_RADII_M = (1.0, 1.5, 2.0, 3.0) # paliers de rayon (m) ; 1 m mini : trous < 2 m comblés GAP_DENSITY_WINDOW_M = 5.0 # fenêtre de mesure de la densité locale GAP_MIN_ISLAND_M2 = 1.0 # îlots de données plus petits : supprimés def _morph_step(mask, square, erode): """Un pas 3×3 (carré ou croix) de dilatation/érosion binaire par tranches numpy (~10× plus rapide que scipy.ndimage). Au bord, le voisin manquant compte comme le pixel lui-même.""" op = np.logical_and if erode else np.logical_or out = mask.copy() op(out[1:], mask[:-1], out=out[1:]) op(out[:-1], mask[1:], out=out[:-1]) src = out.copy() if square else mask # carré : séparable (colonnes puis lignes) op(out[:, 1:], src[:, :-1], out=out[:, 1:]) op(out[:, :-1], src[:, 1:], out=out[:, :-1]) return out def _morph_disk(mask, radius_px, erode=False, first_step=1): """Dilatation/érosion par un disque approché : octogone, pas carrés (impairs) et croix (pairs) alternés. `first_step` poursuit une chaîne.""" for step in range(first_step, first_step + radius_px): mask = _morph_step(mask, step % 2 == 1, erode) return mask def _gap_fill_envelope(valid, resolution): """Masque des pixels à combler ou conserver : enveloppe des points mesurés. Fermeture morphologique (dilatation puis érosion) de rayon variable : les vides ENTRE points voisins sont comblés, rien n'est étendu vers l'extérieur (un point isolé reste un pixel). Le rayon vaut GAP_RADIUS_K × l'espacement local des points, mesuré sur GAP_DENSITY_WINDOW_M dans la zone couverte (le bord d'un trou ne fait pas chuter la densité), arrondi au palier GAP_RADII_M le plus proche. Les îlots de moins de GAP_MIN_ISLAND_M2 sont retirés (retours épars sur l'eau, dans une cour…). """ from scipy import ndimage as nd radii_px = sorted({max(1, int(round(r / resolution))) for r in GAP_RADII_M}) r_max = radii_px[-1] # Masque étendu en miroir : la bordure de la grille n'érode pas les # données et ne comble pas les bandes vides. Une seule chaîne de # dilatations sert tous les paliers (disque de rayon r = r premiers pas). pad = 2 * r_max dilated = {} grown = np.pad(valid, pad, mode="symmetric") for step in range(1, r_max + 1): grown = _morph_disk(grown, 1, first_step=step) if step in radii_px: dilated[step] = grown del grown # Densité locale : part de pixels mesurés dans la zone couverte (points à # moins du plus grand rayon), convertie en espacement moyen des points. # Calcul sur une grille de ~1 m (fenêtre de 5 m : largement suffisant). covered = dilated[r_max][pad:-pad, pad:-pad] h, w = valid.shape f = max(1, int(round(1.0 / resolution))) hc, wc = -(-h // f), -(-w // f) def _block_mean(a): full = np.zeros((hc * f, wc * f), dtype=np.float32) full[:h, :w] = a return full.reshape(hc, f, wc, f).mean(axis=(1, 3)) win = max(3, int(round(GAP_DENSITY_WINDOW_M / (resolution * f))) | 1) frac = nd.uniform_filter(_block_mean(valid), win, mode="nearest") frac_cov = nd.uniform_filter(_block_mean(covered), win, mode="nearest") with np.errstate(divide="ignore", invalid="ignore"): # rayon voulu en pixels = K × espacement / résolution wanted = GAP_RADIUS_K / np.sqrt(frac / np.maximum(frac_cov, 1e-6)) # Rayon voulu agrandi en bilinéaire (sinon changements de palier en # créneaux de 1 m sur le bord des trous), puis palier le plus proche. wanted = np.nan_to_num(np.clip(wanted, 0, 2 * r_max), nan=2 * r_max) if f > 1: wanted = nd.zoom(wanted.astype(np.float32), f, order=1, grid_mode=True, mode="nearest")[:h, :w] mids = [(a + b) / 2 for a, b in zip(radii_px, radii_px[1:])] level = np.digitize(wanted, mids).astype(np.uint8) del covered, frac, frac_cov, wanted # Érosion par palier, restreinte à l'emprise de ses pixels (les grands # rayons ne concernent que les zones clairsemées) envelope = valid.copy() for i, r in enumerate(radii_px): sel = (level == i) & ~valid rows = np.flatnonzero(sel.any(axis=1)) if rows.size == 0: continue cols = np.flatnonzero(sel.any(axis=0)) y0, y1 = rows[0], rows[-1] + 1 x0, x1 = cols[0], cols[-1] + 1 m = 2 * r # marge : l'érosion au bord du découpage reste hors emprise crop = dilated[r][y0 + pad - m:y1 + pad + m, x0 + pad - m:x1 + pad + m] closed = _morph_disk(crop, r, erode=True)[m:-m, m:-m] envelope[y0:y1, x0:x1] |= sel[y0:y1, x0:x1] & closed del dilated, level # Îlots trop petits : retirés (8-connexité, surface en m²) labels, n = nd.label(envelope, structure=np.ones((3, 3), dtype=bool)) if n: area = np.bincount(labels.ravel()) * resolution * resolution small = area < GAP_MIN_ISLAND_M2 small[0] = False if small.any(): envelope &= ~small[labels] return envelope def _fill_small_gaps(dtm, resolution): """Comble les vides entre points dans l'enveloppe des données. Returns: Tuple (dtm, filled_count, removed_count) : pixels comblés et pixels mesurés retirés avec les îlots isolés. """ valid = ~np.isnan(dtm) if valid.all() or not valid.any(): return dtm, 0, 0 envelope = _gap_fill_envelope(valid, resolution) from rasterio.fill import fillnodata # Tout pixel de l'enveloppe est à moins du plus grand rayon d'un point max_px = max(1, int(round(max(GAP_RADII_M) / resolution))) # fillnodata écrit dans le tableau fourni : copie filled = fillnodata(dtm.copy(), mask=valid, max_search_distance=max_px) to_fill = envelope & ~valid & ~np.isnan(filled) removed = valid & ~envelope out = dtm.copy() out[to_fill] = filled[to_fill] out[removed] = np.nan return out, int(to_fill.sum()), int(removed.sum()) # Densité de points sol (couche « densite_sol ») : comptage des points retenus # pour le MNT en mailles de DENSITY_CELL_M, moyenné sur DENSITY_SMOOTH × # DENSITY_SMOOTH mailles (pts/m² sur 9 m² : moins de bruit de petits entiers). # Écrite à côté du DTM (*_dtm*_density.tif) : les visualisations ne reçoivent # que le MNT, plus les points. DENSITY_CELL_M = 1.0 DENSITY_SMOOTH = 3 def density_path(dtm_path): """Fichier annexe de densité d'un DTM.""" dtm_path = Path(dtm_path) return dtm_path.with_name(f"{dtm_path.stem}_density.tif") def _write_density(xs, ys, bounds, dtm_path): """Écrit la densité de points sol (pts/m², GeoTIFF float32 à 1 m).""" from scipy import ndimage as nd min_x, min_y, max_x, max_y = bounds cell = DENSITY_CELL_M nx = max(1, int(np.ceil((max_x - min_x) / cell - 1e-6))) ny = max(1, int(np.ceil((max_y - min_y) / cell - 1e-6))) top = max_y counts, _, _ = np.histogram2d(top - ys, xs - min_x, bins=(ny, nx), range=[[0, ny * cell], [0, nx * cell]]) density = nd.uniform_filter(counts, DENSITY_SMOOTH, mode="nearest") / (cell * cell) out = density_path(dtm_path) with rasterio.open( out, 'w', driver='GTiff', height=ny, width=nx, count=1, dtype='float32', crs='EPSG:2154', transform=from_bounds(min_x, top - ny * cell, min_x + nx * cell, top, nx, ny), compress='deflate', ) as dst: dst.write(density.astype('float32'), 1) return out def read_dtm_gap_fill(dtm_path): """Version du comblement enregistrée dans un DTM (1 si absente/illisible).""" try: with rasterio.open(dtm_path) as src: return int(src.tags().get(GAP_FILL_TAG, 1) or 1) except Exception: return 1 def _interpolate_holes(dtm, downsample=8): """Fill remaining NaN holes with a terrain-aware surface interpolation. In complex / rocky terrain the ground under-classification leaves interior holes far too large for a 1 m gap fill, which otherwise become flat nearest-neighbor patches in the downstream layers. This helper triangulates the valid cells on a downsampled grid (linear, nearest as a fallback for cells outside the data hull) and bilinearly upsamples the result, keeping the operation fast even for large rasters. Args: dtm: 2-D float array (may contain NaN holes). downsample: Coarsening factor for the interpolation grid. Returns: Tuple (filled_array, filled_count). """ holes = np.isnan(dtm) if not holes.any(): return dtm, 0 valid = ~holes if not valid.any(): return dtm, 0 height, width = dtm.shape step = max(1, downsample) coarse = dtm[::step, ::step].astype(np.float64) c_valid = ~np.isnan(coarse) c_holes = np.isnan(coarse) if not c_holes.any() or not c_valid.any(): return dtm, 0 from scipy.interpolate import griddata from scipy.ndimage import map_coordinates cy, cx = np.where(c_valid) c_coords = np.column_stack([cx, cy]).astype(np.float64) c_vals = coarse[c_valid] hy, hx = np.where(c_holes) h_coords = np.column_stack([hx, hy]).astype(np.float64) interp = griddata(c_coords, c_vals, h_coords, method='linear') bad = np.isnan(interp) if bad.any(): interp[bad] = griddata(c_coords, c_vals, h_coords[bad], method='nearest') coarse_filled = coarse.copy() coarse_filled[c_holes] = interp # Coarse cell i represents fine column/row i*step, so fine index c maps to # coarse coordinate c/step (no half-cell offset). rows = np.arange(height) / step cols = np.arange(width) / step grid_y, grid_x = np.meshgrid(rows, cols, indexing='ij') upsampled = map_coordinates(coarse_filled, [grid_y, grid_x], order=1) filled = dtm.copy() filled[holes] = upsampled[holes] return filled, int(holes.sum()) # ------------------------------------------------------------ # 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) # 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)] 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`. 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 # 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 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: 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: 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: (xs, ys, zs) concaténés (tableaux vides si aucune voisine). """ 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 codes = sorted(set(int(c) for c in classes)) or [2] limits = ",".join(f"Classification[{c}:{c}]" for c in codes) 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 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 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, pure=False, strip_align=True, edge_buffer=0.0, neighbor_classes=None): """Create DTM using fast binning method with gap filling. Args: las_file: Path to classified ground LAS file. basename: Base name for output file. dtm_dir: Directory for output DTM GeoTIFF. 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 : 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 : 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. """ output_tif = dtm_dir / f"{basename}_dtm{output_suffix}.tif" if output_tif.exists() and not force: logger.info(f" DTM déjà existant — fichier réutilisé: {output_tif.name}") return output_tif import laspy logger.info(" → Génération DTM...") try: t_read = time.perf_counter() las = laspy.read(str(las_file)) logger.info(f" Lecture {len(las.points):,} points " f"({time.perf_counter() - t_read:.1f}s)") except Exception as e: # laspy can't read COPC v1.1 — try PDAL conversion logger.warning(f" laspy: {e}") logger.info(f" → Conversion via PDAL pour lecture COPC...") las = _read_with_pdal(las_file) if las is None: logger.error(f" ✗ Impossible de lire {las_file.name}") return None if len(las.points) == 0: logger.error(f" ✗ Fichier vide (0 points): {las_file.name}") return None # 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 t_align = time.perf_counter() if strip_align: try: strip_offsets = _strip_offsets_for_file(las_file, las) except Exception as e: logger.warning(f" Mesure du calage faisceaux impossible ({e}) — MNT non calé") strip_offsets = {} if strip_offsets: logger.info(" Calage faisceaux : " + ", ".join( f"PSID {p} {off:+.3f} m" for p, off in sorted(strip_offsets.items()))) 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 # L'ajustement conjoint ligne à ligne (3ᵉ passe) recale déjà chaque # ligne contre les autres faisceaux, à toutes les échelles : la gigue # par fenêtres de temps n'est calculée que si scan_angle manque. lines_possible = _scan_angle(las) is not None if gps_time is not None and len(gps_time) == len(las.points) and not lines_possible: 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") logger.info(f" Calage faisceaux : {time.perf_counter() - t_align:.1f}s") try: min_x, max_x = float(las.header.min[0]), float(las.header.max[0]) min_y, max_y = float(las.header.min[1]), float(las.header.max[1]) 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)...") xs = np.asarray(las.x, dtype=np.float64) ys = np.asarray(las.y, dtype=np.float64) zs = np.asarray(las.z, dtype=np.float64) if strip_offsets: lut = np.zeros(65536) 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) # 3ᵉ passe : décalage ligne à ligne (après les deux premières). strip_lines = {} if strip_align and gps_time is not None and len(gps_time) == len(zs): t_lines = time.perf_counter() try: line_corr, strip_lines = _scan_lines_for_file(las_file, las, zs, gps_time) zs = zs - line_corr except Exception as e: logger.warning(f" Mesure du décalage ligne à ligne impossible ({e}) — non corrigé") strip_lines = {} for p_, (nl_, rms_, max_) in sorted(strip_lines.items()): logger.info(f" Lignes PSID {p_} : rms {rms_ * 100:.1f} cm, max {max_ * 100:.1f} cm " f"({nl_} lignes)") logger.info(f" Décalage ligne à ligne : {time.perf_counter() - t_lines:.1f}s") # 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: t_neigh = time.perf_counter() 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]) logger.info(f" Voisines (raccord) : {len(nx):,} points " f"({time.perf_counter() - t_neigh:.1f}s)") if len(nx): xs = np.concatenate([xs, nx]) ys = np.concatenate([ys, ny]) zs = np.concatenate([zs, nz]) # Densité des points sol retenus (couche densite_sol), best-effort try: _write_density(xs, ys, (min_x, min_y, max_x, max_y), dtm_dir / f"{basename}_dtm{output_suffix}.tif") except Exception as e: logger.warning(f" Densité de points non écrite ({e})") t_raster = time.perf_counter() dtm = bin_mean_2d(xs, ys, zs, width, height, (min_x, max_x), (min_y, max_y)) if dtm is not None: # Sortie déjà en (height, width) : seul le flip Y reste à faire dtm = dtm[::-1, :] # nord en haut logger.info(f" ✓ Rasterisation {width}x{height} GPU " f"({time.perf_counter() - t_raster:.1f}s)") else: stat = binned_statistic_2d( xs, ys, zs, statistic='mean', bins=[width, height], range=[[min_x, max_x], [min_y, max_y]] ) dtm = stat.statistic.T dtm = dtm[::-1, :] # Flip Y so north is at top logger.info(f" ✓ Rasterisation {width}x{height} CPU " f"({time.perf_counter() - t_raster:.1f}s)") # 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. # 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. # Vides entre points comblés dans l'enveloppe des données (rayon # adapté à la densité locale), îlots isolés retirés nan_count = np.count_nonzero(np.isnan(dtm)) if nan_count > 0: total = dtm.size nan_pct = 100.0 * nan_count / total logger.info(f" {nan_count:,} pixels sans données ({nan_pct:.1f}%)") t_fill = time.perf_counter() dtm, filled_count, removed_count = _fill_small_gaps(dtm, resolution) logger.info(f" {filled_count:,} pixels comblés entre points, " f"{removed_count:,} pixels d'îlots isolés retirés " f"({time.perf_counter() - t_fill:.1f}s)") # Save as GeoTIFF output_tif = dtm_dir / f"{basename}_dtm{output_suffix}.tif" transform = from_bounds(min_x, min_y, max_x, max_y, width, height) t_write = time.perf_counter() with rasterio.open( output_tif, 'w', driver='GTiff', height=height, width=width, count=1, dtype='float32', crs='EPSG:2154', transform=transform, nodata=float('nan'), compress='lzw' ) as dst: dst.write(dtm.astype('float32'), 1) # Version du comblement : un DTM d'une version antérieure est # régénéré (pastilles et liserés de l'ancien comblement). dst.update_tags(**{GAP_FILL_TAG: str(GAP_FILL_VERSION)}) 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}"}) logger.info(f" Écriture GeoTIFF : {time.perf_counter() - t_write:.1f}s") if strip_align: _write_strip_align_sidecar(dtm_dir, basename, output_suffix, strip_offsets, strip_jitter, strip_lines) logger.info(f" ✓ DTM créé: {output_tif.name}") return output_tif except Exception as e: logger.error(f" ✗ Erreur DTM: {e}", exc_info=True) return None