5f018fe2b6
Le MNT de chaque dalle s'étend d'une bande de 100 m remplie avec les points sol des 8 tuiles voisines (option « Raccord des bords ») : les rendus à grand noyau (openness, SVF, LRM) deviennent continus d'une dalle à l'autre, les images restent recadrées sur le kilomètre exact. Pendant une génération, la carte encadre les dalles du run : orange pulsant en cours de rendu, rouge en échec — les coins WGS84 sont portés par /api/status pour toutes les dalles non terminées. La file de génération trie les dalles du nord au sud et passe à 10 workers GPU.
1381 lines
55 KiB
Python
1381 lines
55 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. 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
|
||
from pathlib import Path
|
||
|
||
import numpy as np
|
||
import rasterio
|
||
from rasterio.transform import from_bounds
|
||
from scipy.stats import binned_statistic_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).
|
||
STRIP_ALIGN_VERSION = 2
|
||
STRIP_ALIGN_THRESHOLD = 0.005 # m : écart mini pour corriger un faisceau (0,5 cm)
|
||
STRIP_ALIGN_CELL = 1.0 # m : maille de comparaison des faisceaux
|
||
STRIP_ALIGN_MIN_SHARED = 500 # cellules sol communes mini pour valider un offset
|
||
STRIP_JITTER_BIN = 0.1 # s : durée d'une fenêtre de temps GPS (gigue)
|
||
STRIP_JITTER_SMOOTH = 5 # fenêtres : largeur de la médiane glissante
|
||
STRIP_JITTER_MIN_CELLS = 40 # cellules sol partagées mini pour valider une fenêtre
|
||
STRIP_JITTER_MAX = 0.10 # m : amplitude maxi d'une correction de gigue (garde-fou)
|
||
|
||
# Mémo des offsets par fichier : la classification est partagée entre
|
||
# résolutions, le même LAS sol est rasterisé à 0,5 m puis 0,2 m.
|
||
_STRIP_OFFSETS_CACHE = {}
|
||
_STRIP_JITTER_CACHE = {}
|
||
|
||
|
||
def _strip_surface_grid(x, y, z, inv, n_sources, cell=STRIP_ALIGN_CELL):
|
||
"""Surfaces sol par faisceau sur maille 1 m (moyenne des points bas).
|
||
|
||
Pour chaque faisceau et chaque maille : moyenne des points situés à moins
|
||
de 0,5 m du minimum du faisceau dans la maille (robuste à la végétation
|
||
résiduelle), mailles à ≥ 2 points seulement.
|
||
|
||
Returns:
|
||
(allcells, grid, key) : cellules triées, grid[faisceau, cellule] = Z
|
||
moyen (NaN si absent) et clé de maille de CHAQUE point ; None si
|
||
aucune surface n'est peuplée.
|
||
"""
|
||
x0 = np.floor(np.min(x) / cell) * cell
|
||
y0 = np.floor(np.min(y) / cell) * cell
|
||
xi = ((x - x0) / cell).astype(np.int64)
|
||
yi = ((y - y0) / cell).astype(np.int64)
|
||
ny = int(yi.max()) + 1
|
||
key = xi * ny + yi
|
||
|
||
def _surface(k):
|
||
m = inv == k
|
||
kk, zz = key[m], z[m]
|
||
order = np.argsort(kk, kind='stable')
|
||
k_s, z_s = kk[order], zz[order]
|
||
starts = np.flatnonzero(np.r_[True, k_s[1:] != k_s[:-1]])
|
||
mins = np.minimum.reduceat(z_s, starts)
|
||
ukey = k_s[starts]
|
||
thr = mins[np.searchsorted(ukey, kk)]
|
||
sel = np.flatnonzero(zz <= thr + 0.5)
|
||
k2, z2 = kk[sel], zz[sel]
|
||
cnt = np.bincount(k2)
|
||
sums = np.bincount(k2, weights=z2)
|
||
v = np.flatnonzero(cnt >= 2)
|
||
return v, sums[v] / cnt[v]
|
||
|
||
surfaces = [_surface(k) for k in range(n_sources)]
|
||
populated = [c for c, _ in surfaces if len(c)]
|
||
if not populated:
|
||
return None
|
||
allcells = np.unique(np.concatenate(populated))
|
||
grid = np.full((n_sources, len(allcells)), np.nan)
|
||
for k, (c, zs) in enumerate(surfaces):
|
||
grid[k, np.searchsorted(allcells, c)] = zs
|
||
return allcells, grid, key
|
||
|
||
|
||
def _strip_vertical_offsets(x, y, z, psid, cell=STRIP_ALIGN_CELL,
|
||
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 _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):
|
||
"""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,
|
||
}
|
||
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'
|
||
|
||
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
|
||
|
||
|
||
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()
|
||
|
||
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)
|
||
|
||
# Décalages des 8 voisines d'une dalle (col, row) en km
|
||
_NEIGHBOR_OFFSETS = [(-1, -1), (0, -1), (1, -1), (-1, 0),
|
||
(1, 0), (-1, 1), (0, 1), (1, 1)]
|
||
|
||
|
||
def _tile_coords(name):
|
||
"""Coordonnées (col, row) en km d'un nom de fichier LHD, ou None."""
|
||
from .index import parse_basename_coords
|
||
return parse_basename_coords(Path(name).name)
|
||
|
||
|
||
def _neighbor_laz_files(source_laz):
|
||
"""Liste les 8 fichiers LAZ/LAS adjacents à `source_laz` dans son dossier."""
|
||
coords = _tile_coords(source_laz)
|
||
if coords is None:
|
||
return []
|
||
col, row = coords
|
||
directory = Path(source_laz).parent
|
||
neighbors = []
|
||
for dcol, drow in _NEIGHBOR_OFFSETS:
|
||
nc, nr = col + dcol, row + drow
|
||
# Noms LHD : col/row en km sur 4 chiffres complétés (ex. 0638_6628) ;
|
||
# on accepte aussi la variante sans remplissage pour les dalles exotiques.
|
||
names = {f"LHD_FXX_{nc:04d}_{nr:04d}", f"LHD_FXX_{nc}_{nr}"}
|
||
matches = []
|
||
for name in names:
|
||
matches += list(directory.glob(f"{name}_*.las"))
|
||
matches += list(directory.glob(f"{name}_*.laz"))
|
||
if matches:
|
||
neighbors.append(sorted(matches)[0])
|
||
else:
|
||
logger.debug(f" Voisine absente : LHD_FXX_{nc:04d}_{nr:04d} (bande de bord vide)")
|
||
return neighbors
|
||
|
||
|
||
def _neighbor_ground_points(source_laz, bounds, classes):
|
||
"""Points sol des tuiles voisines dans `bounds` (raccord des bords).
|
||
|
||
Lecture PDAL en flux par voisine : découpe sur l'emprise étendue puis
|
||
filtre de classes (mêmes codes que le MNT). Best-effort : une voisine
|
||
illisible ou absente est ignorée — la bande correspondante reste vide.
|
||
|
||
Args:
|
||
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:
|
||
las = laspy.read(str(las_file))
|
||
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
|
||
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
|
||
if gps_time is not None and len(gps_time) == len(las.points):
|
||
try:
|
||
strip_jitter = _strip_jitter_for_file(las_file, las, strip_offsets)
|
||
except Exception as e:
|
||
logger.warning(f" Mesure de la gigue intra-faisceau impossible ({e}) — gigue non corrigée")
|
||
strip_jitter = {}
|
||
if strip_jitter:
|
||
for p_, (jt, jc) in sorted(strip_jitter.items()):
|
||
logger.info(f" Gigue PSID {p_} : ±{np.max(np.abs(jc)) * 100:.1f} cm "
|
||
f"(rms {np.sqrt(np.mean(np.asarray(jc) ** 2)) * 100:.1f} cm, "
|
||
f"{len(jt)} fenêtres de {STRIP_JITTER_BIN:g} s)")
|
||
else:
|
||
logger.debug(" Gigue intra-faisceau : rien à corriger")
|
||
|
||
try:
|
||
|
||
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)
|
||
|
||
# Points sol des tuiles voisines dans la bande de raccord (best-effort,
|
||
# non calés par faisceau : bande de contexte, l'image finale est
|
||
# recadrée sur la dalle avant livraison).
|
||
if used_edge_buffer > 0 and source_laz is not None:
|
||
nx, ny, nz = _neighbor_ground_points(
|
||
source_laz, (min_x, min_y, max_x, max_y),
|
||
neighbor_classes if neighbor_classes is not None else [2])
|
||
if len(nx):
|
||
xs = np.concatenate([xs, nx])
|
||
ys = np.concatenate([ys, ny])
|
||
zs = np.concatenate([zs, nz])
|
||
|
||
stat = binned_statistic_2d(
|
||
xs, ys, zs,
|
||
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
|
||
|
||
# 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))
|
||
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 (< {max_gap_pixels}px)")
|
||
|
||
# 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)
|
||
|
||
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}"})
|
||
|
||
if strip_align:
|
||
_write_strip_align_sidecar(dtm_dir, basename, output_suffix,
|
||
strip_offsets, strip_jitter)
|
||
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 |