Agréger la qualité de zone et ajouter les couleurs de légende de la planche PDF

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

lidar_pipeline/export_pdf.py

@@ -192,3 +192,116 @@ def compose_l93(output_dir, bbox, px_size, layer=LAYER):
     if mask.getextrema() != (255, 255):
         canvas.paste(_hatch(canvas.size), (0, 0), ImageOps.invert(mask))
     return canvas, mask, sorted(cells)
+
+
+ROSE_L = 64.0
+ROSE_CHROMA = 60.0           # = visualizations.RELIEF_CHROMA
+# Classes de densité de points sol (pts/m²) : rouge = donnée faible.
+DENSITY_CLASSES = [
+    (0.0, "moins de 1", "#d7301f"),
+    (1.0, "1 à 3", "#fc8d59"),
+    (3.0, "3 à 6", "#fdcc8a"),
+    (6.0, "6 à 10", "#a1d99b"),
+    (10.0, "10 et plus", "#31a354"),
+]
+
+
+def lab_to_rgb(L, a, b):
+    """CIELAB (D65) → sRGB 8 bits, même formule que la carte (labToRgb)."""
+    fy = (L + 16) / 116
+    fx, fz = fy + a / 500, fy - b / 200
+
+    def finv(t):
+        return t ** 3 if t > 6 / 29 else 3 * (6 / 29) ** 2 * (t - 4 / 29)
+
+    X, Y, Z = 0.95047 * finv(fx), finv(fy), 1.08883 * finv(fz)
+    lin = (3.2406 * X - 1.5372 * Y - 0.4986 * Z,
+           -0.9689 * X + 1.8758 * Y + 0.0415 * Z,
+           0.0557 * X - 0.2040 * Y + 1.0570 * Z)
+    out = []
+    for c in lin:
+        v = 12.92 * c if c <= 0.0031308 else 1.055 * max(c, 0.0) ** (1 / 2.4) - 0.055
+        out.append(int(round(min(1.0, max(0.0, v)) * 255)))
+    return tuple(out)
+
+
+def rose_color(compass_deg):
+    """Couleur d'une orientation de pente (0 = N, sens horaire), comme la rose de la carte."""
+    chroma = ROSE_CHROMA * min(1.0, ROSE_L * (100 - ROSE_L) / 2500)
+    h = math.radians((compass_deg + 90) % 360)
+    return lab_to_rgb(ROSE_L, chroma * math.cos(h), chroma * math.sin(h))
+
+
+def density_color(v):
+    """Couleur de classe d'une densité sol (pts/m²)."""
+    color = DENSITY_CLASSES[0][2]
+    for low, _label, c in DENSITY_CLASSES:
+        if v >= low:
+            color = c
+    return color
+
+
+def _pdf_text(s):
+    """Texte encodable par les polices standard (WinAnsi/cp1252)."""
+    return "".join(ch if ch.encode("cp1252", "ignore") else "?" for ch in str(s))
+
+
+def _cells_of_bbox(bbox):
+    """Dalles LHD 1 km intersectant une emprise L93."""
+    c0, c1 = int(math.floor(bbox[0] / 1000)), int(math.ceil(bbox[2] / 1000))
+    r0, r1 = int(math.floor(bbox[1] / 1000)) + 1, int(math.ceil(bbox[3] / 1000))
+    return [(c, r) for c in range(c0, c1) for r in range(r1, r0 - 1, -1)]
+
+
+def zone_quality(bbox, table, cells_with_relief):
+    """Agrège la qualité des dalles sur une emprise (pondérée par la surface)."""
+    from .index import parse_basename_coords
+    from .quality import DENSITY_CELL_M
+
+    by_cell = {}
+    for base, data in table.items():
+        coords = parse_basename_coords(base)
+        if coords is not None:
+            by_cell[tuple(coords)] = data
+    cells = _cells_of_bbox(bbox)
+    relief = set(map(tuple, cells_with_relief))
+    total_w = dens_w = empty_w = 0.0
+    dmin = None
+    starts, ends, sources = [], [], set()
+    grid_cells, missing_q = [], []
+    for col, row in cells:
+        q = by_cell.get((col, row))
+        if q is None:
+            missing_q.append((col, row))
+            continue
+        x0, y1 = col * 1000.0, row * 1000.0
+        grid = q.get("density_grid") or []
+        for j, line in enumerate(grid):
+            for i, v in enumerate(line):
+                gx0, gx1 = x0 + i * DENSITY_CELL_M, x0 + (i + 1) * DENSITY_CELL_M
+                gy1, gy0 = y1 - j * DENSITY_CELL_M, y1 - (j + 1) * DENSITY_CELL_M
+                ox = min(gx1, bbox[2]) - max(gx0, bbox[0])
+                oy = min(gy1, bbox[3]) - max(gy0, bbox[1])
+                if ox <= 0 or oy <= 0:
+                    continue
+                area = ox * oy
+                total_w += area
+                dens_w += v * area
+                empty_w += float(q.get("empty_fraction") or 0.0) * area
+                dmin = v if dmin is None else min(dmin, v)
+                grid_cells.append((gx0, gy0, gx1, gy1, v))
+        if q.get("acq_start"):
+            starts.append(q["acq_start"]); ends.append(q["acq_end"] or q["acq_start"])
+            sources.add(q.get("acq_source"))
+    return {
+        "density_mean": dens_w / total_w if total_w else None,
+        "density_min": dmin,
+        "empty_fraction": empty_w / total_w if total_w else None,
+        "acq_start": min(starts) if starts else None,
+        "acq_end": max(ends) if ends else None,
+        "acq_sources": sources,
+        "cells": cells,
+        "missing_relief": [c for c in cells if c not in relief],
+        "missing_quality": missing_q,
+        "grid_cells": grid_cells,
+    }

