Ajouter la mesure de qualité par dalle (densité sol, date d'acquisition) et ses sidecars

Co-Authored-By: Claude Haiku 4.5 <noreply@anthropic.com>

lidar_pipeline/quality.py

@@ -0,0 +1,150 @@
+"""Qualité des données LiDAR par dalle : densité de points sol et date d'acquisition.
+
+Un sidecar JSON par dalle (output/quality/{basename}.json) est écrit par le
+pipeline, recopié dans l'inventaire (index_tiles.json, table `quality`) et
+conservé sur les machines légères : l'encart qualité de l'export PDF le lit
+depuis le disque local, amont joignable ou non.
+
+numpy n'est importé que par le calcul : l'image légère (sans numpy) lit et
+écrit les sidecars.
+"""
+
+import json
+import logging
+import os
+import tempfile
+from datetime import datetime, timedelta, timezone
+from pathlib import Path
+
+logger = logging.getLogger("lidar")
+
+QUALITY_VERSION = 1
+QUALITY_DIRNAME = "quality"
+DENSITY_CELL_M = 50.0      # maille de la grille de densité (20 × 20 par dalle)
+EMPTY_PIXEL_M = 0.2        # pixel du MNT : part sans point sol ≈ part interpolée
+
+_GPS_EPOCH = datetime(1980, 1, 6, tzinfo=timezone.utc)
+_GPS_ADJUSTED_OFFSET = 1e9  # temps GPS ajusté standard = secondes GPS − 1e9
+
+
+def gps_adjusted_to_date(t):
+    """Temps GPS ajusté standard → date ISO (UTC ; secondes intercalaires
+    négligées : précision au jour)."""
+    return (_GPS_EPOCH + timedelta(seconds=float(t) + _GPS_ADJUSTED_OFFSET)).date().isoformat()
+
+
+def compute_quality(x, y, gps_time, bounds, gps_adjusted=True, header_date=None):
+    """Mesure la qualité d'une dalle à partir de ses points sol.
+
+    Args:
+        x, y: coordonnées L93 des points sol.
+        gps_time: temps GPS des points (ou None).
+        bounds: emprise nominale (min_x, min_y, max_x, max_y) ; les points
+            hors emprise (bande de raccord) sont ignorés.
+        gps_adjusted: True si gps_time est en temps GPS ajusté standard.
+        header_date: date ISO de l'en-tête LAS (repli des dates).
+
+    Returns:
+        dict JSON-sérialisable (cf. spec, section 1).
+    """
+    import numpy as np
+
+    min_x, min_y, max_x, max_y = bounds
+    x = np.asarray(x, dtype=np.float64)
+    y = np.asarray(y, dtype=np.float64)
+    inside = (x >= min_x) & (x < max_x) & (y >= min_y) & (y < max_y)
+    t = None
+    if gps_time is not None and gps_adjusted:
+        t = np.asarray(gps_time, dtype=np.float64)
+        t = t[inside] if len(t) == len(x) else None
+    x, y = x[inside], y[inside]
+    n = len(x)
+    width, height = max_x - min_x, max_y - min_y
+
+    nx = int(round(width / DENSITY_CELL_M))
+    ny = int(round(height / DENSITY_CELL_M))
+    cx = np.clip(((x - min_x) / DENSITY_CELL_M).astype(np.int64), 0, nx - 1)
+    cy = np.clip(((max_y - y) / DENSITY_CELL_M).astype(np.int64), 0, ny - 1)
+    counts = np.bincount(cy * nx + cx, minlength=nx * ny).reshape(ny, nx)
+    grid = np.round(counts / (DENSITY_CELL_M ** 2), 2)
+
+    px = int(round(width / EMPTY_PIXEL_M))
+    py = int(round(height / EMPTY_PIXEL_M))
+    occupied = np.zeros(px * py, dtype=bool)
+    if n:
+        ix = np.clip(((x - min_x) / EMPTY_PIXEL_M).astype(np.int64), 0, px - 1)
+        iy = np.clip(((max_y - y) / EMPTY_PIXEL_M).astype(np.int64), 0, py - 1)
+        occupied[iy * px + ix] = True
+    empty_fraction = 1.0 - float(occupied.sum()) / (px * py)
+
+    if t is not None and len(t):
+        start, end = gps_adjusted_to_date(t.min()), gps_adjusted_to_date(t.max())
+        source = "gps"
+    elif header_date:
+        start = end = header_date
+        source = "header"
+    else:
+        start = end = source = None
+
+    return {
+        "version": QUALITY_VERSION,
+        "ground_density": round(n / (width * height), 6),
+        "density_grid": grid.tolist(),
+        "empty_fraction": round(empty_fraction, 4),
+        "acq_start": start,
+        "acq_end": end,
+        "acq_source": source,
+    }
+
+
+def quality_path(output_dir, basename):
+    """Chemin du sidecar qualité d'une dalle."""
+    return Path(output_dir) / QUALITY_DIRNAME / f"{basename}.json"
+
+
+def write_quality(output_dir, basename, data):
+    """Écrit le sidecar (atomique) ; False si le contenu est déjà identique."""
+    path = quality_path(output_dir, basename)
+    payload = json.dumps(data, ensure_ascii=False, sort_keys=True)
+    try:
+        if path.read_text(encoding="utf-8") == payload:
+            return False
+    except OSError:
+        pass
+    path.parent.mkdir(parents=True, exist_ok=True)
+    fd, tmp = tempfile.mkstemp(dir=str(path.parent), suffix=".tmp")
+    try:
+        with os.fdopen(fd, "w", encoding="utf-8") as f:
+            f.write(payload)
+        os.replace(tmp, path)
+    except OSError:
+        try:
+            os.unlink(tmp)
+        except OSError:
+            pass
+        raise
+    return True
+
+
+def read_quality(output_dir, basename):
+    """Sidecar d'une dalle, ou None (absent, illisible, autre version)."""
+    try:
+        data = json.loads(quality_path(output_dir, basename).read_text(encoding="utf-8"))
+    except (OSError, ValueError):
+        return None
+    if not isinstance(data, dict) or data.get("version") != QUALITY_VERSION:
+        return None
+    return data
+
+
+def load_quality_table(output_dir):
+    """Tous les sidecars valides : {basename: sidecar}."""
+    folder = Path(output_dir) / QUALITY_DIRNAME
+    table = {}
+    if not folder.is_dir():
+        return table
+    for f in sorted(folder.glob("*.json")):
+        data = read_quality(output_dir, f.stem)
+        if data is not None:
+            table[f.stem] = data
+    return table

lidar_pipeline/tests/test_quality.py

@@ -0,0 +1,91 @@
+"""Tests de la mesure de qualité des dalles (densité sol, date d'acquisition)."""
+
+BOUNDS = (652000.0, 6861000.0, 653000.0, 6862000.0)
+
+
+def _gps_seconds(year, month, day):
+    """Temps GPS ajusté standard (secondes GPS − 1e9) d'un midi UTC."""
+    from datetime import datetime, timezone
+    epoch = datetime(1980, 1, 6, tzinfo=timezone.utc)
+    return (datetime(year, month, day, 12, tzinfo=timezone.utc) - epoch).total_seconds() - 1e9
+
+
+def test_gps_adjusted_to_date():
+    from lidar_pipeline.quality import gps_adjusted_to_date
+    assert gps_adjusted_to_date(_gps_seconds(2023, 3, 15)) == "2023-03-15"
+
+
+def test_compute_quality_density_grid_and_dates():
+    import numpy as np
+    from lidar_pipeline.quality import compute_quality
+    # 1 point par m² sur la moitié nord, rien au sud ; 2 dates de vol.
+    xs, ys = np.meshgrid(np.arange(652000.5, 653000.0, 1.0),
+                         np.arange(6861500.5, 6862000.0, 1.0))
+    x, y = xs.ravel(), ys.ravel()
+    t = np.where(x < 652500, _gps_seconds(2023, 3, 15), _gps_seconds(2023, 3, 17))
+    q = compute_quality(x, y, t, BOUNDS)
+    assert q["version"] == 1
+    assert abs(q["ground_density"] - 0.5) < 1e-6
+    grid = q["density_grid"]
+    assert len(grid) == 20 and all(len(r) == 20 for r in grid)
+    assert grid[0][0] == 1.0          # ligne 0 = nord : couverte
+    assert grid[19][0] == 0.0         # sud : vide
+    assert q["acq_start"] == "2023-03-15" and q["acq_end"] == "2023-03-17"
+    assert q["acq_source"] == "gps"
+    # 1 point par m² : 1 pixel 0,2 m sur 25 occupé au nord, aucun au sud.
+    assert abs(q["empty_fraction"] - (1 - 0.5 / 25)) < 1e-3
+
+
+def test_compute_quality_ignores_points_outside_bounds():
+    import numpy as np
+    from lidar_pipeline.quality import compute_quality
+    x = np.array([652500.0, 660000.0]); y = np.array([6861500.0, 6861500.0])
+    q = compute_quality(x, y, None, BOUNDS, header_date="2024-05-03")
+    assert abs(q["ground_density"] - 1e-6) < 1e-9
+    assert q["acq_start"] == q["acq_end"] == "2024-05-03"
+    assert q["acq_source"] == "header"
+
+
+def test_compute_quality_week_time_falls_back_to_header():
+    import numpy as np
+    from lidar_pipeline.quality import compute_quality
+    x = np.array([652500.0]); y = np.array([6861500.0])
+    q = compute_quality(x, y, np.array([345600.0]), BOUNDS, gps_adjusted=False,
+                        header_date="2024-05-03")
+    assert q["acq_source"] == "header" and q["acq_start"] == "2024-05-03"
+
+
+def test_compute_quality_empty_tile():
+    import numpy as np
+    from lidar_pipeline.quality import compute_quality
+    q = compute_quality(np.array([]), np.array([]), np.array([]), BOUNDS)
+    assert q["ground_density"] == 0.0 and q["empty_fraction"] == 1.0
+    assert q["acq_start"] is None and q["acq_end"] is None and q["acq_source"] is None
+
+
+def test_sidecar_roundtrip_and_table(tmp_path):
+    from lidar_pipeline.quality import (load_quality_table, quality_path,
+                                        read_quality, write_quality)
+    base = "LHD_FXX_0652_6862_PTS_LAMB93_IGN69"
+    data = {"version": 1, "ground_density": 4.2}
+    assert write_quality(tmp_path, base, data) is True
+    assert write_quality(tmp_path, base, data) is False   # contenu identique : pas réécrit
+    assert quality_path(tmp_path, base).is_file()
+    assert read_quality(tmp_path, base) == data
+    assert load_quality_table(tmp_path) == {base: data}
+
+
+def test_read_quality_rejects_other_version(tmp_path):
+    from lidar_pipeline.quality import read_quality, write_quality
+    base = "LHD_FXX_0652_6862_PTS_LAMB93_IGN69"
+    write_quality(tmp_path, base, {"version": 0})
+    assert read_quality(tmp_path, base) is None
+
+
+def test_quality_module_imports_without_numpy():
+    """L'image légère n'a pas numpy : le module ne doit pas l'importer au chargement."""
+    import subprocess, sys
+    code = ("import sys; sys.modules['numpy'] = None; "
+            "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