Raccorder les bords de tuiles voisines et encadrer les rendus en cours
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.
This commit is contained in:
Antoine Jacquin
@ -42,14 +42,74 @@ IGN_CLASS_NAMES = {
|
||||
# selon la campagne, rien n'est codé en dur) et on le retranche avant
|
||||
# rastérisation. Les offsets dérivent le long d'une ligne de vol (signe inversé
|
||||
# entre tuiles voisines mesuré) : le calcul est donc par tuile, jamais global.
|
||||
STRIP_ALIGN_VERSION = 1
|
||||
#
|
||||
# 2ᵉ passe — gigue intra-faisceau : le décalage vertical peut aussi varier au
|
||||
# fil d'une MÊME passe (vibration capteur / bruit haute fréquence de la
|
||||
# trajectoire) : les lignes de balayage successives d'un faisceau apparaissent
|
||||
# alors décalées verticalement de façon aléatoire. Chaque faisceau est découpé
|
||||
# en fenêtres de temps GPS, l'offset robuste de chaque fenêtre est mesuré contre
|
||||
# la surface médiane des autres faisceaux (même maille 1 m), la série est lissée
|
||||
# (médiane glissante) puis interpolée au temps GPS de chaque point. Requiert la
|
||||
# dimension gps_time (ignorée silencieusement sinon).
|
||||
STRIP_ALIGN_VERSION = 2
|
||||
STRIP_ALIGN_THRESHOLD = 0.005 # m : écart mini pour corriger un faisceau (0,5 cm)
|
||||
STRIP_ALIGN_CELL = 1.0 # m : maille de comparaison des faisceaux
|
||||
STRIP_ALIGN_MIN_SHARED = 500 # cellules sol communes mini pour valider un offset
|
||||
STRIP_JITTER_BIN = 0.1 # s : durée d'une fenêtre de temps GPS (gigue)
|
||||
STRIP_JITTER_SMOOTH = 5 # fenêtres : largeur de la médiane glissante
|
||||
STRIP_JITTER_MIN_CELLS = 40 # cellules sol partagées mini pour valider une fenêtre
|
||||
STRIP_JITTER_MAX = 0.10 # m : amplitude maxi d'une correction de gigue (garde-fou)
|
||||
|
||||
# Mémo des offsets par fichier : la classification est partagée entre
|
||||
# résolutions, le même LAS sol est rasterisé à 0,5 m puis 0,2 m.
|
||||
_STRIP_OFFSETS_CACHE = {}
|
||||
_STRIP_JITTER_CACHE = {}
|
||||
|
||||
|
||||
def _strip_surface_grid(x, y, z, inv, n_sources, cell=STRIP_ALIGN_CELL):
|
||||
"""Surfaces sol par faisceau sur maille 1 m (moyenne des points bas).
|
||||
|
||||
Pour chaque faisceau et chaque maille : moyenne des points situés à moins
|
||||
de 0,5 m du minimum du faisceau dans la maille (robuste à la végétation
|
||||
résiduelle), mailles à ≥ 2 points seulement.
|
||||
|
||||
Returns:
|
||||
(allcells, grid, key) : cellules triées, grid[faisceau, cellule] = Z
|
||||
moyen (NaN si absent) et clé de maille de CHAQUE point ; None si
|
||||
aucune surface n'est peuplée.
|
||||
"""
|
||||
x0 = np.floor(np.min(x) / cell) * cell
|
||||
y0 = np.floor(np.min(y) / cell) * cell
|
||||
xi = ((x - x0) / cell).astype(np.int64)
|
||||
yi = ((y - y0) / cell).astype(np.int64)
|
||||
ny = int(yi.max()) + 1
|
||||
key = xi * ny + yi
|
||||
|
||||
def _surface(k):
|
||||
m = inv == k
|
||||
kk, zz = key[m], z[m]
|
||||
order = np.argsort(kk, kind='stable')
|
||||
k_s, z_s = kk[order], zz[order]
|
||||
starts = np.flatnonzero(np.r_[True, k_s[1:] != k_s[:-1]])
|
||||
mins = np.minimum.reduceat(z_s, starts)
|
||||
ukey = k_s[starts]
|
||||
thr = mins[np.searchsorted(ukey, kk)]
|
||||
sel = np.flatnonzero(zz <= thr + 0.5)
|
||||
k2, z2 = kk[sel], zz[sel]
|
||||
cnt = np.bincount(k2)
|
||||
sums = np.bincount(k2, weights=z2)
|
||||
v = np.flatnonzero(cnt >= 2)
|
||||
return v, sums[v] / cnt[v]
|
||||
|
||||
surfaces = [_surface(k) for k in range(n_sources)]
|
||||
populated = [c for c, _ in surfaces if len(c)]
|
||||
if not populated:
|
||||
return None
|
||||
allcells = np.unique(np.concatenate(populated))
|
||||
grid = np.full((n_sources, len(allcells)), np.nan)
|
||||
for k, (c, zs) in enumerate(surfaces):
|
||||
grid[k, np.searchsorted(allcells, c)] = zs
|
||||
return allcells, grid, key
|
||||
|
||||
|
||||
def _strip_vertical_offsets(x, y, z, psid, cell=STRIP_ALIGN_CELL,
|
||||
@ -78,37 +138,10 @@ def _strip_vertical_offsets(x, y, z, psid, cell=STRIP_ALIGN_CELL,
|
||||
us, inv = np.unique(psid, return_inverse=True)
|
||||
if len(us) < 2:
|
||||
return {}
|
||||
x0 = np.floor(np.min(x) / cell) * cell
|
||||
y0 = np.floor(np.min(y) / cell) * cell
|
||||
xi = ((x - x0) / cell).astype(np.int64)
|
||||
yi = ((y - y0) / cell).astype(np.int64)
|
||||
ny = int(yi.max()) + 1
|
||||
key = xi * ny + yi
|
||||
|
||||
def _surface(k):
|
||||
m = inv == k
|
||||
kk, zz = key[m], z[m]
|
||||
order = np.argsort(kk, kind='stable')
|
||||
k_s, z_s = kk[order], zz[order]
|
||||
starts = np.flatnonzero(np.r_[True, k_s[1:] != k_s[:-1]])
|
||||
mins = np.minimum.reduceat(z_s, starts)
|
||||
ukey = k_s[starts]
|
||||
thr = mins[np.searchsorted(ukey, kk)]
|
||||
sel = np.flatnonzero(zz <= thr + 0.5)
|
||||
k2, z2 = kk[sel], zz[sel]
|
||||
cnt = np.bincount(k2)
|
||||
sums = np.bincount(k2, weights=z2)
|
||||
v = np.flatnonzero(cnt >= 2)
|
||||
return v, sums[v] / cnt[v]
|
||||
|
||||
surfaces = [_surface(k) for k in range(len(us))]
|
||||
populated = [c for c, _ in surfaces if len(c)]
|
||||
if not populated:
|
||||
built = _strip_surface_grid(x, y, z, inv, len(us), cell)
|
||||
if built is None:
|
||||
return {}
|
||||
allcells = np.unique(np.concatenate(populated))
|
||||
grid = np.full((len(us), len(allcells)), np.nan)
|
||||
for k, (c, zs) in enumerate(surfaces):
|
||||
grid[k, np.searchsorted(allcells, c)] = zs
|
||||
allcells, grid, _key = built
|
||||
comparable = np.sum(~np.isnan(grid), axis=0) >= 2
|
||||
|
||||
offsets = np.zeros(len(us))
|
||||
@ -123,6 +156,177 @@ def _strip_vertical_offsets(x, y, z, psid, cell=STRIP_ALIGN_CELL,
|
||||
for k, p in enumerate(us) if abs(offsets[k]) >= threshold}
|
||||
|
||||
|
||||
def _rolling_median(values, window):
|
||||
"""Médiane glissante (fenêtre tronquée aux bords) sur un petit vecteur."""
|
||||
n = len(values)
|
||||
if window <= 1 or n == 0:
|
||||
return np.array(values, dtype=np.float64)
|
||||
half = window // 2
|
||||
return np.array([np.median(values[max(0, i - half):i + half + 1])
|
||||
for i in range(n)])
|
||||
|
||||
|
||||
def _strip_jitter_offsets(x, y, z, psid, t, cell=STRIP_ALIGN_CELL,
|
||||
bin_seconds=STRIP_JITTER_BIN,
|
||||
smooth=STRIP_JITTER_SMOOTH,
|
||||
min_cells=STRIP_JITTER_MIN_CELLS,
|
||||
max_corr=STRIP_JITTER_MAX,
|
||||
threshold=STRIP_ALIGN_THRESHOLD):
|
||||
"""Mesure la gigue verticale intra-faisceau par fenêtres de temps GPS.
|
||||
|
||||
Le calage constant retire un offset par faisceau, mais le décalage
|
||||
vertical peut aussi varier au fil d'une même passe (vibration capteur,
|
||||
bruit haute fréquence de la trajectoire). Chaque faisceau est découpé en
|
||||
fenêtres de temps GPS ; l'offset robuste de chaque fenêtre est mesuré
|
||||
contre la surface médiane des AUTRES faisceaux (maille 1 m, même
|
||||
sélection de points bas que le calage constant — en recouvrement à deux,
|
||||
une référence incluant le faisceau testé ne révélerait que la moitié du
|
||||
décalage), puis la série est lissée (médiane glissante) pour ne pas
|
||||
suivre le bruit de mesure. z doit déjà être corrigé des offsets constants.
|
||||
|
||||
Args:
|
||||
x, y, z, psid: coordonnées, Z (constant-calé) et PointSourceId des
|
||||
points sol.
|
||||
t: temps GPS (s) de chaque point.
|
||||
cell: taille de maille de comparaison (m).
|
||||
bin_seconds: durée d'une fenêtre de temps (s).
|
||||
smooth: largeur de la médiane glissante (fenêtres).
|
||||
min_cells: cellules partagées minimales pour valider une fenêtre.
|
||||
max_corr: amplitude maximale d'une correction (garde-fou, m).
|
||||
threshold: amplitude de série sous laquelle un faisceau est
|
||||
considéré comme stable (pas de correction).
|
||||
|
||||
Returns:
|
||||
dict {psid: (times, corrections)} des corrections à SOUSTRAIRE,
|
||||
interpolables linéairement au temps GPS de chaque point ; vide si
|
||||
rien de mesurable (faisceau unique, temps absent, pas de
|
||||
recouvrement).
|
||||
"""
|
||||
us, inv = np.unique(psid, return_inverse=True)
|
||||
if len(us) < 2:
|
||||
return {}
|
||||
t = np.asarray(t, dtype=np.float64)
|
||||
if t.size != z.size or not np.isfinite(t).all():
|
||||
return {}
|
||||
built = _strip_surface_grid(x, y, z, inv, len(us), cell)
|
||||
if built is None:
|
||||
return {}
|
||||
allcells, grid, key = built
|
||||
covered = np.sum(~np.isnan(grid), axis=0) >= 2
|
||||
|
||||
# Référence d'un faisceau = médiane des autres faisceaux.
|
||||
ref = np.full(grid.shape, np.nan)
|
||||
if len(us) == 2:
|
||||
ref[0], ref[1] = grid[1], grid[0]
|
||||
else:
|
||||
import warnings
|
||||
with warnings.catch_warnings():
|
||||
warnings.simplefilter("ignore", RuntimeWarning) # tranches tout-NaN
|
||||
for k in range(len(us)):
|
||||
ref[k] = np.nanmedian(np.delete(grid, k, axis=0), axis=0)
|
||||
|
||||
# Résidu vertical de chaque point contre la référence de son faisceau.
|
||||
idx = np.minimum(np.searchsorted(allcells, key), len(allcells) - 1)
|
||||
ref_point = ref[inv, idx]
|
||||
hit = (allcells[idx] == key) & covered[idx] & np.isfinite(ref_point)
|
||||
residual = np.full(np.asarray(z).shape, np.nan)
|
||||
residual[hit] = np.asarray(z)[hit] - ref_point[hit]
|
||||
|
||||
# Origine de temps propre à chaque faisceau : des passes d'une même tuile
|
||||
# peuvent être espacées de plusieurs heures, les fenêtres restent ainsi
|
||||
# dense autour du vol réel (et un saut de semaine GPS ne crée pas de
|
||||
# géantes plages vides).
|
||||
origin = np.full(len(us), np.inf)
|
||||
np.minimum.at(origin, inv, t)
|
||||
tb = np.floor((t - origin[inv]) / bin_seconds).astype(np.int64)
|
||||
n_bins = int(tb.max()) + 1
|
||||
group = inv * n_bins + tb
|
||||
|
||||
# Médiane robuste par groupe (faisceau, fenêtre) : les résidus finis sont
|
||||
# triés en tête de segment, les NaN (sans référence) sont ignorés.
|
||||
order = np.lexsort((np.where(np.isfinite(residual), residual, np.inf), group))
|
||||
g_s, r_s = group[order], residual[order]
|
||||
starts = np.flatnonzero(np.r_[True, g_s[1:] != g_s[:-1]])
|
||||
ends = np.r_[starts[1:], len(g_s)]
|
||||
|
||||
gids, times_c, med = [], [], []
|
||||
for s, e in zip(starts, ends):
|
||||
finite = np.isfinite(r_s[s:e])
|
||||
if int(finite.sum()) < min_cells:
|
||||
continue
|
||||
gids.append(int(g_s[s]))
|
||||
med.append(float(np.median(r_s[s:e][finite])))
|
||||
if not gids:
|
||||
return {}
|
||||
|
||||
gids = np.asarray(gids, dtype=np.int64)
|
||||
med = np.asarray(med, dtype=np.float64)
|
||||
result = {}
|
||||
for k in range(len(us)):
|
||||
selk = gids // n_bins == k
|
||||
if int(selk.sum()) < 3: # série trop courte : correction non fiable
|
||||
continue
|
||||
b = gids[selk] % n_bins
|
||||
cs = np.clip(_rolling_median(med[selk], smooth), -max_corr, max_corr)
|
||||
if float(np.max(np.abs(cs))) < threshold:
|
||||
continue
|
||||
result[int(us[k])] = (origin[k] + (b + 0.5) * bin_seconds, cs)
|
||||
return result
|
||||
|
||||
|
||||
def _apply_strip_jitter(psid, t, jitter):
|
||||
"""Corrections de gigue interpolées au temps GPS de chaque point.
|
||||
|
||||
Args:
|
||||
psid, t: PointSourceId et temps GPS (s) de chaque point.
|
||||
jitter: dict {psid: (times, corrections)} issu de _strip_jitter_offsets.
|
||||
|
||||
Returns:
|
||||
array des corrections à SOUSTRAIRE (0 pour les points sans série).
|
||||
"""
|
||||
corr = np.zeros(len(t))
|
||||
if not jitter:
|
||||
return corr
|
||||
psid = np.asarray(psid)
|
||||
t = np.asarray(t, dtype=np.float64)
|
||||
for p, (times, cs) in jitter.items():
|
||||
m = psid == p
|
||||
if m.any():
|
||||
corr[m] = np.interp(t[m], times, cs)
|
||||
return corr
|
||||
|
||||
|
||||
def _strip_jitter_for_file(las_file, las, offsets):
|
||||
"""Gigue temporelle d'un LAS sol (offsets constants mesurés), mémoïsée."""
|
||||
try:
|
||||
p = Path(las_file)
|
||||
cache_key = (str(p), p.stat().st_mtime_ns)
|
||||
except OSError:
|
||||
cache_key = (str(las_file), 0)
|
||||
if cache_key in _STRIP_JITTER_CACHE:
|
||||
return _STRIP_JITTER_CACHE[cache_key]
|
||||
jitter = {}
|
||||
try:
|
||||
t = np.asarray(las.gps_time, dtype=np.float64)
|
||||
psid = np.asarray(las.point_source_id)
|
||||
except AttributeError:
|
||||
t, psid = None, None # dimensions absentes : pas de gigue mesurable
|
||||
if (t is not None and psid is not None
|
||||
and len(t) == len(las.points) == len(psid)):
|
||||
z = np.asarray(las.z, dtype=np.float64)
|
||||
if offsets:
|
||||
lut = np.zeros(65536)
|
||||
for p_, off in offsets.items():
|
||||
lut[int(p_) & 0xFFFF] = off
|
||||
z = z - lut[np.asarray(las.point_source_id, dtype=np.int64)]
|
||||
jitter = _strip_jitter_offsets(
|
||||
np.asarray(las.x, dtype=np.float64),
|
||||
np.asarray(las.y, dtype=np.float64),
|
||||
z, psid, t)
|
||||
_STRIP_JITTER_CACHE[cache_key] = jitter
|
||||
return jitter
|
||||
|
||||
|
||||
def _strip_offsets_for_file(las_file, las):
|
||||
"""Offsets de calage d'un LAS sol, mémoïsés par (chemin, mtime)."""
|
||||
try:
|
||||
@ -148,18 +352,31 @@ def _strip_offsets_for_file(las_file, las):
|
||||
return offsets
|
||||
|
||||
|
||||
def _write_strip_align_sidecar(dtm_dir, basename, output_suffix, offsets):
|
||||
"""Consigne les offsets de calage appliqués (version et seuil inclus).
|
||||
def _write_strip_align_sidecar(dtm_dir, basename, output_suffix, offsets,
|
||||
jitter=None):
|
||||
"""Consigne les calages appliqués (version, seuil, offsets et gigue).
|
||||
|
||||
Le sidecar sert de suivi de cache : un DTM sans sidecar, ou produit avec
|
||||
une version/un seuil différents, est régénéré. Il est écrit même quand
|
||||
aucun offset n'a été appliqué, pour ne pas re-mesurer une tuile déjà
|
||||
connue comme bien alignée.
|
||||
une version/un seuil/des paramètres de gigue différents, est régénéré. Il
|
||||
est écrit même quand aucune correction n'a été appliquée, pour ne pas
|
||||
re-mesurer une tuile déjà connue comme bien alignée.
|
||||
"""
|
||||
jitter_payload = {}
|
||||
for p, (times, cs) in (jitter or {}).items():
|
||||
cs = np.asarray(cs, dtype=np.float64)
|
||||
jitter_payload[str(p)] = {
|
||||
"bins": int(len(times)),
|
||||
"rms_m": round(float(np.sqrt(np.mean(cs ** 2))), 4),
|
||||
"max_m": round(float(np.max(np.abs(cs))), 4),
|
||||
"series_m": [round(float(c), 4) for c in cs],
|
||||
}
|
||||
payload = {
|
||||
"version": STRIP_ALIGN_VERSION,
|
||||
"threshold": STRIP_ALIGN_THRESHOLD,
|
||||
"offsets": offsets,
|
||||
"jitter_bin": STRIP_JITTER_BIN,
|
||||
"jitter_smooth": STRIP_JITTER_SMOOTH,
|
||||
"jitter": jitter_payload,
|
||||
}
|
||||
try:
|
||||
sidecar = Path(dtm_dir) / f"{basename}_dtm{output_suffix}_stripalign.json"
|
||||
@ -813,73 +1030,139 @@ def _interpolate_holes(dtm, downsample=8):
|
||||
return filled, int(holes.sum())
|
||||
|
||||
|
||||
def _min_return_grid(laz_file, width, height, bounds, chunk_size=2_000_000,
|
||||
strip_offsets=None):
|
||||
"""Rasterize the per-cell minimum z (lowest return) of the full point cloud.
|
||||
# ------------------------------------------------------------
|
||||
# Raccord des bords (edge buffer avec les tuiles adjacentes)
|
||||
# ------------------------------------------------------------
|
||||
# Les visualisations à grand noyau (openness/SVF : rayons jusqu'à 100 m,
|
||||
# LRM : 15 m) tronquent leur fenêtre au bord de la dalle : les rendus
|
||||
# présentent alors une bande d'artefacts à chaque changement de tuile. Avec
|
||||
# edge_buffer > 0, le MNT est rastérisé sur la dalle nominale 1 km ÉTENDUE
|
||||
# d'une bande de `edge_buffer` mètres remplie avec les points sol des 8
|
||||
# tuiles LAZ voisines (lecture PDAL en flux : découpe + filtre de classes).
|
||||
# Les visualisations voient ainsi le terrain réel au-delà du bord, et les
|
||||
# images finales sont recadrées sur la dalle exacte (cf. rendering.py).
|
||||
EDGE_BUFFER_TAG = "LIDAR_EDGE_BUFFER" # tag GeoTIFF : tampon utilisé (m)
|
||||
|
||||
In complex/forested terrain the ground is under-classified, leaving DTM
|
||||
holes. Filling them with the *lowest measured return* of the cell (Wack &
|
||||
Wimmer 2002) recovers a real ground surface (forest floor, rock, clearing)
|
||||
instead of a pure interpolation. The read is streamed in chunks so memory
|
||||
stays bounded to the output grid regardless of the point count.
|
||||
# Décalages des 8 voisines d'une dalle (col, row) en km
|
||||
_NEIGHBOR_OFFSETS = [(-1, -1), (0, -1), (1, -1), (-1, 0),
|
||||
(1, 0), (-1, 1), (0, 1), (1, 1)]
|
||||
|
||||
|
||||
def _tile_coords(name):
|
||||
"""Coordonnées (col, row) en km d'un nom de fichier LHD, ou None."""
|
||||
from .index import parse_basename_coords
|
||||
return parse_basename_coords(Path(name).name)
|
||||
|
||||
|
||||
def _neighbor_laz_files(source_laz):
|
||||
"""Liste les 8 fichiers LAZ/LAS adjacents à `source_laz` dans son dossier."""
|
||||
coords = _tile_coords(source_laz)
|
||||
if coords is None:
|
||||
return []
|
||||
col, row = coords
|
||||
directory = Path(source_laz).parent
|
||||
neighbors = []
|
||||
for dcol, drow in _NEIGHBOR_OFFSETS:
|
||||
nc, nr = col + dcol, row + drow
|
||||
# Noms LHD : col/row en km sur 4 chiffres complétés (ex. 0638_6628) ;
|
||||
# on accepte aussi la variante sans remplissage pour les dalles exotiques.
|
||||
names = {f"LHD_FXX_{nc:04d}_{nr:04d}", f"LHD_FXX_{nc}_{nr}"}
|
||||
matches = []
|
||||
for name in names:
|
||||
matches += list(directory.glob(f"{name}_*.las"))
|
||||
matches += list(directory.glob(f"{name}_*.laz"))
|
||||
if matches:
|
||||
neighbors.append(sorted(matches)[0])
|
||||
else:
|
||||
logger.debug(f" Voisine absente : LHD_FXX_{nc:04d}_{nr:04d} (bande de bord vide)")
|
||||
return neighbors
|
||||
|
||||
|
||||
def _neighbor_ground_points(source_laz, bounds, classes):
|
||||
"""Points sol des tuiles voisines dans `bounds` (raccord des bords).
|
||||
|
||||
Lecture PDAL en flux par voisine : découpe sur l'emprise étendue puis
|
||||
filtre de classes (mêmes codes que le MNT). Best-effort : une voisine
|
||||
illisible ou absente est ignorée — la bande correspondante reste vide.
|
||||
|
||||
Args:
|
||||
laz_file: Path to the full (unclassified) LAZ/LAS file.
|
||||
width, height: Output grid dimensions (pixels).
|
||||
bounds: (min_x, min_y, max_x, max_y) the grid covers.
|
||||
chunk_size: Points per streaming chunk.
|
||||
strip_offsets: Optionnel : dict {point_source_id: offset} issu du
|
||||
calage des faisceaux, retranché aux Z du nuage complet pour
|
||||
rester cohérent avec le MNT calé.
|
||||
source_laz: LAZ de la tuile traitée (repère pour trouver les voisines).
|
||||
bounds: (min_x, min_y, max_x, max_y) de l'emprise étendue.
|
||||
classes: codes LAS à extraire (ex. [2] = sol).
|
||||
|
||||
Returns:
|
||||
(height, width) float32 array of per-cell min z (NaN where no point).
|
||||
(xs, ys, zs) concaténés (tableaux vides si aucune voisine).
|
||||
"""
|
||||
import laspy
|
||||
import tempfile
|
||||
|
||||
neighbors = _neighbor_laz_files(source_laz)
|
||||
if not neighbors:
|
||||
return (np.empty(0),) * 3
|
||||
|
||||
min_x, min_y, max_x, max_y = bounds
|
||||
grid = np.full((height, width), np.nan, dtype=np.float32)
|
||||
rng = [[min_x, max_x], [min_y, max_y]]
|
||||
codes = sorted(set(int(c) for c in classes)) or [2]
|
||||
limits = ",".join(f"Classification[{c}:{c}]" for c in codes)
|
||||
|
||||
lut = None
|
||||
if strip_offsets:
|
||||
lut = np.zeros(65536)
|
||||
for p, off in strip_offsets.items():
|
||||
lut[int(p) & 0xFFFF] = off
|
||||
xs, ys, zs = [], [], []
|
||||
found = 0
|
||||
for neighbor in neighbors:
|
||||
tmp_path = None
|
||||
try:
|
||||
with tempfile.NamedTemporaryFile(suffix='.las', delete=False) as tmp:
|
||||
tmp_path = tmp.name
|
||||
pipeline = json.dumps({
|
||||
"pipeline": [
|
||||
{"type": "readers.las", "filename": str(neighbor)},
|
||||
{"type": "filters.crop",
|
||||
"bounds": f"([{min_x},{max_x}],[{min_y},{max_y}])"},
|
||||
{"type": "filters.range", "limits": limits},
|
||||
{"type": "writers.las", "filename": tmp_path},
|
||||
]
|
||||
})
|
||||
result = subprocess.run(
|
||||
["pdal", "pipeline", "--stdin"],
|
||||
input=pipeline, capture_output=True, text=True, timeout=300
|
||||
)
|
||||
if result.returncode != 0:
|
||||
raise RuntimeError(result.stderr[:200])
|
||||
import laspy
|
||||
las = laspy.read(tmp_path)
|
||||
if len(las.points) > 0:
|
||||
xs.append(np.asarray(las.x, dtype=np.float64))
|
||||
ys.append(np.asarray(las.y, dtype=np.float64))
|
||||
zs.append(np.asarray(las.z, dtype=np.float64))
|
||||
found += 1
|
||||
except Exception as e:
|
||||
logger.debug(f" Voisine {Path(neighbor).name} ignorée : {e}")
|
||||
finally:
|
||||
if tmp_path:
|
||||
try:
|
||||
Path(tmp_path).unlink(missing_ok=True)
|
||||
except Exception:
|
||||
pass
|
||||
|
||||
def process(points):
|
||||
if len(points) == 0:
|
||||
return
|
||||
x = np.asarray(points.x, dtype=np.float64)
|
||||
y = np.asarray(points.y, dtype=np.float64)
|
||||
z = np.asarray(points.z, dtype=np.float64)
|
||||
if lut is not None:
|
||||
try:
|
||||
z = z - lut[np.asarray(points.point_source_id, dtype=np.int64)]
|
||||
except AttributeError:
|
||||
pass
|
||||
st = binned_statistic_2d(x, y, z, statistic='min',
|
||||
bins=[width, height], range=rng)
|
||||
# Match the DTM convention: .T then flip Y (north at top).
|
||||
cell_min = st.statistic.T[::-1, :].astype(np.float32)
|
||||
# fmin ignores NaN so cells without a point in this chunk stay NaN.
|
||||
np.fmin(grid, cell_min, out=grid)
|
||||
if not found:
|
||||
logger.warning(" Aucune voisine lisible — bande de bord vide")
|
||||
return (np.empty(0),) * 3
|
||||
logger.info(f" Raccord bords : {found} voisine(s), "
|
||||
f"{sum(len(a) for a in xs):,} pts sol")
|
||||
return np.concatenate(xs), np.concatenate(ys), np.concatenate(zs)
|
||||
|
||||
|
||||
def read_dtm_edge_buffer(dtm_path):
|
||||
"""Tampon de raccord enregistré dans un DTM (m ; 0 si absent/illisible)."""
|
||||
try:
|
||||
with laspy.open(str(laz_file)) as las:
|
||||
for chunk in las.chunk_iterator(chunk_size):
|
||||
process(chunk)
|
||||
except Exception as e:
|
||||
logger.warning(f" Lecture streaming impossible ({e}) — lecture complète")
|
||||
las = _read_with_pdal(laz_file)
|
||||
if las is None:
|
||||
return grid
|
||||
process(las)
|
||||
return grid
|
||||
with rasterio.open(dtm_path) as src:
|
||||
tags = src.tags()
|
||||
return float(tags.get(EDGE_BUFFER_TAG, 0.0) or 0.0)
|
||||
except Exception:
|
||||
return 0.0
|
||||
|
||||
|
||||
def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False,
|
||||
output_suffix="", source_laz=None, bare_earth=False,
|
||||
pure=False, strip_align=True):
|
||||
output_suffix="", source_laz=None,
|
||||
pure=False, strip_align=True, edge_buffer=0.0,
|
||||
neighbor_classes=None):
|
||||
"""Create DTM using fast binning method with gap filling.
|
||||
|
||||
Args:
|
||||
@ -889,19 +1172,26 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False,
|
||||
resolution: Grid resolution in meters per pixel.
|
||||
force: If True, regenerate even if DTM already exists.
|
||||
output_suffix: Suffix for output filename (e.g. '_r0p2' for additional resolutions).
|
||||
source_laz: Optionnel : chemin du LAZ complet (non classé). Utilisé
|
||||
uniquement avec bare_earth (plancher au retour le plus bas).
|
||||
bare_earth: If True, pull the DTM down to the lowest measured return of
|
||||
each cell (bare-earth floor). This requalifies the lowest point of
|
||||
every column as terrain, recovering the ground under dense
|
||||
vegetation / steep relief that the ground classifier rejected.
|
||||
source_laz: Optionnel : LAZ complet de la tuile traitée (repère pour
|
||||
lire les 8 tuiles voisines du raccord de bords, cf. edge_buffer).
|
||||
pure: Sans effet (conservé pour compatibilité). Fonctionnement
|
||||
historique rétabli : petits trous comblés par fillnodata, grands
|
||||
trous laissés en nodata (rendus en noir dans les rendus).
|
||||
strip_align: Si True (défaut), mesure et corrige les écarts verticaux
|
||||
entre faisceaux de vol (PointSourceId) avant rastérisation ; les
|
||||
offsets ≥ STRIP_ALIGN_THRESHOLD (0,5 cm) sont consignés dans un
|
||||
sidecar *_dtm*_stripalign.json.
|
||||
entre faisceaux de vol (PointSourceId) avant rastérisation : offsets
|
||||
constants ≥ STRIP_ALIGN_THRESHOLD (0,5 cm) puis gigue intra-faisceau
|
||||
par fenêtres de temps GPS ; tout est consigné dans un sidecar
|
||||
*_dtm*_stripalign.json.
|
||||
edge_buffer: Raccord des bords en mètres (0 = désactivé). Le MNT couvre
|
||||
alors la dalle nominale 1 km étendue de cette bande, remplie avec
|
||||
les points sol des 8 tuiles LAZ voisines (source_laz requis) ; les
|
||||
visualisations calculent sur l'emprise étendue puis les images sont
|
||||
recadrées sur la dalle exacte (rendering.py). Le tampon est inscrit
|
||||
dans le tag GeoTIFF LIDAR_EDGE_BUFFER pour l'invalidation du cache.
|
||||
neighbor_classes: codes LAS extraits chez les voisines (défaut : les
|
||||
classes IGN du MNT, ex. [2]). Les voisines sont lues dans leur
|
||||
pré-classification fournisseur, même si la tuile centrale est
|
||||
classée SMRF/CSF (bande de contexte, quelques cm d'écart au pire).
|
||||
|
||||
Returns:
|
||||
Path to output DTM GeoTIFF, or None on failure.
|
||||
@ -934,6 +1224,8 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False,
|
||||
# Calage vertical des faisceaux de vol avant rastérisation : best-effort,
|
||||
# en cas d'échec de la mesure on continue non calé (jamais d'abort).
|
||||
strip_offsets = {}
|
||||
strip_jitter = {}
|
||||
gps_time = None
|
||||
if strip_align:
|
||||
try:
|
||||
strip_offsets = _strip_offsets_for_file(las_file, las)
|
||||
@ -946,6 +1238,24 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False,
|
||||
else:
|
||||
logger.debug(" Calage faisceaux : aucun écart >= "
|
||||
f"{STRIP_ALIGN_THRESHOLD * 100:.1f} cm, rien à corriger")
|
||||
# 2ᵉ passe : gigue intra-faisceau (temps GPS requis, ignorée sinon).
|
||||
try:
|
||||
gps_time = np.asarray(las.gps_time, dtype=np.float64)
|
||||
except AttributeError:
|
||||
gps_time = None
|
||||
if gps_time is not None and len(gps_time) == len(las.points):
|
||||
try:
|
||||
strip_jitter = _strip_jitter_for_file(las_file, las, strip_offsets)
|
||||
except Exception as e:
|
||||
logger.warning(f" Mesure de la gigue intra-faisceau impossible ({e}) — gigue non corrigée")
|
||||
strip_jitter = {}
|
||||
if strip_jitter:
|
||||
for p_, (jt, jc) in sorted(strip_jitter.items()):
|
||||
logger.info(f" Gigue PSID {p_} : ±{np.max(np.abs(jc)) * 100:.1f} cm "
|
||||
f"(rms {np.sqrt(np.mean(np.asarray(jc) ** 2)) * 100:.1f} cm, "
|
||||
f"{len(jt)} fenêtres de {STRIP_JITTER_BIN:g} s)")
|
||||
else:
|
||||
logger.debug(" Gigue intra-faisceau : rien à corriger")
|
||||
|
||||
try:
|
||||
|
||||
@ -955,6 +1265,32 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False,
|
||||
width = int(np.ceil((max_x - min_x) / resolution))
|
||||
height = int(np.ceil((max_y - min_y) / resolution))
|
||||
|
||||
# Raccord des bords : emprise = dalle nominale 1 km (alignée sur la
|
||||
# grille multi-tuiles) + bande de edge_buffer mètres remplie par les
|
||||
# points sol des voisines. Sinon : bornes de l'en-tête (historique).
|
||||
used_edge_buffer = 0.0
|
||||
if edge_buffer > 0:
|
||||
coords = _tile_coords(source_laz or las_file)
|
||||
if coords is not None:
|
||||
buffer_px = max(1, int(round(edge_buffer / resolution)))
|
||||
buffer_m = buffer_px * resolution
|
||||
col_km, row_km = coords
|
||||
# Grille LHD : (col, row) = coin nord-ouest en km →
|
||||
# X ∈ [col, col+1] km, Y ∈ [row-1, row] km (bord nord = row).
|
||||
min_x = float(col_km) * 1000.0
|
||||
max_x = min_x + 1000.0
|
||||
max_y = float(row_km) * 1000.0
|
||||
min_y = max_y - 1000.0
|
||||
ext_bounds = (min_x - buffer_m, min_y - buffer_m,
|
||||
max_x + buffer_m, max_y + buffer_m)
|
||||
width = int(round(1000.0 / resolution)) + 2 * buffer_px
|
||||
height = width
|
||||
min_x, min_y, max_x, max_y = ext_bounds
|
||||
used_edge_buffer = float(edge_buffer)
|
||||
else:
|
||||
logger.warning(" Raccord des bords impossible : nom de "
|
||||
f"fichier non LHD ({basename}) — tuile seule")
|
||||
|
||||
logger.debug(f" Bounds: X[{min_x:.1f}, {max_x:.1f}] Y[{min_y:.1f}, {max_y:.1f}]")
|
||||
logger.debug(f" Grid: {width}x{height} pixels ({len(las.points):,} points)")
|
||||
logger.info(f" Rasterisation {width}x{height} ({len(las.points):,} points)...")
|
||||
@ -967,6 +1303,20 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False,
|
||||
for p, off in strip_offsets.items():
|
||||
lut[int(p) & 0xFFFF] = off
|
||||
zs = zs - lut[np.asarray(las.point_source_id, dtype=np.int64)]
|
||||
if strip_jitter and gps_time is not None:
|
||||
zs = zs - _apply_strip_jitter(las.point_source_id, gps_time, strip_jitter)
|
||||
|
||||
# Points sol des tuiles voisines dans la bande de raccord (best-effort,
|
||||
# non calés par faisceau : bande de contexte, l'image finale est
|
||||
# recadrée sur la dalle avant livraison).
|
||||
if used_edge_buffer > 0 and source_laz is not None:
|
||||
nx, ny, nz = _neighbor_ground_points(
|
||||
source_laz, (min_x, min_y, max_x, max_y),
|
||||
neighbor_classes if neighbor_classes is not None else [2])
|
||||
if len(nx):
|
||||
xs = np.concatenate([xs, nx])
|
||||
ys = np.concatenate([ys, ny])
|
||||
zs = np.concatenate([zs, nz])
|
||||
|
||||
stat = binned_statistic_2d(
|
||||
xs, ys, zs,
|
||||
@ -981,20 +1331,9 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False,
|
||||
# Comblement « historique » (fonctionnement d'origine, rétabli) :
|
||||
# seuls les petits trous proches des données sont remplis ; les grands
|
||||
# trous restent en nodata et apparaissent en noir dans les rendus.
|
||||
# Le plancher au retour le plus bas n'est appliqué qu'à la demande
|
||||
# explicite (--bare-earth).
|
||||
if bare_earth and source_laz is not None:
|
||||
min_grid = _min_return_grid(source_laz, width, height,
|
||||
(min_x, min_y, max_x, max_y),
|
||||
strip_offsets=strip_offsets)
|
||||
# Cellules sans sol mesuré (NaN) : le retour le plus bas devient
|
||||
# la mesure — sinon la comparaison NaN est fausse et la cellule
|
||||
# retombe sur l'interpolation en fin de passe.
|
||||
has_min = ~np.isnan(min_grid)
|
||||
lower = has_min & (min_grid < dtm)
|
||||
lower |= has_min & np.isnan(dtm)
|
||||
dtm = np.where(lower, min_grid, dtm)
|
||||
logger.info(f" Sol nu : {int(lower.sum()):,} cellules raménées au retour le plus bas")
|
||||
# Volontairement PAS de plancher au retour le plus bas : sous canopée
|
||||
# dense ce retour est la végétation, qui imprimerait les arbres dans
|
||||
# le MNT.
|
||||
|
||||
# Fill small gaps (< 1 m from data) precisely — comme avant
|
||||
nan_count = np.count_nonzero(np.isnan(dtm))
|
||||
@ -1026,10 +1365,14 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False,
|
||||
compress='lzw'
|
||||
) as dst:
|
||||
dst.write(dtm.astype('float32'), 1)
|
||||
if used_edge_buffer > 0:
|
||||
# Tampon de raccord inscrit dans le fichier : changement de
|
||||
# --edge-buffer ⇒ invalidation automatique du cache DTM.
|
||||
dst.update_tags(**{EDGE_BUFFER_TAG: f"{used_edge_buffer:g}"})
|
||||
|
||||
if strip_align:
|
||||
_write_strip_align_sidecar(dtm_dir, basename, output_suffix,
|
||||
strip_offsets)
|
||||
strip_offsets, strip_jitter)
|
||||
logger.info(f" ✓ DTM créé: {output_tif.name}")
|
||||
return output_tif
|
||||
|
||||
|
||||
Reference in New Issue
Block a user