Ajouter la géométrie de planche PDF, l'accès aux sources par emprise L93 et reportlab

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

Dockerfile.maps

@@ -11,12 +11,14 @@ ENV DEBIAN_FRONTEND=noninteractive
 ENV TZ=Europe/Paris
 
 # Pillow >= 11.3 : encodage/décodage AVIF natif (dalles et sous-tuiles du
-# pipeline) ; pyproj : reprojection Lambert 93 → Web Mercator sans GDAL.
+# pipeline) ; pyproj : reprojection Lambert 93 → Web Mercator sans GDAL ;
+# reportlab : planche PDF d'impression, pur Python.
 RUN pip3 install --no-cache-dir \
     fastapi \
     uvicorn \
     "pillow>=11.3" \
-    pyproj
+    pyproj \
+    "reportlab>=4.0"
 
 WORKDIR /app
 

lidar_pipeline/export_pdf.py

@@ -0,0 +1,119 @@
+"""Export PDF d'une zone : planche d'impression terrain du relief orienté.
+
+Composée directement en Lambert 93 depuis les sources de la pyramide
+(tiles.py) — échelle exacte, quadrillage aligné sur les dalles — puis
+dessinée avec reportlab (texte, grille et légende vectoriels). Tourne dans
+l'image légère : Pillow + pyproj + reportlab, sans numpy.
+"""
+
+import logging
+import math
+from dataclasses import dataclass
+
+logger = logging.getLogger("lidar")
+
+LAYER = "relief_oriente"
+PAPERS_MM = {"A4": (210.0, 297.0), "A3": (297.0, 420.0)}
+ORIENTS = ("portrait", "paysage")
+SCALES = (1000, 2000, 5000, 10000)
+DPI = {"A4": 300, "A3": 250}          # borne la mémoire du Pi (~36 Mo en A3)
+NATIVE_RES_M = 0.2
+MARGIN_MM = 10.0                      # bord non imprimable
+ANNOT_MM = 7.0                        # bande des coordonnées autour de la carte
+PANEL_SIDE_MM = 64.0                  # bandeau à droite (paysage)
+PANEL_BOTTOM_MM = 72.0                # bandeau en bas (portrait)
+_GRID_STEPS = {1000: 100, 2000: 100, 5000: 500, 10000: 1000}
+
+
+@dataclass(frozen=True)
+class Layout:
+    """Géométrie de la planche, en mm, origine en bas à gauche (reportlab)."""
+    paper: str
+    orient: str
+    scale: int
+    dpi: int
+    page_w: float
+    page_h: float
+    map_x: float
+    map_y: float
+    map_w: float
+    map_h: float
+    panel_x: float
+    panel_y: float
+    panel_w: float
+    panel_h: float
+
+
+def layout(paper, orient, scale):
+    """Géométrie d'une planche ; ValueError si un réglage est invalide."""
+    if paper not in PAPERS_MM:
+        raise ValueError(f"format inconnu : {paper} (A4 ou A3)")
+    if orient not in ORIENTS:
+        raise ValueError(f"orientation inconnue : {orient} (portrait ou paysage)")
+    try:
+        scale = int(scale)
+    except (TypeError, ValueError):
+        raise ValueError(f"échelle invalide : {scale}") from None
+    if scale not in SCALES:
+        raise ValueError("échelle non proposée : 1:" + str(scale)
+                         + " (1:1000, 1:2000, 1:5000 ou 1:10000)")
+    w, h = PAPERS_MM[paper]
+    if orient == "paysage":
+        w, h = h, w
+    inner = MARGIN_MM + ANNOT_MM
+    if orient == "paysage":
+        map_w = w - 2 * inner - PANEL_SIDE_MM
+        map_h = h - 2 * inner
+        return Layout(paper, orient, scale, DPI[paper], w, h, inner, inner, map_w, map_h,
+                      w - MARGIN_MM - PANEL_SIDE_MM, MARGIN_MM, PANEL_SIDE_MM, h - 2 * MARGIN_MM)
+    map_w = w - 2 * inner
+    map_h = h - 2 * inner - PANEL_BOTTOM_MM
+    return Layout(paper, orient, scale, DPI[paper], w, h, inner, MARGIN_MM + PANEL_BOTTOM_MM + ANNOT_MM,
+                  map_w, map_h, MARGIN_MM, MARGIN_MM, w - 2 * MARGIN_MM, PANEL_BOTTOM_MM)
+
+
+def center_l93(lat, lon):
+    """Centre WGS84 → Lambert 93 ; ValueError si non fini."""
+    from .tiles import wgs84_to_l93
+    lat, lon = float(lat), float(lon)
+    if not (math.isfinite(lat) and math.isfinite(lon)):
+        raise ValueError("coordonnées du centre invalides")
+    cx, cy = wgs84_to_l93(lon, lat)
+    if not (math.isfinite(cx) and math.isfinite(cy)):
+        raise ValueError("centre hors du domaine Lambert 93")
+    return cx, cy
+
+
+def map_bbox(cx, cy, lay):
+    """Emprise terrain L93 de la zone carte (papier × échelle)."""
+    half_w = lay.map_w / 1000.0 * lay.scale / 2
+    half_h = lay.map_h / 1000.0 * lay.scale / 2
+    return (cx - half_w, cy - half_h, cx + half_w, cy + half_h)
+
+
+def pixel_size(lay):
+    """Taille terrain d'un pixel imprimé (m), jamais plus fine que le natif."""
+    return max(NATIVE_RES_M, lay.scale * 0.0254 / lay.dpi)
+
+
+def grid_step(scale):
+    """Pas du quadrillage L93 (m) selon l'échelle."""
+    return _GRID_STEPS[int(scale)]
+
+
+def _to_wgs84(x, y):
+    from .tiles import _transformer
+    lon, lat = _transformer("EPSG:2154", "EPSG:4326").transform(x, y)
+    return lat, lon
+
+
+def frame(lat, lon, paper, orient, scale):
+    """Cadre imprimable pour la carte : emprise L93 et coins WGS84 (NO, NE, SE, SO)."""
+    lay = layout(paper, orient, scale)
+    cx, cy = center_l93(lat, lon)
+    b = map_bbox(cx, cy, lay)
+    corners = [_to_wgs84(b[0], b[3]), _to_wgs84(b[2], b[3]),
+               _to_wgs84(b[2], b[1]), _to_wgs84(b[0], b[1])]
+    return {"cx": cx, "cy": cy, "bbox_l93": list(b),
+            "corners": [[round(a, 7), round(o, 7)] for a, o in corners],
+            "width_m": round(b[2] - b[0]), "height_m": round(b[3] - b[1])}

