Borner le comblement du MNT à l'enveloppe des points, ajouter la couche précision et un affichage relief/précision, encoder les sous-tuiles une seule fois en q75
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
Antoine Jacquin
@ -1446,6 +1446,200 @@ def _repair_laz_with_laspy(input_laz, output_las):
|
||||
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.
|
||||
|
||||
@ -1834,6 +2028,13 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False,
|
||||
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))
|
||||
@ -1861,24 +2062,18 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False,
|
||||
# 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
|
||||
# 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}%)")
|
||||
|
||||
max_gap_pixels = max(1, int(1.0 / resolution))
|
||||
t_fill = time.perf_counter()
|
||||
from rasterio.fill import fillnodata
|
||||
valid_mask = ~np.isnan(dtm)
|
||||
dtm_filled = fillnodata(dtm, mask=valid_mask, max_search_distance=max_gap_pixels)
|
||||
small_gap_mask = np.isnan(dtm) & ~np.isnan(dtm_filled)
|
||||
filled_count = np.count_nonzero(small_gap_mask)
|
||||
if filled_count > 0:
|
||||
dtm = np.where(small_gap_mask, dtm_filled, dtm)
|
||||
logger.info(f" {filled_count:,} petits trous comblés "
|
||||
f"(< {max_gap_pixels}px, {time.perf_counter() - t_fill:.1f}s)")
|
||||
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"
|
||||
@ -1894,6 +2089,9 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False,
|
||||
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.
|
||||
|
||||
Reference in New Issue
Block a user