lidar_pipeline/tests/test_export_pdf.py

@@ -119,3 +119,59 @@ def test_compose_l93_no_data_returns_empty_cells(tmp_path):
     from lidar_pipeline.export_pdf import compose_l93
     img, mask, cells = compose_l93(tmp_path, (0.0, 0.0, 100.0, 100.0), 1.0)
     assert cells == [] and mask.getextrema() == (0, 0)
+
+
+def test_lab_to_rgb_matches_numpy_reference():
+    import numpy as np
+    from lidar_pipeline.export_pdf import lab_to_rgb
+    from lidar_pipeline.visualizations import _lab_to_srgb
+    for L, a, b in [(64, 30, -20), (20, 0, 0), (90, -10, 40), (50, 60, 0)]:
+        ref = tuple(int(round(v * 255)) for v in _lab_to_srgb(
+            np.array(L, float), np.array(a, float), np.array(b, float)))
+        assert all(abs(p - q) <= 1 for p, q in zip(lab_to_rgb(L, a, b), ref))
+
+
+def test_rose_color_distinct_orientations():
+    from lidar_pipeline.export_pdf import rose_color
+    colors = {rose_color(c) for c in (0, 90, 180, 270)}
+    assert len(colors) == 4
+
+
+def test_density_color_classes():
+    from lidar_pipeline.export_pdf import DENSITY_CLASSES, density_color
+    assert density_color(0.0) == DENSITY_CLASSES[0][2]
+    assert density_color(50.0) == DENSITY_CLASSES[-1][2]
+
+
+def test_pdf_text_replaces_unencodable():
+    from lidar_pipeline.export_pdf import _pdf_text
+    assert _pdf_text("Relief orienté — 1:2 000 ≥ 🚀") == "Relief orienté — 1:2 000 ? ?"
+
+
+def _q(density, start="2023-03-15", end="2023-03-17", empty=0.1):
+    return {"version": 1, "ground_density": density,
+            "density_grid": [[density] * 20 for _ in range(20)],
+            "empty_fraction": empty, "acq_start": start, "acq_end": end,
+            "acq_source": "gps"}
+
+
+def test_zone_quality_aggregates_two_cells():
+    from lidar_pipeline.export_pdf import zone_quality
+    table = {"LHD_FXX_1054_6882_PTS_LAMB93_IGN69": _q(4.0, "2023-03-15", "2023-03-15"),
+             "LHD_FXX_1055_6882_PTS_LAMB93_IGN69": _q(8.0, "2023-04-02", "2023-04-03", 0.3)}
+    bbox = (1054500.0, 6881500.0, 1055500.0, 6881600.0)     # moitié de chaque dalle
+    z = zone_quality(bbox, table, [(1054, 6882), (1055, 6882)])
+    assert abs(z["density_mean"] - 6.0) < 1e-6 and z["density_min"] == 4.0
+    assert abs(z["empty_fraction"] - 0.2) < 1e-6
+    assert (z["acq_start"], z["acq_end"]) == ("2023-03-15", "2023-04-03")
+    assert z["missing_relief"] == [] and z["missing_quality"] == []
+    assert len(z["grid_cells"]) == 2 * 10 * 2        # 10 mailles en x × 2 en y par dalle
+
+
+def test_zone_quality_missing_cells():
+    from lidar_pipeline.export_pdf import zone_quality
+    bbox = (1054500.0, 6881500.0, 1055500.0, 6881600.0)
+    z = zone_quality(bbox, {}, [(1054, 6882)])
+    assert z["density_mean"] is None and z["acq_start"] is None
+    assert z["missing_quality"] == [(1054, 6882), (1055, 6882)]
+    assert z["missing_relief"] == [(1055, 6882)]