d32502b74e
Troisième passe du calage vertical : chaque ligne de balayage (décalage et inclinaison due au roulis) est recalée contre le consensus des autres faisceaux, à toutes les échelles, avec un profil d'étalonnage par faisceau et par degré d'angle qui retire les écarts non linéaires en travers de la fauchée. Les lignes sans recouvrement sont corrigées contre leur propre faisceau. Efface les lignes en creux et la marche au bord de fauchée mesurées sur LHD_FXX_0999_6882 (validé sur des blocs jamais vus). Calcul vectorisé, CuPy si GPU ; la gigue par fenêtres de temps devient inutile quand scan_angle existe. Rendu plus rapide : encodage AVIF speed 9 (0,6 s au lieu de 4 s par dalle), classification IGN par extraction directe laspy au lieu de PDAL (4,9 s au lieu de 13,5 s), comblement des trous et gradients sur GPU, cache numba persistant dans l'image. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
1911 lines
80 KiB
Python
1911 lines
80 KiB
Python
"""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
|
||
|
||
|
||
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])
|
||
|
||
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.
|
||
|
||
# Fill small gaps (< 1 m from data) precisely — comme avant
|
||
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}%)")
|
||
|
||
max_gap_pixels = max(1, int(1.0 / resolution))
|
||
t_fill = time.perf_counter()
|
||
from rasterio.fill import fillnodata
|
||
valid_mask = ~np.isnan(dtm)
|
||
dtm_filled = fillnodata(dtm, mask=valid_mask, max_search_distance=max_gap_pixels)
|
||
small_gap_mask = np.isnan(dtm) & ~np.isnan(dtm_filled)
|
||
filled_count = np.count_nonzero(small_gap_mask)
|
||
if filled_count > 0:
|
||
dtm = np.where(small_gap_mask, dtm_filled, dtm)
|
||
logger.info(f" {filled_count:,} petits trous comblés "
|
||
f"(< {max_gap_pixels}px, {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)
|
||
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 |