diff --git a/lidar_pipeline/quality.py b/lidar_pipeline/quality.py new file mode 100644 index 0000000..1922ce6 --- /dev/null +++ b/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 diff --git a/lidar_pipeline/tests/test_quality.py b/lidar_pipeline/tests/test_quality.py new file mode 100644 index 0000000..0403778 --- /dev/null +++ b/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