2109 lines
89 KiB
Python
2109 lines
89 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
|
||
|
||
|
||
# ------------------------------------------------------------
|
||
# 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 |