"""XYZ tile pyramid (EPSG:3857) rendered on demand from the tiles. Same tiling scheme as OpenStreetMap / Google Maps: Web Mercator grid, origin at the north-west corner, `{z}/{x}/{y}`, 256 px tiles (512 px with `scale=2`, `@2x` convention). The pipeline outputs remain 1 km Lambert 93 tiles: each XYZ tile is composed on the fly by reprojecting the source tiles that intersect it (one projective transform per source tile, error well below a pixel), then cached on disk. Deliberately free of GDAL and numpy: Pillow + pyproj are enough, so the lightweight image (Dockerfile.maps) stays small and ARM64-portable. """ 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") # --- Tiling contract (see docs/MAPS.md) ------------------------------------ TILE_SIZE = 256 # canonical size (OSM/XYZ); @2x → 512 TILE_MIN_Z = 5 TILE_MAX_NATIVE_Z = 19 # 0.2 m/px ≈ z19 resolution at latitude 47° TILE_DIRNAME = "index_xyz" # disk cache, under the output directory # Stored levels, in standard OSM numbering (256 px tile; an @2x tile at # level z equals a 256 px tile at z+1). By default ALL levels up to native are # written to disk and generated ahead of time by the background maintenance: # on-the-fly rendering on a Pi (~170 ms per tile, 3 at a time) took 1.5 to 6 s # per screen at high zoom. Reduced storage (small machine, limited disk): # lower LIDAR_TILE_CACHE_MAX_Z and/or LIDAR_TILE_EVEN_LEVELS=1 (even levels # only, the others rendered on the fly with a bounded memory cache). TILE_CACHE_MAX_Z = int(os.environ.get("LIDAR_TILE_CACHE_MAX_Z", str(TILE_MAX_NATIVE_Z)) or TILE_MAX_NATIVE_Z) TILE_EVEN_LEVELS = os.environ.get("LIDAR_TILE_EVEN_LEVELS", "0").strip() == "1" WEBP_QUALITY = 78 AVIF_QUALITY = 60 AVIF_SPEED = 9 # fast encoding (see rendering.AVIF_SPEED: ×7, +3 % size) # Palettized PNG (PNG8 + alpha): ~5× lighter (170 → 32 KB on a real tile) for # a mean error of ~4 levels on a color ramp. DISABLED by default: the # canonical PNG stays lossless, fidelity comes before bandwidth for an # interpretation product. `LIDAR_TILE_PNG_PALETTE=1` enables it when # bandwidth matters (mobile viewing). PNG_PALETTE = os.environ.get("LIDAR_TILE_PNG_PALETTE", "") == "1" # Layers made of flat coded levels, resampled with nearest neighbour NEAREST_LAYERS = frozenset({'densite_sol'}) # Equatorial half-circumference: extent of Web Mercator (EPSG:3857). ORIGIN = 20037508.342789244 # Nominal resolutions of the source tiers reused as-is from the pipeline # (m/px): thumbnail 256 px/km, intermediate thumbnail 640 px/km. _THUMB_RES = 1000.0 / 256 _MID_RES = 1000.0 / 640 _SUBTILE_RE = re.compile(r"_(\d+)_(\d+)\.avif$") # Upstream tile server (the worker's mapserve): when set, the source index # comes from its /api/tiles and missing images are fetched on demand into the # local cache — the map container then needs no local data at startup. 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} # Concurrent source downloads (network) — modest by default: on a small # server, each fetched source is a file of several MB. 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() # Source circuit breaker: after a failure, no further attempt for # _SOURCE_OFFLINE_S — without it, every source in a fetch queue would pay the # network timeout of a powered-off upstream again and the maintenance would # render nothing (the map must stay self-sufficient). _SOURCE_OFFLINE_S = 60.0 _SOURCE_OFFLINE = {"until": 0.0} # Source index rebuilt at most every _INDEX_TTL seconds (or as soon as a # directory's mtime changes: tile added/removed). _INDEX_TTL = 20.0 _index_cache = {} # --------------------------------------------------------------------------- # Grid geometry # --------------------------------------------------------------------------- def tile_bounds_3857(z, x, y): """Extent (west, south, east, north) of an XYZ tile in EPSG:3857 metres.""" 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): """Tile resolution at the equator (m/px); ×cos(lat) on the ground.""" return 2.0 * ORIGIN / (TILE_SIZE * scale * (2 ** z)) def tile_latitude(z, y): """Latitude (degrees) of a tile's centre — used to match resolution.""" 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): """Ground resolution targeted by the tile (m/px), latitude included.""" 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 (coordinate lists).""" return _transformer("EPSG:3857", "EPSG:2154").transform(xs, ys) def to_3857(xs, ys): """EPSG:2154 → EPSG:3857 (coordinate lists).""" return _transformer("EPSG:2154", "EPSG:3857").transform(xs, ys) @lru_cache(maxsize=4096) def wgs84_to_l93(lon, lat): """WGS84 point → Lambert 93 (used by the tile information card).""" return _transformer("EPSG:4326", "EPSG:2154").transform(lon, lat) @lru_cache(maxsize=4096) def tile_bounds_l93(z, x, y, samples=5): """L93 extent enclosing an XYZ tile. A tile's edges are not straight lines in Lambert 93: a `samples`×`samples` grid is sampled rather than the corners alone, otherwise the extent is underestimated at low zoom (tiles spanning hundreds of 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) # --------------------------------------------------------------------------- # Projective transform (output → source mapping, Pillow convention) # --------------------------------------------------------------------------- def _solve(matrix, rhs): """Solve a dense linear system (partial pivoting), without 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): """Pillow `Image.PERSPECTIVE` coefficients mapping output → source. Pillow samples the SOURCE at (x', y') = ((a x + b y + c) / (g x + h y + 1), (d x + e y + f) / (g x + h y + 1)) for each pixel (x, y) of the OUTPUT: the 8 unknowns are therefore solved from 4 correspondences (output point → source point). Returns None if the quadrilateral is degenerate (tile shrunk to a point at extreme zoom-out). """ 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) # --------------------------------------------------------------------------- # Source index: pipeline tiles, per layer and per tier # --------------------------------------------------------------------------- class _Source: """A georeferenced source image (whole tile or quadrant). Non-null `url`: the image lives on the upstream tile server and is only fetched on first use (`ensure()`), into `path`. `version` is then the mtime (ms) announced by the upstream: it serves as the reference date for tile staleness, even before any download. """ __slots__ = ("path", "bounds", "res", "url", "version", "seen") def __init__(self, path, bounds, res, url=None, version=None): self.path = path self.bounds = bounds # (min_x, min_y, max_x, max_y) in L93 self.res = res # nominal resolution (m/px) self.url = url self.version = version self.seen = 0.0 # first appearance of the tile in the inventory def mtime(self): """Reference date for tile staleness: the source's version, pushed forward to its first appearance in the inventory (a tile written before, but inventoried after, an XYZ tile makes that XYZ tile stale).""" if self.version is not None: base = self.version / 1000.0 else: try: base = self.path.stat().st_mtime except OSError: return None return max(base, self.seen) def ensure(self): """Ensure the image is present locally (fetch it if needed).""" if self.url is None or self.path.is_file(): return self.path.is_file() if not _fetch_source(self.url, self.path): return False if self.version is not None: # Dated at its upstream version: the local index (upstream down) # finds the same date and already-rendered tiles stay fresh. try: os.utime(self.path, (self.version / 1000.0, self.version / 1000.0)) except OSError: pass return True def _remote_thumb_px(k): """Side (px) of the thumbnail served by the upstream: 256 per tile, 160 per quadrant.""" return 160.0 if k > 1 else 256.0 def _fetch_source(url, dest): """Download a source image from the upstream (atomic write).""" import urllib.request with _fetch_guard: lock = _fetch_locks.setdefault(str(dest), threading.Lock()) with lock: if dest.is_file(): return True if time.time() < _SOURCE_OFFLINE["until"]: return False 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 — upstream down: partial tile logger.debug(f"Upstream source unavailable ({url}): {e}") _SOURCE_OFFLINE["until"] = time.time() + _SOURCE_OFFLINE_S return False _write_atomic(dest, data) logger.info(f"Source fetched: {dest.name} ({len(data) / 1e6:.1f} MB)") return dest.is_file() def _remote_payload(force=False): """Index from the upstream tile server (/api/tiles), cached for 60 s.""" import json import urllib.request now = time.time() # Failures are remembered too (TTL): with the upstream down and no known # inventory, at most one call per TTL — not a network timeout on every # map request. if not force and _remote_cache["at"] \ 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=10) as r: payload = json.loads(r.read().decode("utf-8")) except Exception as e: # noqa: BLE001 — keep the last known index logger.warning(f"Upstream index unreachable ({REMOTE_SOURCE_URL}): {e}") payload = _remote_cache["payload"] _remote_cache.update(payload=payload, at=now) return payload _SAFE_BASENAME = re.compile(r"^[A-Za-z0-9_.-]+$") def _persist_remote_quality(output_dir, payload): """Copy the upstream quality table into local sidecars. The quality inset of the PDF export always reads the local disk: the self-sufficient map (upstream down) thus keeps the last known values. Returns: number of sidecars (re)written. """ from .quality import QUALITY_VERSION, write_quality written = 0 for base, data in (payload.get("quality") or {}).items(): if (not isinstance(base, str) or not _SAFE_BASENAME.match(base) or base.startswith(".") or not isinstance(data, dict) or data.get("version") != QUALITY_VERSION): continue try: written += bool(write_quality(output_dir, base, data)) except OSError as e: logger.debug(f"Quality sidecar not copied ({base}): {e}") return written def _remote_index(output_dir, force=False): """Inventory built from the upstream: sources not fetched yet. The result is memoized alongside the payload: rebuilding it costs hundreds of milliseconds on a catalogue of several thousand tiles, which would be paid on EVERY tile served. """ 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) _persist_remote_quality(output_dir, payload) 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): """L93 extent of an LHD tile: 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 cell_bounds_wgs84(col, row): """WGS84 extent [west, south, east, north] of an LHD tile (4 L93 corners).""" min_x, min_y, max_x, max_y = _cell_bounds(col, row) xs = (min_x, min_x, max_x, max_x) ys = (min_y, max_y, min_y, max_y) lons, lats = _transformer("EPSG:2154", "EPSG:4326").transform(xs, ys) return (min(lons), min(lats), max(lons), max(lats)) 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), relying on the known keys. Visualization keys contain underscores (`positive_openness`, `hillshade_multi`): the stem cannot be split at the last `_`, so a known suffix is matched instead (longest first). """ 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, resolution) of a tile directory name, or 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): """Inventory {layer: {(col, row): [tiers from coarsest to finest]}}. The three pipeline sources are scanned INDEPENDENTLY — whole tiles (`visualisations/`), quadrants (`index_subtiles/`) and thumbnails (`index_thumbs/`). A partial cache therefore remains usable: on a lightweight machine, only quadrants and thumbnails are fetched, never whole tiles. """ # _SUBTILE_THUMB_PX: size of the quadrant thumbnails, defined by the # index (changing it there must have no effect here). 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, per layer and per tile; sorted into tiers at the end. records = {} def add(layer, col, row, res, source, rank=1): # Tier key = (resolution, rank): at equal resolution, quadrants # (rank 0) come before the whole tile (rank 1) — same rendering, # 4× fewer pixels to decode. records.setdefault(layer, {}).setdefault((col, row), {}) \ .setdefault((res, rank), []).append(source) # 1. Whole tiles (finest tier when present). 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)) # Longest keys first: `positive_openness` before `openness`. known = sorted(known, key=len, reverse=True) # 2. Quadrants (index_subtiles): full-resolution AVIF (lossless WebP for # flat-level layers) + their thumbnails. if sub_dir.is_dir(): quads = {} for f in sub_dir.iterdir(): name = f.name # Plain ".webp" last: thumbnails also end in .webp for suffix, tier in ((".avif", "full"), ("_mid.webp", "mid"), (f"_thumb{_SUBTILE_THUMB_PX}.webp", "thumb"), (".webp", "full")): 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. Tile thumbnails (coarse tiers). 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)) # Tiers sorted from coarsest to finest (the rendering's order of choice). 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 _SEEN_FILE = ".sources_seen.json" _seen_cache = {} _seen_lock = threading.Lock() def _apply_seen(output_dir, layers): """Date each source tile by its first appearance in the inventory. A source tile can enter the inventory AFTER an XYZ tile covering it was rendered while carrying an older date (written before, inventoried after: upstream index TTL, inventory debounce). Comparing version dates alone would leave the XYZ tile "fresh" — a permanent hole at that level, all the way to the browser (unchanged stamp). Persistent registry (index_xyz/.sources_seen.json): a restart makes nothing stale; on the very first inventory, nothing is dated (everything still has to be rendered). """ import json path = Path(output_dir) / TILE_DIRNAME / _SEEN_FILE with _seen_lock: seen = _seen_cache.get(str(path)) first = False if seen is None: try: seen = {k: float(v) for k, v in json.loads(path.read_text(encoding="utf-8")).items()} except (OSError, ValueError, AttributeError): # No registry: very first inventory (nothing to invalidate), # unless a tile cache already exists (upgrade) — its tiles may # have been rendered before some source tiles arrived: they are # all made stale once (the old one keeps being served meanwhile). cache_dir = path.parent upgraded = cache_dir.is_dir() and any( p.is_dir() for p in cache_dir.iterdir()) seen, first = {}, not upgraded _seen_cache[str(path)] = seen now = time.time() added = False for layer, per_cell in layers.items(): for (col, row), tiers in per_cell.items(): key = f"{layer}|{col}|{row}" at = seen.get(key) if at is None: at = seen[key] = 0.0 if first else now added = True if at: for group in tiers: for src in group: src.seen = at if added or first: try: _write_atomic(path, json.dumps(seen).encode("utf-8")) except Exception as e: # noqa: BLE001 — registry lost: staleness by version only logger.debug(f"Tile registry not written ({path}): {e}") def source_index(output_dir, force=False): """Source index, memoized (TTL + mtime of the watched directories).""" 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: # The TTL takes precedence over the directories' mtime: in upstream # mode, each fetched source would change it and trigger a rescan per # tile served. if REMOTE_SOURCE_URL or entry["stamp"] == stamp: return entry["layers"] if REMOTE_SOURCE_URL: # The upstream is authoritative: it knows every tile, the local cache # only holds part of them (and grows with each fetched source — # rescanning it on every tile would cost more than rendering). Sources # already present are served from disk (_Source.ensure). layers = dict(_remote_index(output_dir, force)) if not layers: layers = _build_index(output_dir) # silent upstream: local cache only else: layers = _build_index(output_dir) _apply_seen(output_dir, layers) _index_cache[key] = {"layers": layers, "stamp": stamp, "at": now} return layers def available_layers(output_dir): """Layers present on disk and displayed (PANEL_VIZ), ordered like the map panel.""" from . import index as index_mod from .index import _VIZ_FALLBACK_ORDER found = set(source_index(output_dir)) if index_mod.PANEL_VIZ is not None: found &= set(index_mod.PANEL_VIZ) 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): """L93 extent (min_x, min_y, max_x, max_y) of the available tiles.""" 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): """WGS84 extent [west, south, east, north] of the tiles (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): """Global version of the tile set (max of the finest-tier mtimes) for the client URL.""" 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) # --------------------------------------------------------------------------- # Rendering a tile # --------------------------------------------------------------------------- # Cache of decoded source images: decoding (AVIF above all) dominates the # cost of a tile, and neighbouring tiles share their sources. The budget is # expressed in BYTES, not in number of entries: a 5000² tile weighs ~100 MB # (4 bytes/px) while a thumbnail weighs 0.26 MB — an "N entries" cache would # overflow the memory of a small 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): """Decoded source image, memoized by (path, mtime). `mtime` is part of the key: a regenerated tile invalidates the entry. """ 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 # Copy detached from the file: an open AVIF image keeps its decoder # (libavif/dav1d buffers, ~18 MB per 2500² quadrant) for as long as it # lives — cached, this doubled the real memory footprint and got the Pi's # container (1 GB) killed by the OOM killer when browsing at high zoom. with Image.open(str(path_str)) as opened: opened.load() img = (opened.copy() if opened.mode in ("RGB", "RGBA") else opened.convert("RGB")) # PIL stores RGB like RGBA: 4 bytes per pixel in both cases. size = img.size[0] * img.size[1] * 4 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(): """Empty the source image cache (tests, memory pressure).""" with _source_cache_lock: _source_cache.clear() def sources_to_keep(tiers): """Tiers a self-sufficient machine must keep to render ALL its levels without the upstream: from the coarsest up to the first tier that reaches the finest resolution (a whole tile at the same resolution as its quadrants would be a duplicate).""" if not tiers: return [] finest = min(group[0].res for group in tiers) keep = [] for group in tiers: keep.extend(group) if group[0].res <= finest: break return keep def _pick_tier(tiers, target_res): """Coarsest tier whose resolution is sufficient for the target tile.""" for group in tiers: if group[0].res <= target_res: return group return tiers[-1] def sources_in_bbox(output_dir, layer, bbox, target_res): """Sources of a layer intersecting an L93 extent, at the coarsest tier sufficient for target_res (m/px). Returns [((col, row), source)].""" per_cell = source_index(output_dir).get(layer) if not per_cell: return [] 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_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(((col, row), src)) return out def _contributing(output_dir, layer, z, x, y, scale): """Sources intersecting the tile, tier chosen from the target resolution.""" 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): """Decoded image of a source (fetched if needed), or None if unreadable.""" 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 — unreadable source: partial rendering logger.debug(f"Unreadable source ({src.path.name}): {e}") return None def _paste_source(canvas, src, z, x, y, size, resample): """Reproject a source into the tile (projective transform).""" from PIL import Image img = load_source(src) if img is None: 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 # Useful source window = intersection with the tile's L93 extent, widened # by 2 px so the interpolation has its neighbourhood. 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 # L93 corners of the cropped window → EPSG:3857 → tile pixels. 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): """Render a tile in memory. Returns an RGBA image, or None if empty.""" 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)) # Bicubic in both directions: needed when upscaling (native zoom), and # harmless when downscaling since the source tier already matches the # tile's resolution. # Flat coded levels (density, 1 m/px): nearest neighbour, otherwise # upscaling at zooms 18–19 invents intermediate greys. resample = Image.NEAREST if layer in NEAREST_LAYERS else 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 # --------------------------------------------------------------------------- # Disk cache # --------------------------------------------------------------------------- def tile_cache_path(output_dir, layer, z, x, y, scale=1, fmt="png"): """Path of a tile's cache file.""" 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, speed=AVIF_SPEED) else: raise ValueError(f"unknown tile format: {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) # Dated with the precise clock (time.time), the one used for the # tiles' first-appearance dates (_apply_seen). The kernel's coarse # clock (default mtime) lags by a few ms — a tile written right after # a source tile arrived would otherwise look older than it. now = time.time() os.utime(tmp, (now, now)) os.replace(tmp, path) except OSError as e: logger.debug(f"Cannot write tile ({path}): {e}") tmp.unlink(missing_ok=True) def zoom_cached(z, scale=1): """True if tiles at this level are stored on disk.""" zc = z + (1 if scale == 2 else 0) if zc > TILE_CACHE_MAX_Z: return False return not (TILE_EVEN_LEVELS and zc % 2) # Tiles rendered on the fly (unstored levels): bounded memory cache, so that # panning back and forth does not recompute everything. MEMORY_CACHE_BYTES = int(os.environ.get("LIDAR_TILE_MEMORY_CACHE_MB", "64")) * 1024 * 1024 _mem_tiles = OrderedDict() _mem_tiles_bytes = [0] _mem_tiles_lock = threading.Lock() def _mem_get(key): with _mem_tiles_lock: data = _mem_tiles.get(key) if data is not None: _mem_tiles.move_to_end(key) return data def _mem_put(key, data): with _mem_tiles_lock: if key in _mem_tiles: return _mem_tiles[key] = data _mem_tiles_bytes[0] += len(data) while _mem_tiles_bytes[0] > MEMORY_CACHE_BYTES and _mem_tiles: _, old = _mem_tiles.popitem(last=False) _mem_tiles_bytes[0] -= len(old) def zoom_supported(z, scale=1): """True if the requested zoom is within the range rendered by the server.""" 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"): """Encoded tile (bytes) from the cache, rendered when needed. Returns None when no data intersects the tile; an `.empty` marker records this case so it is never recomputed (stored levels only). A cached tile is reused as long as no contributing source tile is newer than it: a regenerated source tile only invalidates its own XYZ tiles. """ cache = tile_cache_path(output_dir, layer, z, x, y, scale, fmt) empty = cache.with_suffix(".empty") stored = zoom_cached(z, scale) sources = _contributing(output_dir, layer, z, x, y, scale) if not sources: if stored and 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 if not stored: # Unstored level: rendered on the fly, memory cache only. key = (str(output_dir), layer, z, x, y, scale, fmt, newest) data = _mem_get(key) if data is None: img = render_tile(output_dir, layer, z, x, y, scale) if img is None: return None data = _encode(img, fmt) _mem_put(key, data) return data 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"): """True if the tile is already known to be empty (remembered negative).""" return tile_cache_path(output_dir, layer, z, x, y, scale, fmt).with_suffix(".empty").exists() def cached_tile(output_dir, layer, z, x, y, scale=1, fmt="png"): """Tile from the cache ONLY — never rendered (cache-only mode). Returns ``(data, state)``: - ``"fresh"`` : tile cached and newer than its source tiles; - ``"empty"`` : no source intersects it (`.empty` marker set); - ``"pending"``: data present but tile missing or stale — to be rendered elsewhere (background task, warm-up), not while browsing. """ 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, "empty" newest = 0.0 for src in sources: m = src.mtime() if m and m > newest: newest = m # Same freshness rule as for the tile: an .empty marker older than a # source (tile appeared since, newer upstream version) is invalidated — # otherwise an area browsed before it was generated would stay # transparent forever. if empty.exists(): try: if newest and empty.stat().st_mtime >= newest: return None, "empty" except OSError: return None, "empty" empty.unlink(missing_ok=True) try: if cache.stat().st_mtime >= newest: return cache.read_bytes(), "fresh" except OSError: pass return None, "pending" def transparent_tile(scale=1, fmt="png"): """Fully transparent tile (area without data).""" from PIL import Image size = TILE_SIZE * scale return _encode(Image.new("RGBA", (size, size), (0, 0, 0, 0)), fmt) # --------------------------------------------------------------------------- # Warm-up # --------------------------------------------------------------------------- def tiles_in_bounds(bounds_wgs84, z): """Indices (x, y) of the level-z tiles covering a WGS84 bbox.""" 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): """Pre-render the tiles of an extent (returns the report; keys `rendues`/`vides`/`limite` are part of the API payload).""" 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}