From 7d5b5ecd4292a5a4bed5360369460078069ffc20 Mon Sep 17 00:00:00 2001 From: Antoine Jacquin Date: Sun, 27 Sep 2026 15:17:32 +0200 Subject: [PATCH] =?UTF-8?q?Calculer=20le=20sidecar=20qualit=C3=A9=20apr?= =?UTF-8?q?=C3=A8s=20le=20DTM=20et=20ajouter=20--quality-backfill?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Sonnet 5 --- lidar_pipeline/cli.py | 17 ++++++ lidar_pipeline/pipeline.py | 6 ++ lidar_pipeline/quality.py | 91 ++++++++++++++++++++++++++++ lidar_pipeline/tests/test_quality.py | 86 ++++++++++++++++++++++++++ 4 files changed, 200 insertions(+) diff --git a/lidar_pipeline/cli.py b/lidar_pipeline/cli.py index 7ccfac2..d06c289 100644 --- a/lidar_pipeline/cli.py +++ b/lidar_pipeline/cli.py @@ -273,6 +273,13 @@ def main(): help="Conservé pour compatibilité : le catalogue est désormais régénéré après " "chaque tuile terminée par défaut (sauf --no-index)" ) + parser.add_argument( + "--quality-backfill", + action="store_true", + help="Calculer uniquement les sidecars qualité (densité sol, dates de vol) " + "des LAZ du dossier d'entrée qui n'en ont pas, puis régénérer l'inventaire " + "(sans DTM ni visualisation, ~5 s par dalle)" + ) args = parser.parse_args() @@ -325,6 +332,16 @@ def main(): logger.warning("Aucune tuile traitée trouvée — catalogue non généré") return + if args.quality_backfill: + from .dtm import parse_ign_classes + from .index import build_index + from .quality import backfill_quality + n = backfill_quality(args.input, args.output, + codes=tuple(parse_ign_classes(args.ign_classes))) + logger.info(f"{n} sidecar(s) qualité écrit(s)") + build_index(args.output, args.format) + return + # Nouveau run : remise à zéro du journal d'événements (file de # génération). Ici, avant le téléchargement — les événements download # survivent au démarrage du traitement (process_all ne tronque plus). diff --git a/lidar_pipeline/pipeline.py b/lidar_pipeline/pipeline.py index cd50907..43f9777 100644 --- a/lidar_pipeline/pipeline.py +++ b/lidar_pipeline/pipeline.py @@ -649,6 +649,12 @@ class LidarArchaeoPipeline: dtm_rebuilt = True self._write_dtm_method(basename, res_suffix) + if i == 0: + # Sidecar qualité (densité sol, dates de vol) pour l'encart de + # l'export PDF : calculé une fois, même si le DTM vient du cache. + from .quality import ensure_quality + ensure_quality(las_file, basename, self.output_dir) + # Process each resolution: visualizations + PDF # Option de calcul (surcharge le défaut du module) appliquée ICI car les # workers (spawn) réimportent les modules à froid : c'est le seul endroit diff --git a/lidar_pipeline/quality.py b/lidar_pipeline/quality.py index 1922ce6..8e5ea5a 100644 --- a/lidar_pipeline/quality.py +++ b/lidar_pipeline/quality.py @@ -148,3 +148,94 @@ def load_quality_table(output_dir): if data is not None: table[f.stem] = data return table + + +_LAZ_EXTS = (".copc.laz", ".copc.las", ".laz", ".las") + + +def _laz_basename(path): + """Nom de dalle sans extension LiDAR (.copc.laz compris).""" + name = Path(path).name + for ext in _LAZ_EXTS: + if name.lower().endswith(ext): + return name[:-len(ext)] + return Path(path).stem + + +def _cell_bounds_of(basename): + """Emprise nominale 1 km d'une dalle LHD, ou None hors motif LHD.""" + from .index import parse_basename_coords + coords = parse_basename_coords(basename) + if coords is None: + return None + col, row = coords + return (col * 1000.0, (row - 1) * 1000.0, (col + 1) * 1000.0, row * 1000.0) + + +def quality_from_las(las_path, bounds, codes=None): + """Qualité mesurée depuis un LAS/LAZ ; None en cas d'échec (jamais d'exception). + + codes : classes à retenir (retours ≥ 1) ; None = fichier déjà filtré sol. + """ + try: + import laspy + import numpy as np + las = laspy.read(str(las_path)) + keep = np.ones(len(las.points), dtype=bool) + if codes is not None: + keep = ((np.asarray(las.return_number) >= 1) + & (np.asarray(las.number_of_returns) >= 1) + & np.isin(np.asarray(las.classification), np.asarray(sorted(codes)))) + try: + gps = np.asarray(las.gps_time, dtype=np.float64)[keep] + except AttributeError: + gps = None + adjusted = las.header.global_encoding.gps_time_type == laspy.header.GpsTimeType.STANDARD + created = las.header.creation_date + return compute_quality(np.asarray(las.x)[keep], np.asarray(las.y)[keep], gps, + bounds, gps_adjusted=adjusted, + header_date=created.isoformat() if created else None) + except Exception as e: # noqa: BLE001 — qualité best-effort + logger.warning(f" Mesure de qualité impossible ({Path(las_path).name}) : {e}") + return None + + +def ensure_quality(las_path, basename, output_dir, codes=None): + """Garantit le sidecar qualité d'une dalle (calcul si absent ou périmé). + + Returns: + True si un sidecar valide existe à la sortie. + """ + if read_quality(output_dir, basename) is not None: + return True + bounds = _cell_bounds_of(basename) + if bounds is None: + return False + data = quality_from_las(las_path, bounds, codes) + if data is None: + return False + try: + write_quality(output_dir, basename, data) + except OSError as e: + logger.warning(f" Écriture du sidecar qualité impossible : {e}") + return False + return True + + +def backfill_quality(input_dir, output_dir, codes=(2,)): + """Rattrapage : sidecar qualité de chaque LAZ LHD d'input_dir qui n'en a pas. + + Returns: + Nombre de sidecars écrits. + """ + files = sorted(p for p in Path(input_dir).iterdir() + if p.is_file() and p.name.lower().endswith(_LAZ_EXTS)) + written = 0 + for i, laz in enumerate(files, 1): + base = _laz_basename(laz) + if read_quality(output_dir, base) is not None or _cell_bounds_of(base) is None: + continue + logger.info(f"Qualité [{i}/{len(files)}] {base}") + if ensure_quality(laz, base, output_dir, codes): + written += 1 + return written diff --git a/lidar_pipeline/tests/test_quality.py b/lidar_pipeline/tests/test_quality.py index 0403778..c46d4f9 100644 --- a/lidar_pipeline/tests/test_quality.py +++ b/lidar_pipeline/tests/test_quality.py @@ -89,3 +89,89 @@ def test_quality_module_imports_without_numpy(): "import lidar_pipeline.quality as q; print(q.QUALITY_VERSION)") r = subprocess.run([sys.executable, "-c", code], capture_output=True, text=True) assert r.returncode == 0, r.stderr + + +def _write_las(path, x, y, cls, t, adjusted=True, creation=None): + """LAS minimal (format 6 : gps_time) pour les tests.""" + import laspy + import numpy as np + header = laspy.LasHeader(point_format=6, version="1.4") + header.scales = [0.01, 0.01, 0.01] + header.offsets = [652000.0, 6861000.0, 0.0] + if adjusted: + header.global_encoding.gps_time_type = laspy.header.GpsTimeType.STANDARD + if creation is not None: + header.creation_date = creation + las = laspy.LasData(header) + las.x = np.asarray(x, float); las.y = np.asarray(y, float) + las.z = np.zeros(len(x)) + las.classification = np.asarray(cls, np.uint8) + las.return_number = np.ones(len(x), np.uint8) + las.number_of_returns = np.ones(len(x), np.uint8) + las.gps_time = np.asarray(t, float) + las.write(str(path)) + + +def test_quality_from_las_filters_classes(tmp_path): + from lidar_pipeline.quality import quality_from_las + p = tmp_path / "LHD_FXX_0652_6862_PTS_LAMB93_IGN69.laz" + t = _gps_seconds(2022, 6, 1) + _write_las(p, [652100, 652200, 652300], [6861100, 6861200, 6861300], + [2, 2, 5], [t, t, t]) + q = quality_from_las(p, BOUNDS, codes=(2,)) + assert abs(q["ground_density"] - 2e-6) < 1e-9 + assert q["acq_start"] == "2022-06-01" and q["acq_source"] == "gps" + + +def test_quality_from_las_week_time_uses_header_date(tmp_path): + from datetime import date + from lidar_pipeline.quality import quality_from_las + p = tmp_path / "x.las" + _write_las(p, [652100], [6861100], [2], [1000.0], adjusted=False, + creation=date(2021, 9, 2)) + q = quality_from_las(p, BOUNDS) + assert q["acq_source"] == "header" and q["acq_start"] == "2021-09-02" + + +def test_quality_from_las_unreadable_returns_none(tmp_path): + from lidar_pipeline.quality import quality_from_las + p = tmp_path / "bad.laz"; p.write_bytes(b"pas un LAS") + assert quality_from_las(p, BOUNDS) is None + + +def test_ensure_quality_skips_existing_and_unparsable(tmp_path): + from lidar_pipeline.quality import ensure_quality, read_quality, write_quality + base = "LHD_FXX_0652_6862_PTS_LAMB93_IGN69" + p = tmp_path / f"{base}_ground.las" + _write_las(p, [652100], [6861100], [2], [_gps_seconds(2022, 6, 1)]) + assert ensure_quality(p, base, tmp_path) is True + assert read_quality(tmp_path, base)["acq_start"] == "2022-06-01" + write_quality(tmp_path, base, dict(read_quality(tmp_path, base), ground_density=99.0)) + assert ensure_quality(p, base, tmp_path) is True + assert read_quality(tmp_path, base)["ground_density"] == 99.0 # non recalculé + assert ensure_quality(p, "nom_libre", tmp_path) is False + + +def test_backfill_quality_reads_input_laz(tmp_path): + from lidar_pipeline.quality import backfill_quality, read_quality + inp = tmp_path / "input"; inp.mkdir() + out = tmp_path / "output" + base = "LHD_FXX_0652_6862_PTS_LAMB93_IGN69" + _write_las(inp / f"{base}.copc.laz", [652100, 652200], [6861100, 6861200], + [2, 6], [_gps_seconds(2022, 6, 1)] * 2) + assert backfill_quality(inp, out) == 1 + assert abs(read_quality(out, base)["ground_density"] - 1e-6) < 1e-9 + assert backfill_quality(inp, out) == 0 # déjà à jour + + +def test_cli_quality_backfill(tmp_path): + import subprocess, sys + inp = tmp_path / "input"; inp.mkdir() + out = tmp_path / "output" + base = "LHD_FXX_0652_6862_PTS_LAMB93_IGN69" + _write_las(inp / f"{base}.laz", [652100], [6861100], [2], [_gps_seconds(2022, 6, 1)]) + r = subprocess.run([sys.executable, "-m", "lidar_pipeline", str(inp), "-o", str(out), + "--quality-backfill"], capture_output=True, text=True, timeout=180) + assert r.returncode == 0, r.stderr + assert (out / "quality" / f"{base}.json").is_file() + assert not (out / "DTM").exists()