Calculer le sidecar qualité après le DTM et ajouter --quality-backfill

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>

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).

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

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

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()