lidar_pipeline/tests/test_export_pdf.py

@@ -0,0 +1,69 @@
+"""Tests de l'export PDF (planche d'impression terrain)."""
+
+import pytest
+
+
+def test_layout_landscape_a4_dimensions():
+    from lidar_pipeline.export_pdf import layout
+    lay = layout("A4", "paysage", 2000)
+    assert (lay.page_w, lay.page_h) == (297.0, 210.0) and lay.dpi == 300
+    assert lay.map_w > 150 and lay.map_h > 150
+    assert lay.panel_x >= lay.map_x + lay.map_w          # bandeau à droite
+    assert lay.map_x + lay.map_w <= lay.page_w and lay.map_y + lay.map_h <= lay.page_h
+
+
+def test_layout_portrait_a3_panel_below():
+    from lidar_pipeline.export_pdf import layout
+    lay = layout("A3", "portrait", 5000)
+    assert (lay.page_w, lay.page_h) == (297.0, 420.0) and lay.dpi == 250
+    assert lay.panel_y + lay.panel_h <= lay.map_y        # bandeau en bas
+
+
+@pytest.mark.parametrize("args", [("A5", "paysage", 2000), ("A4", "biais", 2000),
+                                  ("A4", "paysage", 1234)])
+def test_layout_rejects_invalid(args):
+    from lidar_pipeline.export_pdf import layout
+    with pytest.raises(ValueError):
+        layout(*args)
+
+
+def test_map_bbox_matches_paper_times_scale():
+    from lidar_pipeline.export_pdf import layout, map_bbox
+    lay = layout("A4", "paysage", 2000)
+    b = map_bbox(652500.0, 6861500.0, lay)
+    assert abs((b[2] - b[0]) - lay.map_w / 1000 * 2000) < 1e-6
+    assert abs((b[3] - b[1]) - lay.map_h / 1000 * 2000) < 1e-6
+    assert abs((b[0] + b[2]) / 2 - 652500.0) < 1e-6
+
+
+def test_pixel_size_and_grid_step():
+    from lidar_pipeline.export_pdf import grid_step, layout, pixel_size
+    assert pixel_size(layout("A4", "paysage", 1000)) == 0.2          # plafonné au natif
+    assert abs(pixel_size(layout("A3", "paysage", 10000)) - 10000 * 0.0254 / 250) < 1e-9
+    assert [grid_step(s) for s in (1000, 2000, 5000, 10000)] == [100, 100, 500, 1000]
+
+
+def test_center_l93_rejects_non_finite():
+    from lidar_pipeline.export_pdf import center_l93
+    with pytest.raises(ValueError):
+        center_l93(float("nan"), 2.0)
+
+
+def test_frame_roundtrip():
+    from lidar_pipeline.export_pdf import center_l93, frame
+    lat, lon = 48.85, 2.35
+    f = frame(lat, lon, "A4", "paysage", 2000)
+    cx, cy = center_l93(lat, lon)
+    assert abs(f["cx"] - cx) < 1e-6 and abs(f["cy"] - cy) < 1e-6
+    assert len(f["corners"]) == 4
+    nw, ne, se, sw = f["corners"]
+    assert nw[0] > sw[0] and ne[1] > nw[1]               # [lat, lon] : nord en haut, est à droite
+    assert f["width_m"] > f["height_m"] > 0
+
+
+def test_export_modules_import_without_numpy():
+    import subprocess, sys
+    code = ("import sys; sys.modules['numpy'] = None; "
+            "import lidar_pipeline.export_pdf, lidar_pipeline.tiles; print('ok')")
+    r = subprocess.run([sys.executable, "-c", code], capture_output=True, text=True)
+    assert r.returncode == 0, r.stderr

