053eb2f529
L'ancienne carte n'est pas une carte à tuiles : une div Leaflet par dalle,
rotée en CSS pour coller la grille Lambert 93 sur le Web Mercator, trois
paliers d'images choisis à la main, un plafond d'images pleine résolution et
une mosaïque d'overview pour boucher les trous au dézoom. Des centaines de
nœuds DOM, des AVIF de 2500² à 5000² décodés dans le navigateur, un LOD
maison — et rien de réutilisable hors de cette page.
Nouvelle image légère (Pillow + pyproj, ni GDAL ni PDAL, port 8975) servant
une pyramide XYZ EPSG:3857 au schéma OpenStreetMap, rendue à la demande
depuis les dalles et mise en cache dans output/index_xyz/ :
- tiles.py : grille XYZ, reprojection par transformation projective dalle par
dalle (calage mesuré < 1 px), choix du palier source parmi ceux que le
pipeline produit déjà (vignette, intermédiaire, quadrant, dalle), cache
périmé dès qu'une dalle contributrice est plus récente, marqueur .empty
pour les zones sans donnée, cache d'images sources à budget mémoire ;
- mapserve.py : /tiles/{couche}/{z}/{x}/{y}.png (256 px canonique) et
@2x.webp (512 px, interface), TileJSON, WMTS, josm.imagery.xml, CORS —
les rendus deviennent un fond d'imagerie pour JOSM, iD, QGIS, uMap ;
- mapui.py : une L.tileLayer par couche dans un pane isolé — LOD, cache et
animation natifs de Leaflet ; pile réordonnable (glisser avec barre
d'insertion, ou boutons ▲▼ au doigt), modes de fusion CSS par couche,
configuration figeable comme défaut de tous les navigateurs.
Deux amonts pour les déploiements en deux machines : LIDAR_SOURCE_URL
(webapp du pipeline — inventaire complet, dalles rapatriées à la demande) et
LIDAR_MAPS_URL (autre instance carte). Charge bornée et réglable, taillée
par défaut pour une petite machine : 2 rendus et 2 téléchargements
simultanés, 192 Mo de cache de sources.
Mesuré sur 849 dalles réelles : première tuile d'une zone 0,9–3,2 s
(téléchargement compris), ~1 ms ensuite ; un chemin OSM se superpose
exactement à la trace du rendu de pente, et les coutures entre tuiles
restent sous le bruit naturel du terrain.
L'interface historique (port 8973) n'est pas touchée : génération et export
y restent, les deux cartes coexistent.
834 lines
32 KiB
Python
834 lines
32 KiB
Python
"""Pyramide de tuiles XYZ (EPSG:3857) rendue à la demande depuis les dalles.
|
||
|
||
Schéma de tuilage identique à celui d'OpenStreetMap / Google Maps : grille
|
||
Web Mercator, origine au coin nord-ouest, `{z}/{x}/{y}`, tuiles de 256 px
|
||
(512 px avec `scale=2`, convention `@2x`). Les rendus du pipeline restent des
|
||
dalles Lambert 93 de 1 km : chaque tuile est composée à la volée en reprojetant
|
||
les dalles qui l'intersectent (transformation projective par dalle, erreur très
|
||
inférieure au pixel), puis mise en cache sur disque.
|
||
|
||
Volontairement sans GDAL ni numpy : Pillow + pyproj suffisent, l'image légère
|
||
(Dockerfile.maps / Dockerfile.webapp) reste petite et portable ARM64.
|
||
"""
|
||
|
||
import logging
|
||
import math
|
||
import os
|
||
import re
|
||
import threading
|
||
import time
|
||
from collections import OrderedDict
|
||
from functools import lru_cache
|
||
from pathlib import Path
|
||
|
||
logger = logging.getLogger("lidar")
|
||
|
||
# --- Contrat de tuilage (cf. docs/MAPS.md) ---------------------------------
|
||
TILE_SIZE = 256 # taille canonique (OSM/XYZ) ; @2x → 512
|
||
TILE_MIN_Z = 5
|
||
TILE_MAX_NATIVE_Z = 19 # 0,2 m/px ≈ résolution du z19 à la latitude 47°
|
||
TILE_DIRNAME = "index_xyz" # cache disque, sous le dossier de sortie
|
||
WEBP_QUALITY = 78
|
||
AVIF_QUALITY = 60
|
||
# PNG palettisé (PNG8 + alpha) : ~5× plus léger (170 → 32 Ko sur une dalle
|
||
# réelle) pour un écart moyen de ~4 niveaux sur une rampe de couleur. Laissé
|
||
# DÉSACTIVÉ par défaut : le PNG canonique reste sans perte, la fidélité prime
|
||
# sur le débit pour un produit d'interprétation. `LIDAR_TILE_PNG_PALETTE=1`
|
||
# l'active quand la bande passante compte (consultation mobile).
|
||
PNG_PALETTE = os.environ.get("LIDAR_TILE_PNG_PALETTE", "") == "1"
|
||
|
||
# Demi-circonférence équatoriale : emprise du Web Mercator (EPSG:3857).
|
||
ORIGIN = 20037508.342789244
|
||
|
||
# Résolutions nominales des paliers de source réutilisés tels quels par le
|
||
# pipeline (m/px) : vignette 256 px/km, vignette intermédiaire 640 px/km.
|
||
_THUMB_RES = 1000.0 / 256
|
||
_MID_RES = 1000.0 / 640
|
||
|
||
_SUBTILE_RE = re.compile(r"_(\d+)_(\d+)\.avif$")
|
||
|
||
# Serveur de dalles amont (webapp du pipeline) : quand il est défini, l'index
|
||
# des sources vient de son /api/tiles et les images manquantes sont rapatriées
|
||
# à la demande dans le cache local — le conteneur carte n'a alors besoin
|
||
# d'aucune donnée locale au démarrage.
|
||
REMOTE_SOURCE_URL = (os.environ.get("LIDAR_SOURCE_URL") or "").rstrip("/")
|
||
REMOTE_SOURCE_TOKEN = os.environ.get("LIDAR_SOURCE_TOKEN") or None
|
||
_REMOTE_TTL = 60.0
|
||
_remote_cache = {"payload": None, "at": 0.0, "index": None, "root": None}
|
||
# Téléchargements de dalles simultanés (réseau) — modeste par défaut : sur un
|
||
# petit serveur, chaque source rapatriée est un fichier de plusieurs Mo.
|
||
FETCH_WORKERS = max(1, int(os.environ.get("LIDAR_TILE_FETCH_WORKERS", "2") or 2))
|
||
_fetch_sem = threading.Semaphore(FETCH_WORKERS)
|
||
_fetch_locks = {}
|
||
_fetch_guard = threading.Lock()
|
||
|
||
# Index des sources reconstruit au plus toutes les _INDEX_TTL secondes (ou dès
|
||
# qu'un dossier change de mtime : ajout/suppression de dalle).
|
||
_INDEX_TTL = 20.0
|
||
_index_cache = {}
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# Géométrie de la grille
|
||
# ---------------------------------------------------------------------------
|
||
|
||
def tile_bounds_3857(z, x, y):
|
||
"""Emprise (ouest, sud, est, nord) d'une tuile XYZ en mètres EPSG:3857."""
|
||
span = 2.0 * ORIGIN / (2 ** z)
|
||
west = -ORIGIN + x * span
|
||
north = ORIGIN - y * span
|
||
return west, north - span, west + span, north
|
||
|
||
|
||
def ground_resolution(z, scale=1):
|
||
"""Résolution d'une tuile à l'équateur (m/px) ; ×cos(lat) sur le terrain."""
|
||
return 2.0 * ORIGIN / (TILE_SIZE * scale * (2 ** z))
|
||
|
||
|
||
def tile_latitude(z, y):
|
||
"""Latitude (degrés) du centre d'une tuile — sert au calage de résolution."""
|
||
n = math.pi - 2.0 * math.pi * (y + 0.5) / (2 ** z)
|
||
return math.degrees(math.atan(math.sinh(n)))
|
||
|
||
|
||
def target_resolution(z, y, scale=1):
|
||
"""Résolution terrain visée par la tuile (m/px), latitude comprise."""
|
||
return ground_resolution(z, scale) * math.cos(math.radians(tile_latitude(z, y)))
|
||
|
||
|
||
@lru_cache(maxsize=4)
|
||
def _transformer(src, dst):
|
||
from pyproj import Transformer
|
||
return Transformer.from_crs(src, dst, always_xy=True)
|
||
|
||
|
||
def to_l93(xs, ys):
|
||
"""EPSG:3857 → EPSG:2154 (listes de coordonnées)."""
|
||
return _transformer("EPSG:3857", "EPSG:2154").transform(xs, ys)
|
||
|
||
|
||
def to_3857(xs, ys):
|
||
"""EPSG:2154 → EPSG:3857 (listes de coordonnées)."""
|
||
return _transformer("EPSG:2154", "EPSG:3857").transform(xs, ys)
|
||
|
||
|
||
@lru_cache(maxsize=4096)
|
||
def wgs84_to_l93(lon, lat):
|
||
"""Point WGS84 → Lambert 93 (utilisé par la fiche d'information dalle)."""
|
||
return _transformer("EPSG:4326", "EPSG:2154").transform(lon, lat)
|
||
|
||
|
||
@lru_cache(maxsize=4096)
|
||
def tile_bounds_l93(z, x, y, samples=5):
|
||
"""Emprise L93 englobant une tuile XYZ.
|
||
|
||
Les bords d'une tuile ne sont pas des droites en Lambert 93 : on
|
||
échantillonne une grille `samples`×`samples` plutôt que les seuls coins,
|
||
sinon l'emprise est sous-estimée aux petits zooms (tuiles de centaines de km).
|
||
"""
|
||
west, south, east, north = tile_bounds_3857(z, x, y)
|
||
xs, ys = [], []
|
||
for i in range(samples):
|
||
for j in range(samples):
|
||
xs.append(west + (east - west) * i / (samples - 1))
|
||
ys.append(south + (north - south) * j / (samples - 1))
|
||
lx, ly = to_l93(xs, ys)
|
||
return min(lx), min(ly), max(lx), max(ly)
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# Transformation projective (mapping sortie → source, convention Pillow)
|
||
# ---------------------------------------------------------------------------
|
||
|
||
def _solve(matrix, rhs):
|
||
"""Résout un système linéaire dense (pivot partiel), sans numpy."""
|
||
n = len(rhs)
|
||
a = [row[:] + [rhs[i]] for i, row in enumerate(matrix)]
|
||
for col in range(n):
|
||
piv = max(range(col, n), key=lambda r: abs(a[r][col]))
|
||
if abs(a[piv][col]) < 1e-12:
|
||
return None
|
||
a[col], a[piv] = a[piv], a[col]
|
||
inv = 1.0 / a[col][col]
|
||
for k in range(col, n + 1):
|
||
a[col][k] *= inv
|
||
for r in range(n):
|
||
if r == col:
|
||
continue
|
||
f = a[r][col]
|
||
if f:
|
||
for k in range(col, n + 1):
|
||
a[r][k] -= f * a[col][k]
|
||
return [a[i][n] for i in range(n)]
|
||
|
||
|
||
def perspective_coeffs(dst_quad, src_quad):
|
||
"""Coefficients Pillow `Image.PERSPECTIVE` mappant sortie → source.
|
||
|
||
Pillow échantillonne la SOURCE en (x', y') = ((a x + b y + c) / (g x + h y + 1),
|
||
(d x + e y + f) / (g x + h y + 1)) pour chaque pixel (x, y) de la SORTIE :
|
||
on résout donc les 8 inconnues à partir de 4 correspondances
|
||
(point de sortie → point source). Retourne None si le quadrilatère est
|
||
dégénéré (dalle réduite à un point au dézoom extrême).
|
||
"""
|
||
matrix, rhs = [], []
|
||
for (dx, dy), (sx, sy) in zip(dst_quad, src_quad):
|
||
matrix.append([dx, dy, 1, 0, 0, 0, -sx * dx, -sx * dy])
|
||
rhs.append(sx)
|
||
matrix.append([0, 0, 0, dx, dy, 1, -sy * dx, -sy * dy])
|
||
rhs.append(sy)
|
||
return _solve(matrix, rhs)
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# Index des sources : dalles du pipeline, par couche et par palier
|
||
# ---------------------------------------------------------------------------
|
||
|
||
class _Source:
|
||
"""Une image source géoréférencée (dalle entière ou quadrant).
|
||
|
||
`url` non nul : l'image vit sur le serveur de dalles amont et n'est
|
||
rapatriée qu'au premier besoin (`ensure()`), dans `path`. `version` est
|
||
alors la mtime (ms) annoncée par l'amont : elle sert de date de référence
|
||
pour la péremption des tuiles, avant même tout téléchargement.
|
||
"""
|
||
|
||
__slots__ = ("path", "bounds", "res", "url", "version")
|
||
|
||
def __init__(self, path, bounds, res, url=None, version=None):
|
||
self.path = path
|
||
self.bounds = bounds # (min_x, min_y, max_x, max_y) en L93
|
||
self.res = res # résolution nominale (m/px)
|
||
self.url = url
|
||
self.version = version
|
||
|
||
def mtime(self):
|
||
if self.version is not None:
|
||
return self.version / 1000.0
|
||
try:
|
||
return self.path.stat().st_mtime
|
||
except OSError:
|
||
return None
|
||
|
||
def ensure(self):
|
||
"""Garantit la présence locale de l'image (rapatriement si besoin)."""
|
||
if self.url is None or self.path.is_file():
|
||
return self.path.is_file()
|
||
return _fetch_source(self.url, self.path)
|
||
|
||
|
||
def _remote_thumb_px(k):
|
||
"""Côté (px) de la vignette servie par l'amont : 256 par dalle, 160 par quadrant."""
|
||
return 160.0 if k > 1 else 256.0
|
||
|
||
|
||
def _fetch_source(url, dest):
|
||
"""Télécharge une image source depuis l'amont (écriture atomique)."""
|
||
import urllib.request
|
||
with _fetch_guard:
|
||
lock = _fetch_locks.setdefault(str(dest), threading.Lock())
|
||
with lock:
|
||
if dest.is_file():
|
||
return True
|
||
try:
|
||
req = urllib.request.Request(
|
||
url, headers={"User-Agent": "lidar-maps-source"})
|
||
if REMOTE_SOURCE_TOKEN:
|
||
req.add_header("X-Lidar-Token", REMOTE_SOURCE_TOKEN)
|
||
with _fetch_sem, urllib.request.urlopen(req, timeout=120) as r:
|
||
data = r.read()
|
||
except Exception as e: # noqa: BLE001 — amont éteint : tuile partielle
|
||
logger.debug(f"Source amont indisponible ({url}) : {e}")
|
||
return False
|
||
_write_atomic(dest, data)
|
||
logger.info(f"Source rapatriée : {dest.name} ({len(data) / 1e6:.1f} Mo)")
|
||
return dest.is_file()
|
||
|
||
|
||
def _remote_payload(force=False):
|
||
"""Index du serveur de dalles amont (/api/tiles), en cache 60 s."""
|
||
import json
|
||
import urllib.request
|
||
now = time.time()
|
||
if not force and _remote_cache["payload"] is not None \
|
||
and now - _remote_cache["at"] < _REMOTE_TTL:
|
||
return _remote_cache["payload"]
|
||
try:
|
||
req = urllib.request.Request(f"{REMOTE_SOURCE_URL}/api/tiles",
|
||
headers={"User-Agent": "lidar-maps-source"})
|
||
with urllib.request.urlopen(req, timeout=30) as r:
|
||
payload = json.loads(r.read().decode("utf-8"))
|
||
except Exception as e: # noqa: BLE001 — on garde le dernier index connu
|
||
logger.warning(f"Index amont injoignable ({REMOTE_SOURCE_URL}) : {e}")
|
||
payload = _remote_cache["payload"]
|
||
_remote_cache.update(payload=payload, at=now)
|
||
return payload
|
||
|
||
|
||
def _remote_index(output_dir, force=False):
|
||
"""Inventaire construit depuis l'amont : sources non encore rapatriées.
|
||
|
||
Le résultat est mémoïsé avec la charge utile : le reconstruire coûte des
|
||
centaines de millisecondes sur un catalogue de plusieurs milliers de
|
||
dalles, ce qui se paierait à CHAQUE tuile servie.
|
||
"""
|
||
payload = _remote_payload(force)
|
||
if not payload or not payload.get("tiles"):
|
||
return {}
|
||
if (_remote_cache["index"] is not None
|
||
and _remote_cache["root"] == str(output_dir)
|
||
and _remote_cache.get("built") is payload):
|
||
return _remote_cache["index"]
|
||
output_dir = Path(output_dir)
|
||
layers = {}
|
||
for entry in payload["tiles"]:
|
||
col, row = entry.get("col"), entry.get("row")
|
||
if col is None or row is None:
|
||
continue
|
||
res = float(entry.get("resolution") or 0.5)
|
||
k = int(entry.get("sub_k") or 1)
|
||
step = 1000.0 / k
|
||
i, j = int(entry.get("sub_i") or 0), int(entry.get("sub_j") or 0)
|
||
base = _cell_bounds(col, row)
|
||
bounds = (base[0] + i * step, base[1] + j * step,
|
||
base[0] + (i + 1) * step, base[1] + (j + 1) * step)
|
||
for viz, info in (entry.get("viz") or {}).items():
|
||
by_res = layers.setdefault(viz, {}).setdefault((col, row), {})
|
||
for tier, nominal in (("thumb", step / _remote_thumb_px(k)),
|
||
("mid", step / 640.0),
|
||
("full", res)):
|
||
url_rel = info.get(tier)
|
||
if not url_rel:
|
||
continue
|
||
rel, _, query = url_rel.partition("?")
|
||
version = None
|
||
if query.startswith("v="):
|
||
try:
|
||
version = int(query[2:])
|
||
except ValueError:
|
||
version = None
|
||
rank = 0 if k > 1 else 1
|
||
by_res.setdefault((nominal, rank), []).append(_Source(
|
||
output_dir / rel, bounds, nominal,
|
||
url=f"{REMOTE_SOURCE_URL}/{url_rel}", version=version))
|
||
out = {}
|
||
for viz, per_cell in layers.items():
|
||
for cell, by_res in per_cell.items():
|
||
out.setdefault(viz, {})[cell] = [
|
||
by_res[key] for key in sorted(by_res, key=lambda t: (-t[0], t[1]))]
|
||
_remote_cache.update(index=out, root=str(output_dir), built=payload)
|
||
return out
|
||
|
||
|
||
def _cell_bounds(col, row):
|
||
"""Emprise L93 d'une dalle LHD : X ∈ [col, col+1] km, Y ∈ [row-1, row] km."""
|
||
return (col * 1000.0, (row - 1) * 1000.0, (col + 1) * 1000.0, row * 1000.0)
|
||
|
||
|
||
def _dir_mtime(path):
|
||
try:
|
||
return path.stat().st_mtime
|
||
except OSError:
|
||
return None
|
||
|
||
|
||
def _split_viz_key(stem, known):
|
||
"""`{dir_name}_{viz}` → (dir_name, viz) en s'appuyant sur les clés connues.
|
||
|
||
Les clés de visualisation contiennent des soulignés (`positive_openness`,
|
||
`hillshade_multi`) : on ne peut pas couper au dernier `_`, on reconnaît un
|
||
suffixe connu (le plus long d'abord).
|
||
"""
|
||
for key in known:
|
||
if stem.endswith("_" + key):
|
||
return stem[:-len(key) - 1], key
|
||
return None, None
|
||
|
||
|
||
def _cell_of_dir(dir_name):
|
||
"""(col, row, résolution) d'un nom de dossier de dalle, ou None."""
|
||
from .index import _strip_res_suffix, parse_basename_coords
|
||
coords = parse_basename_coords(dir_name)
|
||
if coords is None:
|
||
return None
|
||
_base, res = _strip_res_suffix(dir_name)
|
||
return coords[0], coords[1], res
|
||
|
||
|
||
def _build_index(output_dir):
|
||
"""Inventaire {couche: {(col, row): [paliers du plus grossier au plus fin]}}.
|
||
|
||
Les trois sources du pipeline sont scannées INDÉPENDAMMENT — dalles
|
||
(`visualisations/`), quadrants (`index_subtiles/`) et vignettes
|
||
(`index_thumbs/`). Un cache partiel reste donc exploitable : sur une
|
||
machine légère, seuls quadrants et vignettes sont rapatriés, jamais les
|
||
dalles entières.
|
||
"""
|
||
# _SUBTILE_THUMB_PX : taille des vignettes de quadrant, définie par
|
||
# l'index (la changer là-bas doit rester sans effet ici).
|
||
from .index import _SUBTILE_THUMB_PX, VIZ_LABELS, scan_tiles
|
||
|
||
output_dir = Path(output_dir)
|
||
vis_dir = output_dir / "visualisations"
|
||
thumb_dir = output_dir / "index_thumbs"
|
||
sub_dir = output_dir / "index_subtiles"
|
||
|
||
# res → sources, par couche et par dalle ; trié en paliers à la fin.
|
||
records = {}
|
||
|
||
def add(layer, col, row, res, source, rank=1):
|
||
# Clé de palier = (résolution, rang) : à résolution égale, les
|
||
# quadrants (rang 0) passent avant la dalle entière (rang 1) — même
|
||
# rendu, 4× moins de pixels à décoder.
|
||
records.setdefault(layer, {}).setdefault((col, row), {}) \
|
||
.setdefault((res, rank), []).append(source)
|
||
|
||
# 1. Dalles entières (palier le plus fin quand il est présent).
|
||
known = set(VIZ_LABELS)
|
||
for tile in scan_tiles(vis_dir):
|
||
col, row, res = tile["col"], tile["row"], tile["resolution"]
|
||
cell = _cell_bounds(col, row)
|
||
dir_path = Path(tile["dir_path"])
|
||
for viz_key, info in tile["viz"].items():
|
||
known.add(viz_key)
|
||
source_dir = dir_path
|
||
if info.get("dir_name") and info["dir_name"] != dir_path.name:
|
||
source_dir = dir_path.parent / info["dir_name"]
|
||
full = source_dir / info["filename"]
|
||
if full.is_file():
|
||
add(viz_key, col, row, res, _Source(full, cell, res))
|
||
|
||
# Clés les plus longues d'abord : `positive_openness` avant `openness`.
|
||
known = sorted(known, key=len, reverse=True)
|
||
|
||
# 2. Quadrants (index_subtiles) : AVIF pleine résolution + ses vignettes.
|
||
if sub_dir.is_dir():
|
||
quads = {}
|
||
for f in sub_dir.iterdir():
|
||
name = f.name
|
||
for suffix, tier in ((".avif", "full"), ("_mid.webp", "mid"),
|
||
(f"_thumb{_SUBTILE_THUMB_PX}.webp", "thumb")):
|
||
if not name.endswith(suffix):
|
||
continue
|
||
m = _SUBTILE_RE.search(name[:-len(suffix)] + ".avif")
|
||
if not m:
|
||
break
|
||
stem = name[:m.start()]
|
||
i, j = int(m.group(1)), int(m.group(2))
|
||
dir_name, viz = _split_viz_key(stem, known)
|
||
if dir_name is None:
|
||
break
|
||
quads.setdefault((dir_name, viz, tier), []).append((i, j, f))
|
||
break
|
||
for (dir_name, viz, tier), items in quads.items():
|
||
cell = _cell_of_dir(dir_name)
|
||
if cell is None:
|
||
continue
|
||
col, row, res = cell
|
||
k = max(max(i for i, _j, _f in items),
|
||
max(j for _i, j, _f in items)) + 1
|
||
step = 1000.0 / k
|
||
px = {"full": step / res, "mid": 640, "thumb": _SUBTILE_THUMB_PX}[tier]
|
||
tier_res = res if tier == "full" else step / px
|
||
base = _cell_bounds(col, row)
|
||
for i, j, f in items:
|
||
add(viz, col, row, tier_res, _Source(
|
||
f, (base[0] + i * step, base[1] + j * step,
|
||
base[0] + (i + 1) * step, base[1] + (j + 1) * step),
|
||
tier_res), rank=0)
|
||
|
||
# 3. Vignettes de dalle (paliers grossiers).
|
||
if thumb_dir.is_dir():
|
||
for f in thumb_dir.iterdir():
|
||
if f.suffix.lower() not in (".jpg", ".jpeg", ".webp", ".png"):
|
||
continue
|
||
stem = f.name[:-len(f.suffix)]
|
||
res_px = _MID_RES if stem.endswith("_mid") else _THUMB_RES
|
||
if stem.endswith("_mid"):
|
||
stem = stem[:-4]
|
||
dir_name, viz = _split_viz_key(stem, known)
|
||
if dir_name is None:
|
||
continue
|
||
cell = _cell_of_dir(dir_name)
|
||
if cell is None:
|
||
continue
|
||
col, row, _res = cell
|
||
add(viz, col, row, res_px, _Source(f, _cell_bounds(col, row), res_px))
|
||
|
||
# Paliers triés du plus grossier au plus fin (l'ordre de choix du rendu).
|
||
layers = {}
|
||
for viz, per_cell in records.items():
|
||
for (col, row), by_res in per_cell.items():
|
||
tiers = [by_res[k] for k in sorted(by_res, key=lambda k: (-k[0], k[1]))]
|
||
layers.setdefault(viz, {})[(col, row)] = tiers
|
||
return layers
|
||
|
||
|
||
def source_index(output_dir, force=False):
|
||
"""Index des sources, mémoïsé (TTL + mtime des dossiers surveillés)."""
|
||
output_dir = Path(output_dir)
|
||
key = str(output_dir)
|
||
stamp = tuple(_dir_mtime(output_dir / d) for d in
|
||
("visualisations", "index_thumbs", "index_subtiles"))
|
||
entry = _index_cache.get(key)
|
||
now = time.time()
|
||
if entry and not force and now - entry["at"] < _INDEX_TTL:
|
||
# Le TTL prime sur la mtime des dossiers : en mode amont, chaque source
|
||
# rapatriée la modifierait et provoquerait un rescan par tuile servie.
|
||
if REMOTE_SOURCE_URL or entry["stamp"] == stamp:
|
||
return entry["layers"]
|
||
if REMOTE_SOURCE_URL:
|
||
# L'amont fait autorité : il connaît toutes les dalles, le cache local
|
||
# n'en détient qu'une partie (et grossit à chaque source rapatriée —
|
||
# le rescanner à chaque tuile coûterait plus cher que le rendu). Les
|
||
# sources déjà présentes sont servies depuis le disque (_Source.ensure).
|
||
layers = dict(_remote_index(output_dir, force))
|
||
if not layers:
|
||
layers = _build_index(output_dir) # amont muet : cache local seul
|
||
else:
|
||
layers = _build_index(output_dir)
|
||
_index_cache[key] = {"layers": layers, "stamp": stamp, "at": now}
|
||
return layers
|
||
|
||
|
||
def available_layers(output_dir):
|
||
"""Couches présentes sur disque, ordonnées comme le panneau de la carte."""
|
||
from .index import _VIZ_FALLBACK_ORDER
|
||
found = set(source_index(output_dir))
|
||
ordered = [v for v in _VIZ_FALLBACK_ORDER if v in found]
|
||
ordered += sorted(found - set(ordered))
|
||
return ordered
|
||
|
||
|
||
def grid_bounds_l93(output_dir, layer=None):
|
||
"""Emprise L93 (min_x, min_y, max_x, max_y) des dalles disponibles."""
|
||
layers = source_index(output_dir)
|
||
cells = set()
|
||
for key, per_cell in layers.items():
|
||
if layer and key != layer:
|
||
continue
|
||
cells.update(per_cell)
|
||
if not cells:
|
||
return None
|
||
xs = [c for c, _r in cells]
|
||
ys = [r for _c, r in cells]
|
||
return (min(xs) * 1000.0, (min(ys) - 1) * 1000.0,
|
||
(max(xs) + 1) * 1000.0, max(ys) * 1000.0)
|
||
|
||
|
||
def grid_bounds_wgs84(output_dir, layer=None):
|
||
"""Emprise WGS84 [ouest, sud, est, nord] des dalles (TileJSON, WMTS)."""
|
||
l93 = grid_bounds_l93(output_dir, layer)
|
||
if l93 is None:
|
||
return None
|
||
min_x, min_y, max_x, max_y = l93
|
||
tr = _transformer("EPSG:2154", "EPSG:4326")
|
||
xs = [min_x, max_x, min_x, max_x]
|
||
ys = [min_y, min_y, max_y, max_y]
|
||
lons, lats = tr.transform(xs, ys)
|
||
return [min(lons), min(lats), max(lons), max(lats)]
|
||
|
||
|
||
def tiles_stamp(output_dir):
|
||
"""Version globale du jeu de tuiles (max des mtimes) pour l'URL du client."""
|
||
newest = 0.0
|
||
for per_cell in source_index(output_dir).values():
|
||
for tiers in per_cell.values():
|
||
for src in tiers[-1]:
|
||
m = src.mtime()
|
||
if m and m > newest:
|
||
newest = m
|
||
return int(newest * 1000)
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# Rendu d'une tuile
|
||
# ---------------------------------------------------------------------------
|
||
|
||
# Cache des images sources décodées : le décodage (AVIF surtout) domine le
|
||
# coût d'une tuile et les tuiles voisines partagent leurs dalles. Le budget est
|
||
# exprimé en OCTETS, pas en nombre d'entrées : une dalle 5000² pèse ~75 Mo
|
||
# quand une vignette en pèse 0,2 — un cache « N entrées » ferait déborder la
|
||
# mémoire d'une petite machine.
|
||
SOURCE_CACHE_BYTES = int(os.environ.get("LIDAR_TILE_SOURCE_CACHE_MB", "192")) * 1024 * 1024
|
||
_source_cache = OrderedDict()
|
||
_source_cache_lock = threading.Lock()
|
||
|
||
|
||
def _open_source(path_str, mtime):
|
||
"""Image source décodée, mémoïsée par (chemin, mtime).
|
||
|
||
`mtime` fait partie de la clé : une dalle régénérée invalide l'entrée.
|
||
"""
|
||
from PIL import Image
|
||
key = (str(path_str), mtime)
|
||
with _source_cache_lock:
|
||
img = _source_cache.get(key)
|
||
if img is not None:
|
||
_source_cache.move_to_end(key)
|
||
return img
|
||
img = Image.open(str(path_str))
|
||
img.load()
|
||
if img.mode not in ("RGB", "RGBA"):
|
||
img = img.convert("RGB")
|
||
size = img.size[0] * img.size[1] * len(img.getbands())
|
||
with _source_cache_lock:
|
||
_source_cache[key] = img
|
||
_source_cache[key].info["_bytes"] = size
|
||
total = sum(i.info.get("_bytes", 0) for i in _source_cache.values())
|
||
while total > SOURCE_CACHE_BYTES and len(_source_cache) > 1:
|
||
_k, old_img = _source_cache.popitem(last=False)
|
||
total -= old_img.info.get("_bytes", 0)
|
||
return img
|
||
|
||
|
||
def clear_source_cache():
|
||
"""Vide le cache d'images sources (tests, pression mémoire)."""
|
||
with _source_cache_lock:
|
||
_source_cache.clear()
|
||
|
||
|
||
def _pick_tier(tiers, target_res):
|
||
"""Palier le plus grossier dont la résolution suffit à la tuile visée."""
|
||
for group in tiers:
|
||
if group[0].res <= target_res:
|
||
return group
|
||
return tiers[-1]
|
||
|
||
|
||
def _contributing(output_dir, layer, z, x, y, scale):
|
||
"""Sources intersectant la tuile, palier choisi selon la résolution visée."""
|
||
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)
|
||
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):
|
||
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)
|
||
return out
|
||
|
||
|
||
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}")
|
||
return False
|
||
|
||
w, h = img.size
|
||
min_x, min_y, max_x, max_y = src.bounds
|
||
res_x = (max_x - min_x) / w
|
||
res_y = (max_y - min_y) / h
|
||
|
||
# Fenêtre source utile = intersection avec l'emprise L93 de la tuile,
|
||
# élargie de 2 px pour que l'interpolation dispose de son voisinage.
|
||
t_min_x, t_min_y, t_max_x, t_max_y = tile_bounds_l93(z, x, y)
|
||
c0 = max(0, int(math.floor((max(t_min_x, min_x) - min_x) / res_x)) - 2)
|
||
c1 = min(w, int(math.ceil((min(t_max_x, max_x) - min_x) / res_x)) + 2)
|
||
r0 = max(0, int(math.floor((max_y - min(t_max_y, max_y)) / res_y)) - 2)
|
||
r1 = min(h, int(math.ceil((max_y - max(t_min_y, min_y)) / res_y)) + 2)
|
||
if c1 - c0 < 1 or r1 - r0 < 1:
|
||
return False
|
||
|
||
crop = img.crop((c0, r0, c1, r1))
|
||
if crop.mode != "RGBA":
|
||
crop = crop.convert("RGBA")
|
||
cw, ch = crop.size
|
||
|
||
# Coins L93 de la fenêtre découpée → EPSG:3857 → pixels de la tuile.
|
||
wx0 = min_x + c0 * res_x
|
||
wx1 = min_x + c1 * res_x
|
||
wy1 = max_y - r0 * res_y
|
||
wy0 = max_y - r1 * res_y
|
||
mx, my = to_3857([wx0, wx1, wx1, wx0], [wy1, wy1, wy0, wy0])
|
||
west, south, east, north = tile_bounds_3857(z, x, y)
|
||
px = [(m - west) / (east - west) * size for m in mx]
|
||
py = [(north - m) / (north - south) * size for m in my]
|
||
dst_quad = list(zip(px, py)) # NW, NE, SE, SW
|
||
src_quad = [(0, 0), (cw, 0), (cw, ch), (0, ch)]
|
||
coeffs = perspective_coeffs(dst_quad, src_quad)
|
||
if coeffs is None:
|
||
return False
|
||
|
||
warped = crop.transform((size, size), Image.PERSPECTIVE, coeffs,
|
||
resample=resample, fillcolor=(0, 0, 0, 0))
|
||
canvas.alpha_composite(warped)
|
||
return True
|
||
|
||
|
||
def render_tile(output_dir, layer, z, x, y, scale=1):
|
||
"""Rend une tuile en mémoire. Retourne une image RGBA, ou None si vide."""
|
||
from PIL import Image
|
||
|
||
sources = _contributing(output_dir, layer, z, x, y, scale)
|
||
if not sources:
|
||
return None
|
||
size = TILE_SIZE * scale
|
||
canvas = Image.new("RGBA", (size, size), (0, 0, 0, 0))
|
||
# Agrandissement (zoom natif) : bicubique ; réduction : bilinéaire suffit
|
||
# puisque le palier source est déjà calé sur la résolution de la tuile.
|
||
resample = Image.BICUBIC
|
||
painted = False
|
||
for src in sources:
|
||
painted |= _paste_source(canvas, src, z, x, y, size, resample)
|
||
if not painted:
|
||
return None
|
||
return canvas
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# Cache disque
|
||
# ---------------------------------------------------------------------------
|
||
|
||
def tile_cache_path(output_dir, layer, z, x, y, scale=1, fmt="png"):
|
||
"""Chemin du fichier de cache d'une tuile."""
|
||
suffix = "@2x" if scale == 2 else ""
|
||
return (Path(output_dir) / TILE_DIRNAME / layer / str(z) / str(x)
|
||
/ f"{y}{suffix}.{fmt}")
|
||
|
||
|
||
def _encode(img, fmt):
|
||
import io
|
||
buf = io.BytesIO()
|
||
if fmt == "png":
|
||
if PNG_PALETTE:
|
||
from PIL import Image
|
||
img = img.quantize(colors=255, method=Image.Quantize.FASTOCTREE)
|
||
img.save(buf, format="PNG", optimize=False)
|
||
elif fmt == "webp":
|
||
img.save(buf, format="WEBP", quality=WEBP_QUALITY)
|
||
elif fmt == "avif":
|
||
img.save(buf, format="AVIF", quality=AVIF_QUALITY)
|
||
else:
|
||
raise ValueError(f"format de tuile inconnu : {fmt}")
|
||
return buf.getvalue()
|
||
|
||
|
||
def _write_atomic(path, data):
|
||
path.parent.mkdir(parents=True, exist_ok=True)
|
||
tmp = path.with_name(path.name + f".{os.getpid()}.part")
|
||
try:
|
||
tmp.write_bytes(data)
|
||
os.replace(tmp, path)
|
||
except OSError as e:
|
||
logger.debug(f"Écriture de tuile impossible ({path}) : {e}")
|
||
tmp.unlink(missing_ok=True)
|
||
|
||
|
||
def zoom_supported(z, scale=1):
|
||
"""Vrai si le zoom demandé est dans la plage rendue par le serveur."""
|
||
max_z = TILE_MAX_NATIVE_Z - (1 if scale == 2 else 0)
|
||
return TILE_MIN_Z <= z <= max_z
|
||
|
||
|
||
def get_tile(output_dir, layer, z, x, y, scale=1, fmt="png"):
|
||
"""Tuile encodée (bytes) depuis le cache, rendue au besoin.
|
||
|
||
Retourne None quand aucune donnée n'intersecte la tuile ; un marqueur
|
||
`.empty` mémorise ce cas pour ne jamais le recalculer. Une tuile en cache
|
||
est réutilisée tant qu'aucune dalle contributrice n'est plus récente
|
||
qu'elle : une dalle régénérée n'invalide que ses propres tuiles.
|
||
"""
|
||
cache = tile_cache_path(output_dir, layer, z, x, y, scale, fmt)
|
||
empty = cache.with_suffix(".empty")
|
||
sources = _contributing(output_dir, layer, z, x, y, scale)
|
||
if not sources:
|
||
if not empty.exists():
|
||
_write_atomic(empty, b"")
|
||
return None
|
||
|
||
newest = 0.0
|
||
for src in sources:
|
||
m = src.mtime()
|
||
if m and m > newest:
|
||
newest = m
|
||
try:
|
||
if cache.stat().st_mtime >= newest:
|
||
return cache.read_bytes()
|
||
except OSError:
|
||
pass
|
||
|
||
img = render_tile(output_dir, layer, z, x, y, scale)
|
||
if img is None:
|
||
_write_atomic(empty, b"")
|
||
return None
|
||
data = _encode(img, fmt)
|
||
_write_atomic(cache, data)
|
||
empty.unlink(missing_ok=True)
|
||
return data
|
||
|
||
|
||
def empty_marker_exists(output_dir, layer, z, x, y, scale=1, fmt="png"):
|
||
"""Vrai si la tuile est déjà connue comme vide (négatif mémorisé)."""
|
||
return tile_cache_path(output_dir, layer, z, x, y,
|
||
scale, fmt).with_suffix(".empty").exists()
|
||
|
||
|
||
def transparent_tile(scale=1, fmt="png"):
|
||
"""Tuile entièrement transparente (zone sans donnée)."""
|
||
from PIL import Image
|
||
size = TILE_SIZE * scale
|
||
return _encode(Image.new("RGBA", (size, size), (0, 0, 0, 0)), fmt)
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# Pré-chauffage
|
||
# ---------------------------------------------------------------------------
|
||
|
||
def tiles_in_bounds(bounds_wgs84, z):
|
||
"""Indices (x, y) des tuiles du niveau z couvrant une bbox WGS84."""
|
||
west, south, east, north = bounds_wgs84
|
||
n = 2 ** z
|
||
|
||
def _x(lon):
|
||
return int((lon + 180.0) / 360.0 * n)
|
||
|
||
def _y(lat):
|
||
lat = max(-85.0511, min(85.0511, lat))
|
||
rad = math.radians(lat)
|
||
return int((1.0 - math.log(math.tan(rad) + 1 / math.cos(rad)) / math.pi) / 2.0 * n)
|
||
|
||
x0, x1 = sorted((_x(west), _x(east)))
|
||
y0, y1 = sorted((_y(north), _y(south)))
|
||
return [(x, y)
|
||
for x in range(max(0, x0), min(n - 1, x1) + 1)
|
||
for y in range(max(0, y0), min(n - 1, y1) + 1)]
|
||
|
||
|
||
def warm(output_dir, layers, z_min, z_max, bounds_wgs84=None, scale=1,
|
||
fmt="png", limit=20000):
|
||
"""Pré-calcule les tuiles d'une emprise (retourne le compte rendu)."""
|
||
if bounds_wgs84 is None:
|
||
bounds_wgs84 = grid_bounds_wgs84(output_dir)
|
||
if bounds_wgs84 is None:
|
||
return {"rendues": 0, "vides": 0, "limite": False}
|
||
done = empty = 0
|
||
for layer in layers:
|
||
for z in range(max(TILE_MIN_Z, z_min), min(TILE_MAX_NATIVE_Z, z_max) + 1):
|
||
if not zoom_supported(z, scale):
|
||
continue
|
||
for x, y in tiles_in_bounds(bounds_wgs84, z):
|
||
if done + empty >= limit:
|
||
return {"rendues": done, "vides": empty, "limite": True}
|
||
if get_tile(output_dir, layer, z, x, y, scale, fmt) is None:
|
||
empty += 1
|
||
else:
|
||
done += 1
|
||
return {"rendues": done, "vides": empty, "limite": False}
|