fb892ea9f2
Comments, docstrings, logs, CLI help, map UI, legends, PDF sheet, scripts, compose files and AGENTS.md are now English. Data keys stay unchanged (relief_oriente, densite_sol, visualisations/, API JSON keys, link params). Wrong comments and help defaults found along the way are corrected. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
1247 lines
53 KiB
Python
1247 lines
53 KiB
Python
"""Catalog of processed tiles: shared registries, thumbnails, inventory.
|
||
|
||
This module no longer holds any user interface (the map UI lives in
|
||
mapserve.py/mapui.py and web/map.{html,css,js}, image lidar-maps). It produces
|
||
the artifacts the pipeline and the map need:
|
||
|
||
- shared registries: VIZ_LABELS/VIZ_LEGENDS (labels, legends), display
|
||
defaults (DEFAULT_VIZ, PRECISION_VIZ, VIEW_MODES), output keyword ↔
|
||
pipeline step mappings (KEYWORD_TO_STEP);
|
||
- thumbnails (index_thumbs/) and sub-tiles (index_subtiles/) used as source
|
||
levels of the XYZ pyramid (tiles.py);
|
||
- inventory output/index_tiles.json: tiles, layers and versioned URLs,
|
||
served by mapserve's /api/tiles to the lightweight machines
|
||
(LIDAR_SOURCE_URL) and rebuilt after every tile by default (unless
|
||
--no-index).
|
||
|
||
Integration:
|
||
- called automatically at the end of process_all() in pipeline.py;
|
||
- standalone rebuild via --rebuild-index in cli.py.
|
||
"""
|
||
|
||
import json
|
||
import logging
|
||
import re
|
||
import time
|
||
from datetime import datetime
|
||
from pathlib import Path
|
||
|
||
logger = logging.getLogger("lidar")
|
||
|
||
|
||
# Display names of the layers (map panel, inventory, TileJSON, WMTS, JOSM).
|
||
# Key = keyword in the output file name (after the basename).
|
||
VIZ_LABELS = {
|
||
'hillshade_multi': 'Multidirectional hillshade',
|
||
'slope': 'Slope',
|
||
'aspect': 'Aspect',
|
||
'mslrm': 'MSRM (multi-scale relief)',
|
||
'sailore': 'SAILORE (adaptive LRM)',
|
||
'positive_openness': 'Positive openness',
|
||
'negative_openness': 'Negative openness',
|
||
'svf': 'Sky-View Factor',
|
||
'roughness': 'Roughness',
|
||
'wavelet': 'Wavelet',
|
||
'flow_acc': 'Flow accumulation',
|
||
'solar': 'Solar illumination',
|
||
'anomaly': 'Anomaly map',
|
||
'relief_oriente': 'Oriented relief',
|
||
'densite_sol': 'Precision (ground point density)',
|
||
'ortho': 'IGN orthophoto',
|
||
'topo': 'IGN topographic map',
|
||
}
|
||
|
||
# Layer legends: title, how to read the rendering (colors), computation method
|
||
# and colormap gradient. Single source shared by rendering.py (merged into
|
||
# COLORMAPS), mapserve.py (/api/map/meta, TileJSON) and export_pdf.py (PDF
|
||
# legend) — this module deliberately has no heavy dependency (the lightweight
|
||
# lidar-maps image has neither matplotlib nor GDAL).
|
||
# 'gradient': 9 stops sampled from the matplotlib colormap (plt.get_cmap at
|
||
# i/8) to draw a gradient bar without matplotlib; it must be kept in sync by
|
||
# hand with the 'cmap' of rendering.COLORMAPS (no test checks it).
|
||
# 'ticks': labels of the gradient ends (None = no bar).
|
||
# 'reading' (optional): "How to read" sentences shown by the map and the PDF.
|
||
VIZ_LEGENDS = {
|
||
'hillshade_multi': {
|
||
'title': 'Multidirectional Hillshade',
|
||
'legend': 'Combined illumination from 8 directions (fixed 0–1 scale)\nWhite = lit face | Black = shadow\nConsistent colors across tiles',
|
||
'description': 'Cast shadows revealing micro-relief (walls, ditches, terraces)',
|
||
'cmap': 'gray',
|
||
'gradient': ('#000000', '#202020', '#404040', '#606060', '#808080',
|
||
'#a0a0a0', '#c0c0c0', '#e0e0e0', '#ffffff'),
|
||
'ticks': ('0', '1'),
|
||
},
|
||
'slope': {
|
||
'title': 'Slope (terrain steepness)',
|
||
'legend': 'Steepness in degrees\nFixed 0–30° scale — consistent colors across tiles\nYellow = steep slope | Dark purple = flat ground',
|
||
'description': 'Walls, banks and edges stand out in yellow — flat ground is dark',
|
||
'cmap': 'inferno',
|
||
'gradient': ('#000004', '#210c4a', '#57106e', '#8a226a', '#bc3754',
|
||
'#e45a31', '#f98e09', '#f9cb35', '#fcffa4'),
|
||
'ticks': ('0°', '30°'),
|
||
},
|
||
'aspect': {
|
||
'title': 'Aspect (slope direction)',
|
||
'legend': 'Direction in which the ground slopes down\nContinuous cycle: North→East→South→West→North\nPerceptually uniform colors (no hue jump)',
|
||
'description': 'Slope orientation — helps tell structures from natural landforms',
|
||
'cmap': 'twilight',
|
||
'gradient': ('#e2d9e2', '#95b5c7', '#6276ba', '#592a8f', '#2f1436',
|
||
'#741e4f', '#b25652', '#cca389', '#e2d9e2'),
|
||
'ticks': ('North 0°', 'North 360°'),
|
||
},
|
||
'mslrm': {
|
||
'title': 'MSRM - Multi-Scale Relief Model (adaptive scales)',
|
||
'legend': 'Combined multi-scale relief (local σ, fixed ±3σ scale)\nRed = raised (wall, mound, embankment)\nBlue = depression (ditch, moat)\n\nScales: 2 to 200 m, weighted towards 5–20 m\nConsistent colors across tiles\nDetects from micro to macro',
|
||
'description': 'Combines LRM at 5–7 scales — detects structures from 5 m to 100 m at once',
|
||
'cmap': 'seismic',
|
||
'gradient': ('#00004c', '#0000a6', '#0101ff', '#8181ff', '#fffdfd',
|
||
'#ff7d7d', '#fe0000', '#be0000', '#800000'),
|
||
'ticks': ('-3σ', '+3σ'),
|
||
},
|
||
'sailore': {
|
||
'title': 'SAILORE - Self-Adaptive LRM',
|
||
'legend': 'Adaptive local relief (local σ, fixed ±3σ scale)\nRed = raised | Blue = depression\nConsistent colors across tiles\n\nKernel adapted to the local slope\nFlat = large kernel (25 m) | Slope = small kernel (2 m)',
|
||
'description': 'Kernel that adapts to the local slope — flat ground = large kernel, slope = small kernel',
|
||
'cmap': 'seismic',
|
||
'gradient': ('#00004c', '#0000a6', '#0101ff', '#8181ff', '#fffdfd',
|
||
'#ff7d7d', '#fe0000', '#be0000', '#800000'),
|
||
'ticks': ('-3σ', '+3σ'),
|
||
},
|
||
'positive_openness': {
|
||
'title': 'Positive Openness (upward openness)',
|
||
'legend': 'Opening angle towards the sky (deviation from a fixed national reference)\nLight = open view of the sky (summits, plateaus)\nDark = blocked view (deep valleys)\nSame angle = same color on every tile',
|
||
'description': 'Ray tracing in 8 directions, multi-radius — detects ridges and summits',
|
||
'cmap': 'YlOrBr',
|
||
'gradient': ('#ffffe5', '#fff7bc', '#fee390', '#fec34f', '#fe9829',
|
||
'#eb6f14', '#cb4b02', '#983404', '#662506'),
|
||
'ticks': ('-3σ', '+3σ'),
|
||
},
|
||
'negative_openness': {
|
||
'title': 'Negative Openness (downward openness)',
|
||
'legend': 'Opening angle downwards (deviation from a fixed national reference)\nLight = overhang (ditch edges, caves)\nDark = flat ground (valley floors)\nSame angle = same color on every tile\nBest detector of cavities and sinkholes',
|
||
'description': 'Ray tracing in 8 directions, multi-radius — detects ditches, sinkholes, underground features',
|
||
'cmap': 'PuBu',
|
||
'gradient': ('#fff7fb', '#ece7f2', '#d0d1e6', '#a5bddb', '#73a9cf',
|
||
'#358fc0', '#056faf', '#04598c', '#023858'),
|
||
'ticks': ('-3σ', '+3σ'),
|
||
},
|
||
'svf': {
|
||
'title': 'Sky-View Factor (visible sky fraction)',
|
||
'legend': 'Share of visible sky (fixed physical 0–1 scale)\nWhite/yellow = sky hidden (valley, ditch, trench)\nBlack = open sky (summit, plateau)\nConsistent colors across tiles\nDitches stand out brightly — excellent for linear features',
|
||
'description': 'Micro-relief detection — ditches in yellow/white, banks in dark',
|
||
'cmap': 'hot_r',
|
||
'gradient': ('#ffffff', '#ffff81', '#ffff03', '#ffad00', '#ff5900',
|
||
'#ff0500', '#b00000', '#5c0000', '#0b0000'),
|
||
'ticks': ('0', '1'),
|
||
},
|
||
'roughness': {
|
||
'title': 'Multi-Scale Roughness (3 m + 15 m)',
|
||
'legend': 'Combined fine + broad terrain irregularity\nDark purple = smooth surface (road, wall, flat ground)\nBright yellow = rough surface (vegetation, ruins, stones)\nCombines fine 3 m roughness (70%) + broad 15 m (30%)\nPhysical scale shared by all tiles (fixed references\nmeasured on real tiles): seamless mosaic,\nsame value = same color',
|
||
'description': 'Measures local variability — smooth man-made surfaces vs rough natural ones',
|
||
'cmap': 'plasma',
|
||
'gradient': ('#0d0887', '#4c02a1', '#7e03a8', '#aa2395', '#cc4778',
|
||
'#e66c5c', '#f89540', '#fdc527', '#f0f921'),
|
||
'ticks': ('smooth', 'rough'),
|
||
},
|
||
'wavelet': {
|
||
'title': 'Mexican Hat Wavelet (multi-scale CWT)',
|
||
'legend': 'Multi-scale RMS index, centered on the tile\nmedian (1 = average level, higher = structure)\n\nLarge volumes removed (35 m local mean):\na ditch on a summit or a slope stands out\nno more than a ditch on flat ground\n\nFixed global quantile stretch (calibrated on a\nsample of tiles): same value = same color\non every tile and resolution\n\nTuned for small structures:\npaths, ditches, ramparts',
|
||
'description': '2D wavelet transform: detection of small structures (paths, ditches, ramparts)',
|
||
'cmap': 'inferno',
|
||
'gradient': ('#000004', '#210c4a', '#57106e', '#8a226a', '#bc3754',
|
||
'#e45a31', '#f98e09', '#f9cb35', '#fcffa4'),
|
||
'ticks': ('noise', 'structure'),
|
||
},
|
||
'flow_acc': {
|
||
'title': 'Flow Accumulation',
|
||
'legend': 'Log10 of the number of upstream cells\nDark green = high accumulation (ditch, channel, drainage)\nYellow = low accumulation (flat ground)\n\nDetects ditches and linear hydrological features',
|
||
'description': 'Priority-flood + D8 — detects archaeological ditches and drainage',
|
||
'cmap': 'YlGn',
|
||
'gradient': ('#ffffe5', '#f7fcb9', '#d9f0a3', '#acdd8e', '#77c679',
|
||
'#40aa5c', '#228343', '#006737', '#004529'),
|
||
'ticks': ('low', 'high'),
|
||
},
|
||
'solar': {
|
||
'title': 'Solar Illumination',
|
||
'legend': "Solar illumination (azimuth 90°, altitude 30°)\nLight = lit face | Dark = shadow",
|
||
'description': "Simulated morning sunlight",
|
||
'cmap': 'gray',
|
||
'gradient': ('#000000', '#202020', '#404040', '#606060', '#808080',
|
||
'#a0a0a0', '#c0c0c0', '#e0e0e0', '#ffffff'),
|
||
'ticks': ('0', '1'),
|
||
},
|
||
'anomaly': {
|
||
'title': 'Anomaly Map (automatic detection)',
|
||
'legend': 'Composite anomaly score (0–1)\nRed = strong anomaly (suspected structures)\nYellow = moderate anomaly\nWhite = no signal (natural ground)\n\nAuto threshold: pixels > 2σ from the local mean\nCombined: MSRM, SVF, wavelet, openness, roughness',
|
||
'description': 'Automatic detection — targets to be checked in the field',
|
||
'cmap': 'YlOrRd',
|
||
'gradient': ('#ffffcc', '#ffeda0', '#fed976', '#feb24c', '#fd8c3c',
|
||
'#fc4d2a', '#e2191c', '#bb0026', '#800026'),
|
||
'ticks': ('0', '1'),
|
||
},
|
||
'relief_oriente': {
|
||
'title': 'Oriented relief (local openness × orientation)',
|
||
'legend': 'Lightness = micro-relief (local openness 5–20 m + shading)\nLight = bump, ridge | Dark = hollow, ditch\nHue = slope orientation\nFixed scale — consistent colors across tiles',
|
||
'description': 'Openness on the detrended DTM (σ 10 m), radii 5/10/20 m, 16 directions; CIELAB hue = aspect',
|
||
'cmap': None,
|
||
'gradient': None,
|
||
'ticks': None,
|
||
# "How to read": single text source for the map (/api/map/meta) and
|
||
# the PDF sheet legend.
|
||
'reading': (
|
||
'Lightness = local openness of the terrain: light = bump, ridge; dark = hollow, ditch.',
|
||
'Hue = slope orientation (see the rose).',
|
||
'Black gaps = no ground point: buildings, water, dense cover.',
|
||
"Between points, the relief is filled in only inside the envelope "
|
||
"of the points, over a radius of 1.5 × the local spacing (at least "
|
||
"1 m): enough to avoid holes, without inventing relief or "
|
||
"amplifying noise.",
|
||
),
|
||
},
|
||
'densite_sol': {
|
||
'title': 'Geometric precision (ground point density)',
|
||
'legend': 'Ground points kept for the DTM, per m² (3 × 3 m mean)\n16 grays, fixed log scale: 2 levels = density doubled\nBlack = ≤ 0.35 pt/m² or no point (relief interpolated or missing)\nWhite = ≥ 45 pts/m²',
|
||
'description': 'Where the relief is measured (light) and where it is interpolated (dark)',
|
||
'cmap': 'gray',
|
||
'gradient': ('#000000', '#202020', '#404040', '#606060', '#808080',
|
||
'#a0a0a0', '#c0c0c0', '#e0e0e0', '#ffffff'),
|
||
'ticks': ('≤ 0.35 pt/m²', '≥ 45 pts/m²'),
|
||
'reading': (
|
||
'Ground points per m², 16 grays: light = measured relief, dark = interpolated relief.',
|
||
'Black = less than 0.35 point per m², or no point at all.',
|
||
),
|
||
},
|
||
'ortho': {
|
||
'title': 'IGN Aerial Photograph',
|
||
'legend': 'Orthophoto\nAerial image',
|
||
'description': 'IGN aerial photograph (orthophoto)',
|
||
'cmap': None,
|
||
'gradient': None,
|
||
'ticks': None,
|
||
},
|
||
'topo': {
|
||
'title': 'IGN Topographic Map',
|
||
'legend': 'IGN map\nTopographic map',
|
||
'description': 'IGN topographic map (Plan IGN)',
|
||
'cmap': None,
|
||
'gradient': None,
|
||
'ticks': None,
|
||
},
|
||
}
|
||
|
||
# Map display: ONE main layer (DEFAULT_VIZ) and the "precision" layer
|
||
# (ground point density), shown alone or compared with the relief on either
|
||
# side of a sliding bar ("compare" mode).
|
||
PRECISION_VIZ = 'densite_sol'
|
||
VIEW_MODES = ('relief', 'precision', 'compare')
|
||
DEFAULT_VIEW_MODE = 'relief'
|
||
|
||
# "Main" layer (shown by default on the map): the oriented relief, which
|
||
# merges openness and aspect.
|
||
DEFAULT_VIZ = 'relief_oriente'
|
||
|
||
# Layers produced and displayed: only this selection is generated by
|
||
# default (pipeline without --only, generation from the map) and served by
|
||
# the map (panel, XYZ tiles, TileJSON, WMTS, JOSM). The other visualizations
|
||
# can still be computed with --only but are no longer offered.
|
||
# None = every visualization present on disk.
|
||
PANEL_VIZ = ('relief_oriente', 'densite_sol')
|
||
|
||
# Output file keyword → pipeline --only step name (the three visualizations
|
||
# whose output name differs from the step name, see _expected_output_path in
|
||
# pipeline.py).
|
||
KEYWORD_TO_STEP = {
|
||
'hillshade_multi': 'hillshade',
|
||
'positive_openness': 'pos_open',
|
||
'negative_openness': 'neg_open',
|
||
}
|
||
|
||
# Reverse mapping: --only step name → output file keyword.
|
||
STEP_TO_KEYWORD = {step: kw for kw, step in KEYWORD_TO_STEP.items()}
|
||
|
||
|
||
def panel_steps():
|
||
"""--only step names of the layers produced by default (PANEL_VIZ)."""
|
||
if PANEL_VIZ is None:
|
||
return None
|
||
return [KEYWORD_TO_STEP.get(k, k) for k in PANEL_VIZ]
|
||
|
||
|
||
def step_to_keyword(step):
|
||
"""Pipeline step name (e.g. 'pos_open') → file keyword ('positive_openness')."""
|
||
return STEP_TO_KEYWORD.get(step, step)
|
||
|
||
|
||
def default_main_layer(all_viz_keys):
|
||
"""Default main layer present on disk (DEFAULT_VIZ, otherwise the first
|
||
one that is not the precision layer), or None."""
|
||
keys = [k for k in all_viz_keys if k != PRECISION_VIZ]
|
||
if DEFAULT_VIZ in keys:
|
||
return DEFAULT_VIZ
|
||
return keys[0] if keys else None
|
||
|
||
|
||
def cells_with_all_viz(vis_dir, viz_keys, resolutions=(0.5,)):
|
||
"""Cells (col, row) that have ALL the requested visualizations.
|
||
|
||
A cell is complete if, for every resolution in `resolutions`, a
|
||
visualization directory matches it and contains every keyword of
|
||
`viz_keys` (e.g. 'aspect', 'hillshade_multi'). Used to tell tiles that
|
||
are really finished from those still to be completed: an existing but
|
||
incomplete tile (missing visualization or resolution) is still to process.
|
||
|
||
Returns:
|
||
Set of complete (col, row).
|
||
"""
|
||
by_cell = {}
|
||
for t in scan_tiles(vis_dir):
|
||
per_res = by_cell.setdefault((t['col'], t['row']), {})
|
||
per_res.setdefault(t['resolution'], set()).update(t['viz'].keys())
|
||
return {cell for cell, per_res in by_cell.items()
|
||
if all(kw in per_res.get(res, ()) for res in resolutions for kw in viz_keys)}
|
||
|
||
# Every visualization is cut into sub-tiles (500 m quadrants) to lighten the
|
||
# map. Set a tuple to restrict the cutting — excluded visualizations fall
|
||
# back to the whole tile.
|
||
_CARTO_SUBTILED_VIZ = ()
|
||
|
||
# Sub-tile AVIF encoding: q75 in 4:2:0, encoded ONCE from the original
|
||
# raster (write_subtiles called by tif_to_crop). Measured on the oriented
|
||
# relief (2 real tiles, compression gallery): the former chain tile q60 →
|
||
# sub-tile q55 gave 18.1 dB / SSIM 0.84 for 3.9 MB per tile; a single q75
|
||
# gives 19.5 dB / SSIM 0.93 for 8.8 MB. Beyond that, 4:2:0 hits a ceiling
|
||
# (the relief hue, pixel by pixel, is averaged over 2 × 2): only 4:4:4 would
|
||
# go further (q75: 28 dB, 14 MB).
|
||
# speed 9: fast encoding (see rendering.AVIF_SPEED).
|
||
_SUBTILE_AVIF_QUALITY = 75
|
||
_SUBTILE_AVIF_SPEED = 9
|
||
|
||
# Sub-tile thumbnail (px): a small source level of the XYZ pyramid — 160 px
|
||
# covers display up to ~220 px on screen and cuts decoded memory by a factor
|
||
# of ~2.5 vs 256 px. The size is encoded in the file name: changing it
|
||
# invalidates the cache.
|
||
_SUBTILE_THUMB_PX = 160
|
||
|
||
# Layers with photographic content or thin lines (orthophoto, topo map):
|
||
# their own quality (currently equal to that of the color ramps).
|
||
_SUBTILE_AVIF_QUALITY_DETAIL = 75
|
||
_SUBTILE_DETAIL_VIZ = frozenset({'ortho', 'topo'})
|
||
|
||
# Layers made of flat coded levels (see rendering.LOSSLESS_GRAY_KEYWORDS):
|
||
# sub-tiles in lossless WebP (grayscale; 3× lighter than AVIF q100, the only
|
||
# exact AVIF setting) and a lossless intermediate thumbnail.
|
||
_SUBTILE_LOSSLESS_VIZ = frozenset({'densite_sol'})
|
||
|
||
# Intermediate thumbnail (px): a level between the 256 px thumbnail and the
|
||
# full-resolution image, so the pyramid neither stretches the thumbnail nor
|
||
# decodes the full AVIF as soon as a tile exceeds ~300 px on screen.
|
||
_MID_THUMB_SIZE = 640
|
||
|
||
# Preferred layer order (inventory viz_meta order) and choice of each tile's
|
||
# fallback display layer (_pick_display_viz).
|
||
_VIZ_FALLBACK_ORDER = [
|
||
'hillshade_multi', 'svf', 'slope', 'mslrm', 'positive_openness',
|
||
'negative_openness', 'aspect', 'sailore', 'roughness',
|
||
'wavelet', 'flow_acc', 'solar', 'anomaly', 'relief_oriente', 'ortho', 'topo',
|
||
]
|
||
|
||
# Regex parsing the tile coordinates in the LHD basename.
|
||
# LHD_FXX_{COL}_{ROW}_PTS_LAMB93_IGN69 (COL/ROW in km, Lambert 93)
|
||
_RE_LHD_COORDS = re.compile(r'^LHD_FXX_(\d+)_(\d+)_PTS_LAMB93')
|
||
|
||
|
||
def parse_basename_coords(name):
|
||
"""Extract the tile coordinates (col, row in km) from a basename.
|
||
|
||
Args:
|
||
name: candidate basename (e.g. 'LHD_FXX_1000_6881_PTS_LAMB93_IGN69')
|
||
or directory name with a resolution suffix ('..._r0p2').
|
||
|
||
Returns:
|
||
(col_km, row_km), or None if the name does not match the LHD pattern.
|
||
"""
|
||
m = _RE_LHD_COORDS.match(name)
|
||
if not m:
|
||
return None
|
||
return int(m.group(1)), int(m.group(2))
|
||
|
||
|
||
def _strip_res_suffix(dirname):
|
||
"""Split a visualization directory name into base basename and resolution.
|
||
|
||
'LHD_FXX_1000_6881_PTS_LAMB93_IGN69' → (basename, 0.5)
|
||
'LHD_FXX_1000_6881_PTS_LAMB93_IGN69_r0p2' → (basename, 0.2)
|
||
|
||
Returns:
|
||
(basename_without_suffix, resolution_float), or (dirname, 0.5) without a suffix.
|
||
"""
|
||
m = re.match(r'^(.+?)_r(\d+p\d+)$', dirname)
|
||
if m:
|
||
res_str = m.group(2).replace('p', '.')
|
||
try:
|
||
return m.group(1), float(res_str)
|
||
except ValueError:
|
||
pass
|
||
return dirname, 0.5
|
||
|
||
|
||
def _res_suffix_str(resolution):
|
||
"""Naming suffix of a resolution (0.5 m = primary resolution, no suffix).
|
||
|
||
Same convention as pipeline.LidarArchaeoPipeline._res_suffix, but
|
||
reimplemented locally: importing the pipeline would pull in dtm→numpy,
|
||
absent from the lightweight lidar-maps image (Dockerfile.maps), and would
|
||
break build_index there. Any change to the format must stay in sync with
|
||
pipeline.py (_res_suffix) and the reverse decoding (_strip_res_suffix).
|
||
"""
|
||
if resolution == 0.5:
|
||
return ""
|
||
return f"_r{f'{resolution}'.replace('.', 'p')}"
|
||
|
||
|
||
def scan_tiles(vis_dir):
|
||
"""Scan the visualization directory to inventory the processed tiles.
|
||
|
||
Args:
|
||
vis_dir: Path to output/visualisations/
|
||
|
||
Returns:
|
||
List of dicts:
|
||
{basename, col, row, resolution, dir_path,
|
||
viz: {viz_key: {filename, ext}}, dir_name}
|
||
Sorted by (resolution, decreasing row, col).
|
||
"""
|
||
vis_dir = Path(vis_dir)
|
||
if not vis_dir.is_dir():
|
||
return []
|
||
|
||
tiles = []
|
||
for entry in sorted(vis_dir.iterdir()):
|
||
if not entry.is_dir():
|
||
continue
|
||
coords = parse_basename_coords(entry.name)
|
||
if coords is None:
|
||
continue
|
||
col, row = coords
|
||
basename, resolution = _strip_res_suffix(entry.name)
|
||
|
||
# List the visualization image files in the directory.
|
||
viz = {}
|
||
for f in sorted(entry.iterdir()):
|
||
if not f.is_file():
|
||
continue
|
||
# Detect the AVIF/WebP extension
|
||
ext = None
|
||
low = f.name.lower()
|
||
for e in ('.avif', '.webp'):
|
||
if low.endswith(e):
|
||
ext = e.lstrip('.')
|
||
break
|
||
if ext is None:
|
||
continue
|
||
# viz_key = name without the basename_ prefix and the extension
|
||
stem = f.name[:-len('.' + ext)]
|
||
prefix = basename + '_'
|
||
if not stem.startswith(prefix):
|
||
continue
|
||
viz_key = stem[len(prefix):]
|
||
viz[viz_key] = {'filename': f.name, 'ext': ext}
|
||
|
||
if not viz:
|
||
# Empty directory or no valid image → skipped
|
||
continue
|
||
|
||
tiles.append({
|
||
'basename': basename,
|
||
'col': col,
|
||
'row': row,
|
||
'resolution': resolution,
|
||
'dir_name': entry.name,
|
||
'dir_path': str(entry),
|
||
'viz': viz,
|
||
})
|
||
|
||
tiles.sort(key=lambda t: (t['resolution'], -t['row'], t['col']))
|
||
return tiles
|
||
|
||
|
||
def compute_bbox(tiles):
|
||
"""Compute the bounding box (in km) covered by the tiles.
|
||
|
||
Returns:
|
||
Dict {min_col, max_col, min_row, max_row}, or None if there is no tile.
|
||
"""
|
||
if not tiles:
|
||
return None
|
||
cols = [t['col'] for t in tiles]
|
||
rows = [t['row'] for t in tiles]
|
||
return {
|
||
'min_col': min(cols),
|
||
'max_col': max(cols),
|
||
'min_row': min(rows),
|
||
'max_row': max(rows),
|
||
}
|
||
|
||
|
||
def compute_zones(tiles, proximity_threshold=15):
|
||
"""Group the tiles into geographic zones by proximity clustering.
|
||
|
||
Tiles within proximity_threshold km of each other are grouped in the
|
||
same zone. Zones are sorted by decreasing size (largest zone first).
|
||
|
||
Args:
|
||
tiles: List of tile dicts.
|
||
proximity_threshold: Maximum distance in km to group two tiles.
|
||
|
||
Returns:
|
||
List of zone dicts:
|
||
{label, tiles, bbox}
|
||
"""
|
||
if not tiles:
|
||
return []
|
||
|
||
# Union-Find for the clustering
|
||
parent = list(range(len(tiles)))
|
||
|
||
def find(x):
|
||
while parent[x] != x:
|
||
parent[x] = parent[parent[x]]
|
||
x = parent[x]
|
||
return x
|
||
|
||
def union(x, y):
|
||
px, py = find(x), find(y)
|
||
if px != py:
|
||
parent[px] = py
|
||
|
||
# Group nearby tiles
|
||
for i in range(len(tiles)):
|
||
for j in range(i + 1, len(tiles)):
|
||
dc = abs(tiles[i]['col'] - tiles[j]['col'])
|
||
dr = abs(tiles[i]['row'] - tiles[j]['row'])
|
||
if max(dc, dr) <= proximity_threshold:
|
||
union(i, j)
|
||
|
||
# Build the zones
|
||
zone_members = {}
|
||
for i in range(len(tiles)):
|
||
root = find(i)
|
||
if root not in zone_members:
|
||
zone_members[root] = []
|
||
zone_members[root].append(tiles[i])
|
||
|
||
zones = [{'tiles': zone_tiles, 'bbox': compute_bbox(zone_tiles)}
|
||
for zone_tiles in zone_members.values()]
|
||
# Sort BEFORE numbering: the "Zone N" labels follow the display order
|
||
# (decreasing size)
|
||
zones.sort(key=lambda z: len(z['tiles']), reverse=True)
|
||
for i, z in enumerate(zones, 1):
|
||
z['label'] = f'Zone {i} ({len(z["tiles"])} tiles)'
|
||
return zones
|
||
|
||
|
||
def _approx_l93_to_wgs84(x_m, y_m):
|
||
"""Affine approximation Lambert 93 → WGS84 (fallback without rasterio/pyproj).
|
||
|
||
Exact origin: (700000, 6600000) L93 ↔ (3.0°E, 46.5°N).
|
||
Accuracy of the order of a km — only used when neither rasterio nor
|
||
pyproj is available (never the case in the Docker images).
|
||
"""
|
||
import math
|
||
lat = 46.5 + (y_m - 6600000.0) / 111320.0
|
||
lon = 3.0 + (x_m - 700000.0) / (111320.0 * math.cos(math.radians(47.0)))
|
||
return lon, lat
|
||
|
||
|
||
def _approx_wgs84_to_l93(lon, lat):
|
||
"""Affine approximation WGS84 → Lambert 93 (exact inverse of the previous one).
|
||
|
||
Accuracy of the order of a km — fallback without pyproj (see
|
||
bbox_to_cells / point_to_cell in mapserve.py).
|
||
"""
|
||
import math
|
||
y = (lat - 46.5) * 111320.0 + 6600000.0
|
||
x = (lon - 3.0) * (111320.0 * math.cos(math.radians(47.0))) + 700000.0
|
||
return x, y
|
||
|
||
|
||
def attach_gps_bounds(tiles):
|
||
"""Attach to each tile its GPS corners for the Leaflet display.
|
||
|
||
Each 1×1 km tile is defined by its north-west corner in L93 km
|
||
(col, row) → X ∈ [col, col+1] km, Y ∈ [row-1, row] km.
|
||
(Checked against the DTM bounds: X_min = col×1000, Y_max = row×1000.)
|
||
|
||
Uses rasterio.warp (exact PROJ conversion) if available, otherwise
|
||
pyproj (lightweight image without GDAL), otherwise the affine
|
||
approximation _approx_l93_to_wgs84 (accuracy ~km).
|
||
|
||
Adds to each tile:
|
||
corners: [[lat, lon] × 4] in the order SW, SE, NE, NW
|
||
bounds : [[lat_south, lon_west], [lat_north, lon_east]]
|
||
"""
|
||
# Corners SW, SE, NE, NW — south edge Y = (row-1)×1000, north edge = row×1000
|
||
xs = []
|
||
ys = []
|
||
for t in tiles:
|
||
xs.extend([t['col'] * 1000, (t['col'] + 1) * 1000,
|
||
(t['col'] + 1) * 1000, t['col'] * 1000])
|
||
ys.extend([(t['row'] - 1) * 1000, (t['row'] - 1) * 1000,
|
||
t['row'] * 1000, t['row'] * 1000])
|
||
try:
|
||
from rasterio.warp import transform as warp_transform
|
||
lons, lats = warp_transform('EPSG:2154', 'EPSG:4326', xs, ys)
|
||
ok = True
|
||
except ImportError:
|
||
try:
|
||
from pyproj import Transformer
|
||
transformer = Transformer.from_crs('EPSG:2154', 'EPSG:4326',
|
||
always_xy=True)
|
||
lons, lats = transformer.transform(xs, ys)
|
||
ok = True
|
||
except ImportError:
|
||
logger.debug("Approximate GPS corners (rasterio and pyproj unavailable)")
|
||
lons = None
|
||
lats = None
|
||
ok = False
|
||
except Exception as e:
|
||
logger.debug(f"Approximate GPS corners (rasterio unavailable: {e})")
|
||
lons = None
|
||
lats = None
|
||
ok = False
|
||
|
||
for i, t in enumerate(tiles):
|
||
if ok:
|
||
corners = [[lats[4 * i], lons[4 * i]],
|
||
[lats[4 * i + 1], lons[4 * i + 1]],
|
||
[lats[4 * i + 2], lons[4 * i + 2]],
|
||
[lats[4 * i + 3], lons[4 * i + 3]]]
|
||
else:
|
||
corners = []
|
||
for cx in (t['col'] * 1000, (t['col'] + 1) * 1000):
|
||
for cy in ((t['row'] - 1) * 1000, t['row'] * 1000):
|
||
lon, lat = _approx_l93_to_wgs84(cx, cy)
|
||
corners.append([lat, lon])
|
||
# Reorder to SW, SE, NE, NW (the loop yields SW, NW, SE, NE)
|
||
corners = [corners[0], corners[2], corners[3], corners[1]]
|
||
t['corners'] = corners
|
||
t['bounds'] = [[min(c[0] for c in corners), min(c[1] for c in corners)],
|
||
[max(c[0] for c in corners), max(c[1] for c in corners)]]
|
||
return ok
|
||
|
||
|
||
def _mtime(path):
|
||
"""Mtime of a file, or None if inaccessible."""
|
||
try:
|
||
return Path(path).stat().st_mtime
|
||
except OSError:
|
||
return None
|
||
|
||
|
||
def _cached_file_fresh(path, src_mtime):
|
||
"""True if a cached file exists and is newer than its source.
|
||
|
||
Used to invalidate thumbnails and sub-tiles when a tile is recomputed:
|
||
the source image (AVIF/WebP) being rewritten, its mtime becomes newer
|
||
than the cache's, which must then be regenerated.
|
||
"""
|
||
cached_mtime = _mtime(path)
|
||
if cached_mtime is None:
|
||
return False
|
||
return src_mtime is None or cached_mtime >= src_mtime
|
||
|
||
|
||
def _url_version(mtime):
|
||
"""Cache-busting suffix for an image URL, or '' if unknown.
|
||
|
||
The map images are served with an immutable cache when the URL carries
|
||
?v=: the suffix MUST therefore identify the content of the served file,
|
||
not that of its source (a thumbnail recomputed later, or a file fetched
|
||
from the upstream, changes content without its source moving). Each URL
|
||
is versioned by the mtime of ITS file. A recomputation changes the URL
|
||
and forces a reload — including live during a run, when the inventory is
|
||
rewritten after each tile.
|
||
"""
|
||
return f"?v={int(mtime * 1000)}" if mtime is not None else ""
|
||
|
||
|
||
def generate_thumbnail(src_path, thumb_path, max_size=256, mid_path=None,
|
||
mid_size=640):
|
||
"""Generate a JPEG thumbnail from an existing AVIF/WebP image.
|
||
|
||
Args:
|
||
src_path: path of the source image (AVIF/WebP).
|
||
thumb_path: JPEG output path.
|
||
max_size: maximum size (longest side) in pixels.
|
||
mid_path: optional intermediate thumbnail (JPEG, size mid_size).
|
||
mid_size: maximum size of the intermediate thumbnail.
|
||
|
||
Returns:
|
||
True if the main thumbnail is OK, False on failure.
|
||
"""
|
||
try:
|
||
from PIL import Image as PILImage
|
||
except ImportError:
|
||
logger.warning("PIL unavailable — cannot generate thumbnails")
|
||
return False
|
||
|
||
try:
|
||
try:
|
||
resample = PILImage.Resampling.LANCZOS
|
||
except AttributeError:
|
||
resample = getattr(PILImage, 'LANCZOS', 1)
|
||
|
||
def resized(source, target):
|
||
scale = min(1.0, target / max(source.size))
|
||
if scale >= 1.0:
|
||
return source
|
||
return source.resize((max(1, int(source.size[0] * scale)),
|
||
max(1, int(source.size[1] * scale))), resample)
|
||
|
||
with PILImage.open(str(src_path)) as _src_img:
|
||
img = _src_img.convert('RGB')
|
||
Path(thumb_path).parent.mkdir(parents=True, exist_ok=True)
|
||
if mid_path is not None:
|
||
try:
|
||
resized(img, mid_size).save(str(mid_path), format='JPEG', quality=82)
|
||
except Exception as e:
|
||
logger.debug(f"Intermediate thumbnail skipped {src_path}: {e}")
|
||
resized(img, max_size).save(str(thumb_path), format='JPEG', quality=80)
|
||
return True
|
||
except Exception as e:
|
||
logger.debug(f"Thumbnail skipped {src_path}: {e}")
|
||
return False
|
||
|
||
|
||
def _pick_display_viz(viz_keys):
|
||
"""Choose the default visualization of a tile.
|
||
|
||
Prefers hillshade_multi, otherwise the first available in the fallback order.
|
||
"""
|
||
for v in _VIZ_FALLBACK_ORDER:
|
||
if v in viz_keys:
|
||
return v
|
||
return sorted(viz_keys)[0]
|
||
|
||
|
||
def _subdivision_k(resolution, tile_m=1000, target_px=2500):
|
||
"""Split factor k (k×k grid) to lighten map rendering.
|
||
|
||
At 0.2 m/px a 1 km tile is 5000×5000 px (~100 MB decoded): it is cut
|
||
into 500 m quadrants (k=2, 2500×2500 px). At 0.5 m/px (2000 px) the tile
|
||
stays whole (k=1).
|
||
"""
|
||
px = max(1, int(round(tile_m / resolution)))
|
||
return max(1, int(round(px / target_px)))
|
||
|
||
|
||
def _subtile_corners(corners, i, j, k):
|
||
"""WGS84 corners [SW, SE, NE, NW] of sub-tile (i, j) of a k×k split.
|
||
|
||
i: index towards the east (0..k-1), j: index towards the north (0..k-1).
|
||
Bilinear interpolation of the tile corners — the projected quadrilateral
|
||
is almost a parallelogram at this scale (screen error < 1 px).
|
||
|
||
Edges are shared between neighboring sub-tiles: every grid point is
|
||
computed from the same integer grid indices, so two adjacent sub-tiles
|
||
get exactly the same point (evaluating the interpolation on slightly
|
||
different fractions would make edges miss by less than a pixel — a
|
||
broken seam at medium zoom).
|
||
"""
|
||
sw, se, ne, nw = corners
|
||
|
||
def lerp(p, q, u):
|
||
return [p[0] + (q[0] - p[0]) * u, p[1] + (q[1] - p[1]) * u]
|
||
|
||
# Shared edges: each point of the (k+1)×(k+1) grid is derived from its
|
||
# integer indices (i0, j0) only, so two neighboring sub-tiles get exactly
|
||
# the same point.
|
||
def at(i0, j0):
|
||
u = i0 / k
|
||
v = j0 / k
|
||
bottom = lerp(sw, se, u) # along the south edge, position u
|
||
top = lerp(nw, ne, u) # along the north edge, position u
|
||
return lerp(bottom, top, v) # north-south interpolation
|
||
|
||
# Corners SW, SE, NE, NW of sub-tile (i, j).
|
||
return [at(i, j), at(i + 1, j), at(i + 1, j + 1), at(i, j + 1)]
|
||
|
||
|
||
def _fallback_full_dalle(entries, viz_key, info):
|
||
"""Whole-tile fallback for a layer that cannot be cut: it still works as
|
||
a layer (heavier images, but functional)."""
|
||
for entry in entries.values():
|
||
entry['viz'][viz_key] = dict(info)
|
||
|
||
|
||
def _subtile_ext(viz_key):
|
||
"""Extension of a layer's full-resolution sub-tiles."""
|
||
return '.webp' if viz_key in _SUBTILE_LOSSLESS_VIZ else '.avif'
|
||
|
||
|
||
def _save_subtile_thumb(quad, path):
|
||
"""Write the sub-tile thumbnail (resized to _SUBTILE_THUMB_PX)."""
|
||
from PIL import Image as PILImage
|
||
scale = min(1.0, _SUBTILE_THUMB_PX / max(quad.size))
|
||
out = quad
|
||
if scale < 1.0:
|
||
out = quad.resize((max(1, int(quad.size[0] * scale)),
|
||
max(1, int(quad.size[1] * scale))), PILImage.LANCZOS)
|
||
out.save(str(path), format='WEBP', quality=80)
|
||
|
||
|
||
def write_subtiles(output_dir, dir_name, viz_key, img, k, sub_dir_name='index_subtiles'):
|
||
"""Cut a tile image (PIL, north up) into k × k sub-tiles: full
|
||
resolution, intermediate thumbnail and thumbnail.
|
||
|
||
Called by build_index from the tile image, and by the pipeline
|
||
(rendering.tif_to_crop) directly from the original raster: the sub-tile
|
||
then undergoes a single lossy encoding (re-encoding the AVIF tile
|
||
compounded two losses). Encoding: AVIF 4:2:0 q75 (_SUBTILE_AVIF_QUALITY,
|
||
q75 for ortho/topo too), lossless WebP for flat level layers
|
||
(_SUBTILE_LOSSLESS_VIZ). Raises on failure.
|
||
"""
|
||
from PIL import Image as PILImage
|
||
output_dir = Path(output_dir)
|
||
out_dir = output_dir / sub_dir_name
|
||
out_dir.mkdir(parents=True, exist_ok=True)
|
||
lossless = viz_key in _SUBTILE_LOSSLESS_VIZ
|
||
ext = _subtile_ext(viz_key)
|
||
other_ext = '.avif' if ext == '.webp' else '.webp'
|
||
quality = (_SUBTILE_AVIF_QUALITY_DETAIL if viz_key in _SUBTILE_DETAIL_VIZ
|
||
else _SUBTILE_AVIF_QUALITY)
|
||
if img.mode not in ('RGB', 'L'):
|
||
img = img.convert('RGB')
|
||
if lossless:
|
||
img = img.convert('L')
|
||
W, H = img.size
|
||
for j in range(k):
|
||
for i in range(k):
|
||
stem = f"{dir_name}_{viz_key}_{i}_{j}"
|
||
# Image: row 0 = north → the northern sub-tile j is at the top
|
||
left, right = round(W * i / k), round(W * (i + 1) / k)
|
||
top = round(H * (1 - (j + 1) / k))
|
||
bottom = round(H * (1 - j / k))
|
||
quad = img.crop((left, top, right, bottom))
|
||
if lossless:
|
||
quad.save(str(out_dir / (stem + ext)), format='WEBP', lossless=True)
|
||
else:
|
||
quad.save(str(out_dir / (stem + ext)), format='AVIF', quality=quality,
|
||
subsampling='4:2:0', speed=_SUBTILE_AVIF_SPEED)
|
||
# Old files from a previous generation (other format, 256 px
|
||
# thumbnail without the size in the name)
|
||
(out_dir / (stem + other_ext)).unlink(missing_ok=True)
|
||
(out_dir / (stem + '_thumb.webp')).unlink(missing_ok=True)
|
||
mid_scale = min(1.0, _MID_THUMB_SIZE / max(quad.size))
|
||
mid_img = quad
|
||
if mid_scale < 1.0:
|
||
mid_img = quad.resize((max(1, int(quad.size[0] * mid_scale)),
|
||
max(1, int(quad.size[1] * mid_scale))), PILImage.LANCZOS)
|
||
mid_img.save(str(out_dir / (stem + '_mid.webp')), format='WEBP',
|
||
**({'lossless': True} if lossless else {'quality': 82}))
|
||
_save_subtile_thumb(quad, out_dir / (stem + f"_thumb{_SUBTILE_THUMB_PX}.webp"))
|
||
|
||
|
||
def _build_subtiles(tile, offered_viz_keys, output_dir, sub_dir_name):
|
||
"""Cut a tile into sub-tiles (AVIF crops) for the interactive map.
|
||
|
||
Only cuts the visualizations in offered_viz_keys. Returns the list of
|
||
display entries (one per sub-tile), or None if cutting is not
|
||
needed/possible (the whole tile is then displayed).
|
||
"""
|
||
k = _subdivision_k(tile['resolution'])
|
||
if k <= 1:
|
||
return None
|
||
try:
|
||
from PIL import Image as PILImage
|
||
except ImportError:
|
||
return None
|
||
|
||
sub_dir = output_dir / sub_dir_name
|
||
sub_dir.mkdir(parents=True, exist_ok=True)
|
||
|
||
entries = {}
|
||
for j in range(k):
|
||
for i in range(k):
|
||
corners = _subtile_corners(tile['corners'], i, j, k)
|
||
entries[(i, j)] = {
|
||
'col': tile['col'], 'row': tile['row'],
|
||
'name': tile['name'], 'dir_name': tile['dir_name'],
|
||
'resolution': tile['resolution'],
|
||
'bounds': [[min(c[0] for c in corners), min(c[1] for c in corners)],
|
||
[max(c[0] for c in corners), max(c[1] for c in corners)]],
|
||
'corners': corners,
|
||
'display_viz': tile['display_viz'],
|
||
'viz': {},
|
||
'meta': tile.get('meta'),
|
||
'sub_i': i, 'sub_j': j, 'sub_k': k,
|
||
'size_km': round(1.0 / k, 3),
|
||
}
|
||
|
||
try:
|
||
thumb_suffix = f"_thumb{_SUBTILE_THUMB_PX}.webp"
|
||
|
||
for viz_key in offered_viz_keys:
|
||
info = tile['viz'].get(viz_key)
|
||
if not info:
|
||
continue
|
||
ext = _subtile_ext(viz_key)
|
||
stems = {key: f"{tile['dir_name']}_{viz_key}_{key[0]}_{key[1]}"
|
||
for key in entries}
|
||
# Regenerate if at least one file is missing or stale (source
|
||
# tile recomputed since — like the thumbnails). The full URL
|
||
# carries a ?v= cache-busting suffix: strip it for the path.
|
||
# Sub-tiles written by the pipeline from the original raster
|
||
# (tif_to_crop) are newer than the tile: they are kept as is
|
||
# (a single lossy encoding).
|
||
src = output_dir / info['full'].split('?')[0]
|
||
src_mtime = _mtime(src)
|
||
|
||
def _fresh(name):
|
||
return _cached_file_fresh(output_dir / sub_dir_name / name, src_mtime)
|
||
|
||
# Full resolution + intermediate thumbnail on one side, thumbnails
|
||
# on the other: a thumbnail size change (source tile unchanged)
|
||
# re-cuts from the existing sub-tiles, without re-encoding them
|
||
# (minutes of CPU per rebuild).
|
||
heavy_fresh = all(_fresh(stem + ext) and _fresh(stem + '_mid.webp')
|
||
for stem in stems.values())
|
||
thumbs_fresh = all(_fresh(stem + thumb_suffix) for stem in stems.values())
|
||
if heavy_fresh and not thumbs_fresh:
|
||
logger.info(f" Sub-tile thumbnails recomputed: "
|
||
f"{tile['dir_name']}/{viz_key} ({len(stems)} crops)")
|
||
try:
|
||
for stem in stems.values():
|
||
with PILImage.open(str(output_dir / sub_dir_name / (stem + ext))) as q:
|
||
q.load()
|
||
if q.mode not in ('RGB', 'L'):
|
||
q = q.convert('RGB')
|
||
_save_subtile_thumb(q, output_dir / sub_dir_name / (stem + thumb_suffix))
|
||
# 256 px thumbnail from a previous generation
|
||
(output_dir / sub_dir_name / (stem + '_thumb.webp')).unlink(missing_ok=True)
|
||
except Exception as e:
|
||
logger.debug(f"Could not re-cut the thumbnails ({viz_key}): {e}")
|
||
heavy_fresh = False
|
||
if not heavy_fresh:
|
||
logger.info(f" Sub-tiles recomputed: {tile['dir_name']}/{viz_key} "
|
||
f"({len(stems)} crops)")
|
||
img = None
|
||
for attempt in range(2):
|
||
try:
|
||
img = PILImage.open(str(src))
|
||
img.load()
|
||
break
|
||
except Exception as e:
|
||
# Incremental rebuild during a run: the source tile may
|
||
# be being written by another worker (partial read).
|
||
# One retry after a short pause is almost always
|
||
# enough.
|
||
if attempt == 0:
|
||
time.sleep(2.0)
|
||
continue
|
||
logger.warning(f"Could not cut {viz_key} into sub-tiles after "
|
||
f"a retry ({src.name}): {e}")
|
||
_fallback_full_dalle(entries, viz_key, info)
|
||
if img is None:
|
||
continue
|
||
try:
|
||
write_subtiles(output_dir, tile['dir_name'], viz_key, img, k, sub_dir_name)
|
||
except Exception as e:
|
||
logger.warning(f"Could not encode the sub-tiles ({tile['dir_name']}), "
|
||
f"whole-tile fallback for {viz_key}: {e}")
|
||
_fallback_full_dalle(entries, viz_key, info)
|
||
continue
|
||
# Each URL is versioned by the mtime of ITS file (see
|
||
# _url_version): the "thumbnails only recomputed" shortcut above
|
||
# does not change the URLs of the AVIF/mid files not rewritten —
|
||
# the browser's immutable cache stays valid for them.
|
||
def _file_v(name):
|
||
return _url_version(_mtime(output_dir / sub_dir_name / name))
|
||
|
||
for (i, j), stem in stems.items():
|
||
entries[(i, j)]['viz'][viz_key] = {
|
||
'thumb': f"{sub_dir_name}/{stem}{thumb_suffix}{_file_v(stem + thumb_suffix)}",
|
||
'mid': f"{sub_dir_name}/{stem}_mid.webp{_file_v(stem + '_mid.webp')}",
|
||
'full': f"{sub_dir_name}/{stem}{ext}{_file_v(stem + ext)}",
|
||
}
|
||
except Exception as e:
|
||
logger.warning(f"Sub-tiling abandoned for {tile['dir_name']}: {e}")
|
||
return None
|
||
|
||
usable = [e for e in entries.values() if e['viz']]
|
||
if not usable:
|
||
return None
|
||
# Whole-tile fallback for the visualizations outside the selection: they
|
||
# still work as layers (heavier images, but functional).
|
||
sub_keys = set(offered_viz_keys)
|
||
for viz_key, info in tile['viz'].items():
|
||
if viz_key not in sub_keys:
|
||
_fallback_full_dalle(entries, viz_key, info)
|
||
return usable
|
||
|
||
|
||
def _viz_src_dir(tile, info):
|
||
"""Actual directory of a visualization's file.
|
||
|
||
Layers merged from another resolution (run interrupted between the two
|
||
passes) live in their original directory — info['dir_name'] — and not
|
||
in the dir_path of the displayed tile.
|
||
"""
|
||
d = info.get('dir_name') or tile.get('dir_name')
|
||
if d and d != Path(tile['dir_path']).name:
|
||
return Path(tile['dir_path']).parent / d
|
||
return Path(tile['dir_path'])
|
||
|
||
|
||
def _collect_tile_metadata(tile, dtm_dir):
|
||
"""Gather the generation metadata of a tile.
|
||
|
||
Reads the ground classification method from the DTM sidecar
|
||
(output/DTM/{basename}_dtm{suffix}_method.txt, written by pipeline.py),
|
||
and the dates/sizes of the visualization files.
|
||
|
||
Returns:
|
||
{method: str|None, generated: str|None,
|
||
viz: {viz_key: {date: str, size: int}}}
|
||
"""
|
||
meta = {'method': None, 'generated': None, 'viz': {}}
|
||
suffix = _res_suffix_str(tile['resolution'])
|
||
method_file = Path(dtm_dir) / f"{tile['basename']}_dtm{suffix}_method.txt"
|
||
if not method_file.exists() and suffix:
|
||
# Ground classification is shared across resolutions: fall back to
|
||
# the primary resolution's sidecar if the specific one is missing.
|
||
method_file = Path(dtm_dir) / f"{tile['basename']}_dtm_method.txt"
|
||
try:
|
||
if method_file.exists():
|
||
method = method_file.read_text(encoding='utf-8').strip()
|
||
if method:
|
||
meta['method'] = method
|
||
# The sidecar is written right after the DTM is created:
|
||
# its date ≈ the tile's generation date.
|
||
meta['generated'] = datetime.fromtimestamp(
|
||
method_file.stat().st_mtime).strftime('%Y-%m-%d %H:%M')
|
||
except OSError as e:
|
||
logger.debug(f"Unreadable metadata {method_file.name}: {e}")
|
||
|
||
for viz_key, info in tile['viz'].items():
|
||
try:
|
||
viz_dir = _viz_src_dir(tile, info)
|
||
st = (viz_dir / info['filename']).stat()
|
||
meta['viz'][viz_key] = {
|
||
'date': datetime.fromtimestamp(st.st_mtime).strftime('%Y-%m-%d %H:%M'),
|
||
'size': st.st_size,
|
||
}
|
||
except OSError:
|
||
continue
|
||
|
||
if meta['generated'] is None and meta['viz']:
|
||
dates = [v['date'] for v in meta['viz'].values()]
|
||
meta['generated'] = min(dates)
|
||
return meta
|
||
|
||
|
||
def build_index(output_dir, output_format='avif'):
|
||
"""Rebuild the catalog of processed tiles: thumbnails + inventory.
|
||
|
||
Scans output_dir/visualisations/, collects the generation metadata,
|
||
generates the JPEG thumbnails (index_thumbs/) and the sub-tiles
|
||
(index_subtiles/) — source levels of the XYZ pyramid (tiles.py) — then
|
||
writes the inventory output/index_tiles.json (served by /api/tiles).
|
||
|
||
Args:
|
||
output_dir: root output directory (contains visualisations/).
|
||
output_format: image format ('avif' or 'webp') — unused, kept for
|
||
call compatibility.
|
||
|
||
Returns:
|
||
Path to index_tiles.json on success, None on failure or if there is no tile.
|
||
"""
|
||
output_dir = Path(output_dir)
|
||
vis_dir = output_dir / 'visualisations'
|
||
dtm_dir = output_dir / 'DTM'
|
||
t_start = time.time()
|
||
tiles = scan_tiles(vis_dir)
|
||
|
||
if not tiles:
|
||
logger.info("No processed tile found — global index not generated")
|
||
return None
|
||
|
||
# GPS bounds per tile (exact georeferencing for the Leaflet map)
|
||
attach_gps_bounds(tiles)
|
||
|
||
# A single tile per position (col, row): keep the finest available
|
||
# resolution. Otherwise the 0.5 m and 0.2 m versions of the same tile
|
||
# would overlap exactly on the map and one would hide the other.
|
||
best_by_pos = {}
|
||
tiles_by_pos = {}
|
||
for t in tiles:
|
||
key = (t['col'], t['row'])
|
||
tiles_by_pos.setdefault(key, []).append(t)
|
||
if key not in best_by_pos or t['resolution'] < best_by_pos[key]['resolution']:
|
||
best_by_pos[key] = t
|
||
# Complete each displayed tile (finest resolution) with the
|
||
# visualizations produced only at the other resolution: otherwise a layer
|
||
# being generated (0.5 m pass done, 0.2 m not yet) would stay invisible
|
||
# on the map and missing from the layer list.
|
||
# dir_name records the original directory: thumbnails, URLs and metadata
|
||
# must read the file where it actually exists.
|
||
for key, best in best_by_pos.items():
|
||
for other in tiles_by_pos[key]:
|
||
if other is best:
|
||
continue
|
||
for viz_key, viz_info in other['viz'].items():
|
||
if viz_key not in best['viz']:
|
||
best['viz'][viz_key] = dict(viz_info, dir_name=other['dir_name'])
|
||
tiles = sorted(best_by_pos.values(),
|
||
key=lambda t: (t['resolution'], -t['row'], t['col']))
|
||
|
||
# Detect the geographic zones (for information only)
|
||
zones = compute_zones(tiles)
|
||
logger.info(f" {len(zones)} zone(s) detected")
|
||
|
||
thumb_dir = output_dir / 'index_thumbs'
|
||
thumb_dir.mkdir(parents=True, exist_ok=True)
|
||
|
||
# Collect every available visualization (for the layer list).
|
||
all_viz_keys = set()
|
||
for t in tiles:
|
||
all_viz_keys.update(t['viz'].keys())
|
||
|
||
# Restrict the sub-tile cutting to the chosen visualizations
|
||
if _CARTO_SUBTILED_VIZ:
|
||
sub_viz = [v for v in _CARTO_SUBTILED_VIZ if v in all_viz_keys]
|
||
if sub_viz:
|
||
logger.info(f" Sub-tiling limited to: {', '.join(sub_viz)}")
|
||
else:
|
||
sub_viz = list(all_viz_keys)
|
||
|
||
# Generate the thumbnails and build the inventory records.
|
||
zone_records = []
|
||
thumbs_generated = 0
|
||
thumbs_failed = 0
|
||
tile_idx = 0
|
||
n_tiles = len(tiles)
|
||
logger.info(f" Thumbnails: {n_tiles} tile(s) × {len(all_viz_keys)} visualization(s)")
|
||
for zone in zones:
|
||
zone_tile_records = []
|
||
for t in zone['tiles']:
|
||
tile_idx += 1
|
||
viz_thumbs = {}
|
||
regen = 0
|
||
for viz_key, info in t['viz'].items():
|
||
# Layers merged from another resolution live in their
|
||
# original directory (info['dir_name']), not dir_path.
|
||
src = _viz_src_dir(t, info) / info['filename']
|
||
src_mtime = _mtime(src)
|
||
thumb_name = f"{t['dir_name']}_{viz_key}.jpg"
|
||
thumb_path = thumb_dir / thumb_name
|
||
mid_name = f"{t['dir_name']}_{viz_key}_mid.jpg"
|
||
mid_path = thumb_dir / mid_name
|
||
# Regenerate if missing or stale (tile recomputed since)
|
||
if (not _cached_file_fresh(thumb_path, src_mtime)
|
||
or not _cached_file_fresh(mid_path, src_mtime)):
|
||
if generate_thumbnail(src, thumb_path, mid_path=mid_path,
|
||
mid_size=_MID_THUMB_SIZE):
|
||
thumbs_generated += 1
|
||
regen += 1
|
||
else:
|
||
thumbs_failed += 1
|
||
continue
|
||
else:
|
||
thumbs_generated += 1
|
||
viz_dir_name = _viz_src_dir(t, info).name
|
||
# Each URL versioned by the mtime of ITS file (see
|
||
# _url_version): immutable browser cache possible.
|
||
viz_thumbs[viz_key] = {
|
||
'thumb': f"index_thumbs/{thumb_name}"
|
||
f"{_url_version(_mtime(thumb_path))}",
|
||
# The source tile IS the served file: its mtime is enough.
|
||
'full': f"visualisations/{viz_dir_name}/{info['filename']}"
|
||
f"{_url_version(src_mtime)}",
|
||
}
|
||
if mid_path.is_file():
|
||
viz_thumbs[viz_key]['mid'] = (
|
||
f"index_thumbs/{mid_name}{_url_version(_mtime(mid_path))}")
|
||
|
||
if regen:
|
||
logger.info(f" [{tile_idx}/{n_tiles}] {t['dir_name']} — "
|
||
f"{regen} thumbnail(s) regenerated")
|
||
|
||
if not viz_thumbs:
|
||
continue
|
||
|
||
display_viz = _pick_display_viz(viz_thumbs.keys())
|
||
tile_meta = _collect_tile_metadata(t, dtm_dir)
|
||
zone_tile_records.append({
|
||
'col': t['col'],
|
||
'row': t['row'],
|
||
'name': t['basename'],
|
||
'dir_name': t['dir_name'],
|
||
'resolution': t['resolution'],
|
||
'bounds': t.get('bounds'),
|
||
'corners': t.get('corners'),
|
||
'display_viz': display_viz,
|
||
'viz': viz_thumbs,
|
||
'meta': tile_meta,
|
||
})
|
||
|
||
if zone_tile_records:
|
||
zone_records.append({
|
||
'label': zone['label'],
|
||
'tiles': zone_tile_records,
|
||
'bbox': zone['bbox'],
|
||
})
|
||
|
||
if not zone_records:
|
||
logger.warning("No thumbnail generated — global index abandoned")
|
||
return None
|
||
|
||
# Compute the global bbox (for the summary log below)
|
||
global_bbox = compute_bbox(tiles)
|
||
|
||
# Flat list of displayable quads for the inventory.
|
||
# 0.2 m tiles (5000×5000 px) are cut into 500 m sub-tiles
|
||
# (2500×2500 px quadrants) to lighten memory use and loading.
|
||
sub_dir_name = 'index_subtiles'
|
||
display_tiles = []
|
||
n_dalles = 0
|
||
n_sous = 0
|
||
for zr in zone_records:
|
||
for t in zr['tiles']:
|
||
n_dalles += 1
|
||
subs = _build_subtiles(t, sub_viz, output_dir, sub_dir_name)
|
||
if subs:
|
||
display_tiles.extend(subs)
|
||
n_sous += len(subs)
|
||
else:
|
||
display_tiles.append(t)
|
||
if n_sous:
|
||
logger.info(f" {n_sous} sub-tile(s) generated for {n_dalles} tile(s)")
|
||
|
||
# Tile inventory (index_tiles.json): mapserve serves it via /api/tiles
|
||
# to the lightweight machines (LIDAR_SOURCE_URL); the pipeline rewrites
|
||
# it after every tile by default (unless --no-index).
|
||
viz_meta = {k: {'label': VIZ_LABELS.get(k, k)}
|
||
for k in _VIZ_FALLBACK_ORDER if k in all_viz_keys}
|
||
for k in sorted(all_viz_keys):
|
||
viz_meta.setdefault(k, {'label': VIZ_LABELS.get(k, k)})
|
||
from .quality import load_quality_table
|
||
inventory_path = output_dir / 'index_tiles.json'
|
||
inventory_path.write_text(json.dumps({
|
||
'tiles': display_tiles,
|
||
'viz_meta': viz_meta,
|
||
'stats': {'n_tiles': len(display_tiles)},
|
||
# Per-tile quality (PDF export inset): copied to the lightweight
|
||
# machines so it stays available when the upstream is down.
|
||
'quality': load_quality_table(output_dir),
|
||
}, ensure_ascii=False), encoding='utf-8')
|
||
|
||
logger.info(f"Inventory generated: {inventory_path} "
|
||
f"({time.time() - t_start:.1f}s)")
|
||
logger.info(f" {len(tiles)} tile(s) • {thumbs_generated} thumbnail(s) generated"
|
||
+ (f" • {thumbs_failed} failure(s)" if thumbs_failed else ""))
|
||
logger.info(f" Grid: {global_bbox['min_col']}-{global_bbox['max_col']} km E × "
|
||
f"{global_bbox['min_row']}-{global_bbox['max_row']} km N")
|
||
|
||
return inventory_path
|