lidar_pipeline/tests/test_tiles.py

@@ -166,6 +166,20 @@ def test_pick_tier_by_resolution(tmp_path):
     assert tiles._pick_tier(tiers, 0.01)[0].res == tiers[-1][0].res
 
 
+def test_sources_in_bbox_picks_tier(tmp_path):
+    from lidar_pipeline import tiles
+    _make_dalle(tmp_path, 1054, 6882, ["aspect"], px=64)
+    tiles.source_index(tmp_path, force=True)
+    bbox = (1054000.0, 6881000.0, 1055000.0, 6882000.0)
+    fine = tiles.sources_in_bbox(tmp_path, "aspect", bbox, 0.2)
+    coarse = tiles.sources_in_bbox(tmp_path, "aspect", bbox, 50.0)
+    assert fine and coarse and fine[0][0] == (1054, 6882)
+    assert min(s.res for _c, s in fine) <= min(s.res for _c, s in coarse)
+    assert tiles.sources_in_bbox(tmp_path, "aspect", (0, 0, 10, 10), 0.2) == []
+    img = tiles.load_source(fine[0][1])
+    assert img is not None and img.mode in ("RGB", "RGBA")
+
+
 def test_subtiles_preferred_over_full_dalle(tmp_path):
     """Les quadrants index_subtiles servent de palier fin (4× moins à décoder)."""
     import pytest

lidar_pipeline/tiles.py

@@ -753,39 +753,52 @@ def _pick_tier(tiers, target_res):
     return tiers[-1]
 
 
-def _contributing(output_dir, layer, z, x, y, scale):
-    """Sources intersectant la tuile, palier choisi selon la résolution visée."""
+def sources_in_bbox(output_dir, layer, bbox, target_res):
+    """Sources d'une couche intersectant une emprise L93, au palier le plus
+    grossier suffisant pour target_res (m/px). Retourne [((col, row), source)]."""
     per_cell = source_index(output_dir).get(layer)
     if not per_cell:
         return []
-    min_x, min_y, max_x, max_y = tile_bounds_l93(z, x, y)
-    target = target_resolution(z, y, scale)
+    min_x, min_y, max_x, max_y = bbox
     out = []
     for (col, row), tiers in per_cell.items():
         b = _cell_bounds(col, row)
         if b[2] <= min_x or b[0] >= max_x or b[3] <= min_y or b[1] >= max_y:
             continue
-        for src in _pick_tier(tiers, target):
+        for src in _pick_tier(tiers, target_res):
             s = src.bounds
             if s[2] <= min_x or s[0] >= max_x or s[3] <= min_y or s[1] >= max_y:
                 continue
-            out.append(src)
+            out.append(((col, row), src))
     return out
 
 
+def _contributing(output_dir, layer, z, x, y, scale):
+    """Sources intersectant la tuile, palier choisi selon la résolution visée."""
+    return [src for _cell, src in sources_in_bbox(
+        output_dir, layer, tile_bounds_l93(z, x, y), target_resolution(z, y, scale))]
+
+
+def load_source(src):
+    """Image décodée d'une source (rapatriée si besoin), ou None si illisible."""
+    if not src.ensure():
+        return None
+    mtime = src.mtime()
+    if mtime is None:
+        return None
+    try:
+        return _open_source(str(src.path), mtime)
+    except Exception as e:  # noqa: BLE001 — source illisible : rendu partiel
+        logger.debug(f"Source illisible ({src.path.name}) : {e}")
+        return None
+
+
 def _paste_source(canvas, src, z, x, y, size, resample):
     """Reprojette une source dans la tuile (transformation projective)."""
     from PIL import Image
 
-    if not src.ensure():
-        return False
-    mtime = src.mtime()
-    if mtime is None:
-        return False
-    try:
-        img = _open_source(str(src.path), mtime)
-    except Exception as e:  # noqa: BLE001 — source illisible : tuile partielle
-        logger.debug(f"Source de tuile illisible ({src.path.name}) : {e}")
+    img = load_source(src)
+    if img is None:
         return False
 
     w, h = img.size

requirements.txt

@@ -16,6 +16,7 @@ fastapi>=0.110
 uvicorn>=0.29
 pydantic>=2.6
 pyproj>=3.6
+reportlab>=4.0  # Export PDF de la carte (mapserve)
 
 # Tests (./run.sh --test)
 pytest>=7.4