41206f65de
Trois chantiers liés à la qualité et au coût des rendus : - Calage vertical des faisceaux : les passes d'une tuile peuvent être biaisées de quelques cm (±2,5 cm mesurés sur 1000_6882), créant des marches aux recouvrements. Les offsets par PointSourceId sont mesurés sur les points sol de la tuile (réf. médiane itérée) et retranchés ≥ 0,5 cm avant rastérisation, avec sidecar de cache et application au plancher bare-earth. - Openness décimée ×2 : lancé de rayons sur grille par blocs (max/min) puis rééchantillonnage bilinéaire — 532 s → 40 s par tuile à 0,2 m sur CPU, signal archéologique préservé. Réglable --openness-downsample. - Génération webapp concentrée sur les couches affichées : les défauts /api/generate, /api/preview et le sélecteur génèrent le panneau complet (slope, aspect, pos_open) au lieu d'aspect seul.
1038 lines
40 KiB
Python
1038 lines
40 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.
|
||
STRIP_ALIGN_VERSION = 1
|
||
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
|
||
|
||
# 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 = {}
|
||
|
||
|
||
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 {}
|
||
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:
|
||
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
|
||
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 _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):
|
||
"""Consigne les offsets de calage appliqués (version et seuil inclus).
|
||
|
||
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.
|
||
"""
|
||
payload = {
|
||
"version": STRIP_ALIGN_VERSION,
|
||
"threshold": STRIP_ALIGN_THRESHOLD,
|
||
"offsets": offsets,
|
||
}
|
||
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())
|
||
|
||
|
||
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.
|
||
|
||
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.
|
||
|
||
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é.
|
||
|
||
Returns:
|
||
(height, width) float32 array of per-cell min z (NaN where no point).
|
||
"""
|
||
import laspy
|
||
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]]
|
||
|
||
lut = None
|
||
if strip_offsets:
|
||
lut = np.zeros(65536)
|
||
for p, off in strip_offsets.items():
|
||
lut[int(p) & 0xFFFF] = off
|
||
|
||
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)
|
||
|
||
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
|
||
|
||
|
||
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):
|
||
"""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 : 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.
|
||
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.
|
||
|
||
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 = {}
|
||
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")
|
||
|
||
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))
|
||
|
||
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)]
|
||
|
||
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.
|
||
# 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")
|
||
|
||
# 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 strip_align:
|
||
_write_strip_align_sidecar(dtm_dir, basename, output_suffix,
|
||
strip_offsets)
|
||
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 |