From 90dbe61f3b2b0f7334229d25a727c5f33344e544 Mon Sep 17 00:00:00 2001 From: Antoine Jacquin Date: Sun, 27 Sep 2026 15:07:04 +0200 Subject: [PATCH] =?UTF-8?q?Borner=20le=20comblement=20du=20MNT=20=C3=A0=20?= =?UTF-8?q?l'enveloppe=20des=20points,=20ajouter=20la=20couche=20pr=C3=A9c?= =?UTF-8?q?ision=20et=20un=20affichage=20relief/pr=C3=A9cision,=20encoder?= =?UTF-8?q?=20les=20sous-tuiles=20une=20seule=20fois=20en=20q75?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Opus 5.5 --- AGENTS.md | 8 +- docs/MAPS.md | 61 +-- lidar_pipeline/dtm.py | 222 ++++++++- lidar_pipeline/index.py | 213 +++++---- lidar_pipeline/mapserve.py | 77 ++- lidar_pipeline/mapui.py | 501 ++++++++------------ lidar_pipeline/pipeline.py | 18 +- lidar_pipeline/rendering.py | 49 +- lidar_pipeline/tests/test_dtm.py | 99 ++++ lidar_pipeline/tests/test_index.py | 18 +- lidar_pipeline/tests/test_mapserve.py | 155 +++--- lidar_pipeline/tests/test_pipeline.py | 12 +- lidar_pipeline/tests/test_rendering.py | 49 ++ lidar_pipeline/tests/test_tiles.py | 25 + lidar_pipeline/tests/test_visualizations.py | 27 ++ lidar_pipeline/tiles.py | 11 +- lidar_pipeline/visualizations.py | 42 ++ 17 files changed, 994 insertions(+), 593 deletions(-) diff --git a/AGENTS.md b/AGENTS.md index bb54896..50aa880 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -19,10 +19,10 @@ ## Conventions - **Un seul serveur web : `mapserve.py`** (image `lidar-maps`, port 8975 léger / 8973 worker). L'ancienne webapp (`webapp.py`, `export.py`, index.html/`_APP_JS`) a été supprimée : la génération de tuiles (portée de la webapp historique — `/api/preview`, `/api/generate`, `/api/status`, `/api/stop`, `/api/queue/clear`, `/api/cell`) vit dans `mapserve.py`, l'interface dans `mapui.py` (constantes `_MAP_HTML`/`_MAP_CSS`/`_MAP_JS`, écrites par `write_map_assets()` et bâchées dans les images). Sur l'image légère sans `LIDAR_GENERATION_URL`, `/api/status` répond `available: false` et l'interface masque les boutons. -- **`index.py` = catalogue + registres partagés** (plus d'interface) : `VIZ_LABELS`/`VIZ_LEGENDS`, défauts d'affichage (`DEFAULT_LAYERS`/`DEFAULT_OPACITY`/`DEFAULT_BLEND`), `PANEL_VIZ`/`KEYWORD_TO_STEP`, `scan_tiles`/`cells_with_all_viz`, vignettes + sous-tuiles + inventaire `index_tiles.json` (`build_index`). L'inventaire est servi par `/api/tiles` de mapserve aux machines légères (`LIDAR_SOURCE_URL`). +- **`index.py` = catalogue + registres partagés** (plus d'interface) : `VIZ_LABELS`/`VIZ_LEGENDS`, défauts d'affichage (`DEFAULT_VIZ`/`PRECISION_VIZ`/`VIEW_MODES`), `PANEL_VIZ`/`KEYWORD_TO_STEP`, `scan_tiles`/`cells_with_all_viz`, vignettes + sous-tuiles + inventaire `index_tiles.json` (`build_index`). L'inventaire est servi par `/api/tiles` de mapserve aux machines légères (`LIDAR_SOURCE_URL`). - **Generation is 0.2 m only** (policy): `/api/generate` (`GENERATE_RESOLUTIONS` in `mapserve.py`), the compose `process` command and the CLI `-r` default all produce 0.2 m exclusively; 0.5 m stays available via explicit `-r 0.5`. Completeness detection (`complete_cells`) requires the viz at 0.2 m only. - **Génération du nord au sud** : les tuiles sont traitées par ligne décroissante (row = nord en km), colonnes croissantes — `find_laz_files` (pipeline.py) pour les passes batch et `_resolve_request` (mapserve.py) pour les runs lancés depuis la carte. Les workers prennent les fichiers dans l'ordre de soumission : la carte se remplit de haut en bas pendant un run (`--file` explicite au CLI = ordre utilisateur préservé). Parallélisme de génération : `LIDAR_WORKERS` (10 dans les compose ; `auto` sinon). -- **Sub-tuilage intégral** : `_CARTO_SUBTILED_VIZ` (vide dans `index.py`) découpe TOUTES les couches en quadrants 500 m à 0,2 m ; ortho/topo sont encodées en AVIF q75 (`_SUBTILE_DETAIL_VIZ`) contre q55 pour les rampes de couleur. Une couche qui échoue à la découpe retombe en dalle entière (`_fallback_full_dalle`) sans pénaliser les autres. +- **Sub-tuilage intégral** : `_CARTO_SUBTILED_VIZ` (vide dans `index.py`) découpe TOUTES les couches en quadrants 500 m à 0,2 m ; toutes en AVIF 4:2:0 q75 (`_SUBTILE_AVIF_QUALITY`), encodées UNE fois depuis le raster d'origine : `tif_to_crop` appelle `index.write_subtiles` juste après la dalle (plus récentes qu'elle, `build_index` ne les réencode pas ; l'ancienne chaîne dalle q60 → sous-tuile q55 cumulait deux pertes : 18,1 dB/SSIM 0,84 contre 19,5 dB/0,93). Le 4:2:0 plafonne vers 19,5 dB sur le relief (teinte pixel par pixel moyennée par 2 × 2) : seul le 4:4:4 irait plus loin (q75 : 28 dB, ~1,6× plus lourd). Exception : les aplats de niveaux (`_SUBTILE_LOSSLESS_VIZ`, la précision) en WebP sans perte (`.webp`, accepté par `tiles.py`). **Pillow ignore `lossless=True` en AVIF** (q75 avec perte) ; en RGB aucun réglage AVIF n'est exact. Une couche qui échoue à la découpe retombe en dalle entière (`_fallback_full_dalle`) sans pénaliser les autres. - **Tuiles XYZ réutilisables hors du projet** : `/tiles/{layer}/{z}/{x}/{y}.png` suit le schéma OpenStreetMap (256 px, EPSG:3857, PNG RGBA transparent hors emprise, CORS `*`) — JOSM, iD, QGIS, uMap et MapLibre le consomment tel quel, avec découverte via TileJSON / WMTS / `josm.imagery.xml`. Le 512 px (`@2x.webp`) est réservé à l'interface interne (moitié moins de requêtes en HTTP/1.1). Toute évolution du gabarit d'URL casse des configurations clientes : la changer demande une décision explicite. - **Bilingual naming**: all code identifiers are English; every user-facing string, log message, argparse help, and comment is French. - **Adding a visualization requires 4 edits**: (1) `generate_X()` in `visualizations.py`, (2) entry in `VIZ_STEPS` in `pipeline.py`, (3) entry in `COLORMAPS` in `rendering.py` (or `RGB_LEGENDS` for an RGB output), (4) entry in `VIZ_LEGENDS` in `index.py` (title/legend/description + sampled cmap gradient — single text source merged into `COLORMAPS` at import, also served in `/api/map/meta` and le TileJSON). Missing any one breaks the pipeline. @@ -30,6 +30,7 @@ - **Return `None` on failure, never raise**: `dtm.py`, `visualizations.py`, and `ign.py` all return `None` to let the pipeline continue. Raising aborts the entire file. - **Logger is always `logging.getLogger("lidar")`**, never `__name__`. All modules route through this single logger so worker processes can configure it. - **Ajustement conjoint des lignes de balayage (3ᵉ passe du calage, `STRIP_ALIGN_VERSION` 3)** : chaque ligne de balayage (~150 lignes/s, découpée par `_scan_line_ids` : retour de la dent de scie de `scan_angle`) peut être décalée ET inclinée (roulis) de 1 à 15 cm — lignes en creux isolées, passe entière basculée visible en bande (mesuré sur LHD_FXX_0999_6882). `_joint_line_corrections` (`dtm.py`) recale chaque ligne (a + b·u) contre le consensus des AUTRES faisceaux de sa maille 1 m (plan local par la pente), moindres carrés tronqués à 5 MAD par `bincount`, pas amorti 0,5, arrêt sous 1 mm, recentrage (pas de dérive d'ensemble), estimation sur 1 point sur 3 ; CuPy si GPU actif, repli numpy. S'y ajoute un **profil d'étalonnage par faisceau et par degré d'angle** (`STRIP_ANGLE_BIN`, commun à toutes les lignes du faisceau, moyenne retirée, pente conservée car l'inclinaison par ligne est indéterminable quand la fauchée ne traverse la dalle qu'en partie ; classes pauvres = valeur voisine) : il retire les écarts non linéaires en travers de la fauchée (±1,5 cm mesurés au bord), sources de marches parallèles au vol au bord de fauchée. Lignes sans recouvrement : `_scan_line_corrections_beam` (contre la surface de leur propre faisceau, composante ligne à ligne seule, σ 3 lignes). Validé sur blocs jamais vus (damier 10 m). La gigue par fenêtres de temps (2ᵉ passe) n'est plus calculée quand `scan_angle` existe (redondante). Coût CPU ~26 s par dalle pour tout le calage (73 s avant). Sidecar : `lines`, `line_window`, `line_cell`, `line_model`. +- **Comblement des vides entre points borné à l'enveloppe (`GAP_FILL_VERSION` 2, `_fill_small_gaps` dans `dtm.py`)** : à 0,2 m ~80 % des pixels n'ont aucun point. L'ancien `fillnodata` à 1 m comblait tout pixel à moins de 1 m d'un point — chaque point isolé devenait une pastille plate de 2 m (teinte `atan2(0,0)` = rose saturé dans le relief orienté) et chaque trou recevait une bande extrapolée de 1 m (liseré coloré). Désormais : fermeture morphologique des pixels mesurés (rien n'est étendu vers l'extérieur) dont le rayon suit l'espacement local des points (`GAP_RADIUS_K` 1,5 × espacement mesuré sur 5 m, paliers `GAP_RADII_M` 1/1,5/2/3 m, 1 m mini pour boucher les trous de voitures), puis îlots < `GAP_MIN_ISLAND_M2` (1 m²) retirés. Morphologie par tranches numpy (`_morph_step`, ~10× scipy) ; ~5–9 s par dalle 5000². Version dans le tag GeoTIFF `LIDAR_GAP_FILL` (3 : + fichier annexe de densité) : absent/différent ⇒ DTM régénéré (`_gap_fill_matches` dans `pipeline.py`). `rasterio.fill.fillnodata` écrit dans son entrée : toujours lui passer une copie. - **Openness à échelle fixe** : `generate_openness` normalise par des références figées `OPENNESS_POS_REF` / `OPENNESS_NEG_REF` (degrés, médianes mesurées sur 15 dalles réparties sur le territoire) et plus par z-score de dalle — même ouverture = même couleur, mosaïque jointive. - **Filename special-cases** in `_expected_output_path()`: `pos_open` → `positive_openness`, `neg_open` → `negative_openness`, `hillshade` → `hillshade_multi`. - **Génération imposée** : la carte ne propose plus de réglage de classification ni de raccord — `_build_command` (`mapserve.py`) passe toujours `--ground-classification ign --ign-classes sol --edge-buffer 100` ; les champs `ground_class`/`ign_classes`/`reclassify`/`edge_buffer` éventuellement envoyés sont ignorés. Défauts CLI alignés (`ign`, `100`). @@ -47,7 +48,8 @@ - **Default output is AVIF**, not WebP. Use `--format webp` for WebP. Quality default is 60 (visually lossless on smooth color ramps, ~÷3 vs q98). - **Calage vertical des faisceaux de vol** : chaque tuile mélange plusieurs passes (1-2 `PointSourceId` par passe) parfois biaisées verticalement de quelques cm (±2,5 cm mesurés sur 1000_6882). `create_dtm_fast` mesure l'offset robuste de chaque faisceau (points sol, maille 1 m, surface médiane itérée 3×) et retranche les offsets ≥ 0,5 cm (`STRIP_ALIGN_THRESHOLD` dans `dtm.py`) avant rastérisation. Offsets calculés **par tuile** (ils dérivent le long d'une ligne de vol : jamais de table globale), mémoïsés par LAS sol, consignés dans `DTM/*_dtm*_stripalign.json` (sidecar de cache : absent, ou version/seuil/paramètres différents ⇒ régénération du DTM). Désactivable : `--no-strip-align`. - **Gigue intra-faisceau (2ᵉ passe du calage)** : les lignes de balayage successives d'une MÊME passe peuvent être décalées verticalement de façon aléatoire (vibration capteur / bruit haute fréquence de trajectoire) — un offset constant par faisceau n'y suffit pas. `_strip_jitter_offsets` découpe chaque faisceau en fenêtres de temps GPS (`STRIP_JITTER_BIN` = 0,1 s, origine de temps propre à chaque faisceau), mesure l'offset robuste de chaque fenêtre contre la surface médiane des AUTRES faisceaux (maille 1 m partagée, ≥ `STRIP_JITTER_MIN_CELLS` = 40 cellules), lisse la série (médiane glissante `STRIP_JITTER_SMOOTH` = 5 fenêtres), la borne à ± `STRIP_JITTER_MAX` (10 cm) puis l'interpole au temps GPS de chaque point (`_apply_strip_jitter`) ; en recouvrement à deux faisceaux, chacun reçoit une série (chacun absorbe sa part). Requiert la dimension `gps_time` (silencieusement ignorée sinon). Sidecar version 2 (séries dans `jitter`), couverte par `--no-strip-align`. -- **Couches produites et servies = `PANEL_VIZ` (`index.py`, aujourd'hui `('relief_oriente',)`)** : pipeline sans `--only`/`--skip` (`panel_steps()`), génération lancée depuis la carte (`_panel_viz_steps` dans `mapserve.py`) et couches servies (`tiles.available_layers` filtre : panneau, `/tiles/…`, TileJSON, WMTS, JOSM). Les autres visualisations restent calculables avec `--only`. Les tests de la carte lèvent la restriction via `_setup(..., panel=None)`. +- **Couches produites et servies = `PANEL_VIZ` (`index.py`, aujourd'hui `('relief_oriente', 'densite_sol')`)** : pipeline sans `--only`/`--skip` (`panel_steps()`), génération lancée depuis la carte (`_panel_viz_steps` dans `mapserve.py`) et couches servies (`tiles.available_layers` filtre : panneau, `/tiles/…`, TileJSON, WMTS, JOSM). Les autres visualisations restent calculables avec `--only`. Les tests de la carte lèvent la restriction via `_setup(..., panel=None)`. +- **Précision (`densite_sol`) et affichage de la carte** : `create_dtm_fast` écrit la densité des points sol retenus (bande de raccord comprise) dans `DTM/*_dtm*_density.tif` (mailles 1 m, moyenne 3 × 3, pts/m²) ; `generate_densite_sol` la quantifie en 16 niveaux (`density_levels` : échelle log fixe, niveau k dès 0,25 × 2^(k/2) pts/m²) et la garde à 1 m (dalle 1000², ~200 Ko en WebP sans perte, `rendering.LOSSLESS_GRAY_KEYWORDS`) ; tuiles au plus proche voisin (`tiles.NEAREST_LAYERS`) pour garder 16 gris nets aux zooms 18–19. Interface (`mapui.py`) : plus de pile de couches — une couche principale (`DEFAULT_VIZ`) et trois modes relief / précision / les deux (`VIEW_MODES`, précision en `multiply` à `DEFAULT_PRECISION_OPACITY` 0,6), touche P, lien `&M=…&P=mode:opacité`, défauts figés `{main, mode, precision_opacity, base}`. - **Relief orienté (`relief_oriente`, seule couche de la carte)** : image RGB unique (GeoTIFF uint8 3 bandes, rendue telle quelle comme ortho/topo via `RGB_KEYWORDS` dans `rendering.py`) — clarté CIELAB = openness positive **locale** (MNT − gaussienne `RELIEF_DETREND_M` = 10 m, rayons `RELIEF_RADII_M` = 5/10/20 m, 16 directions) 65 % + ombrage 35 % ; teinte = aspect, chroma fixe (`RELIEF_CHROMA`). Échelle log fixe `RELIEF_OPEN_RANGE` (pas de statistique par dalle) et support total 40 m < bande de raccord 100 m : dalles jointives. Rapide : détendance + rayons sur grille décimée ~0,8 m (`RELIEF_GRID_M`), noyau dédié qui n'accumule que la moyenne des angles (`_mean_horizon_*` : CuPy `RawKernel` sur GPU, numba parallèle sur CPU, numpy en repli), colorisation fusionnée (numba) ou vectorisée sans trigonométrie (CuPy) via une table L* × teinte (`_relief_lut`). ~5 s de calcul hors préparation sur CPU 12 cœurs. Tout changement de constante change le rendu : régénérer les dalles (`--only relief_oriente --force`). - **Openness sous-échantillonnée** : `generate_openness` calcule le lancé de rayons (l'étape la plus coûteuse : 532 s/tuile à 0,2 m sur CPU) sur une grille décimée par blocs (`OPENNESS_DOWNSAMPLE = 2` : max par bloc en positive, min en négative — préserve les reliefs qui bornent l'horizon) puis rééchantillonne en bilinéaire. Coût ÷ facteur³ : 532 s → 40 s (×13). Signal archéologique préservé (corr. 0,93 après lissage) ; la texture de bruit sub-métrique disparaît. `--openness-downsample 1` = pleine résolution. SVF et openness anisotrope ne sont PAS concernées. - **Raccord des bords entre tuiles** : les rendus à grand noyau (openness/SVF : rayons 100 m ; LRM : 15 m) tronquent leur fenêtre au bord de dalle — bandes d'artefacts à chaque changement de tuile. `--edge-buffer N` (défaut 100 ; toujours appliqué par la génération depuis la carte, `EDGE_BUFFER_METERS` = 100 m dans `mapserve.py`, plus de case à cocher ; 0 = off en ligne de commande) fait rastériser le DTM sur la **dalle nominale 1 km alignée sur la grille** plus une bande de N m remplie avec les points sol des 8 LAZ voisines (`_neighbor_ground_points` dans `dtm.py` : PDAL en flux, découpe + filtre de classes IGN ; voisine absente = téléchargement automatique depuis le catalogue IGN avant le run, **isolée dans `input/edge_neighbors/`** pour ne pas gonfler le corpus des passes globales (dédupliqué sur tout le lot, `_fetch_edge_neighbors` dans `pipeline.py`) ; introuvable ou échec = bande vide). Les visualisations calculent sur l'emprise étendue puis `rendering.py` (`_core_tile_window`, via `tif_to_crop`/`tif_to_png`) recadre les sorties sur la dalle 1 km exacte lue dans le nom LHD — les AVIF restent des carrés 1 km alignés dans la mosaïque. Tampon consigné dans le tag GeoTIFF `LIDAR_EDGE_BUFFER` du DTM : changer `--edge-buffer` invalide le cache DTM automatiquement (tag absent = 0). Bandes voisines non calées par faisceaux (contexte seul, recadrée hors image finale). Coût : ~7 s de lecture par voisine + ~44 % de pixels en plus à 100 m/0,2 m. Nom hors pattern LHD : option ignorée (bornes d'en-tête, pas de recadrage). diff --git a/docs/MAPS.md b/docs/MAPS.md index 18f3278..b7772e9 100644 --- a/docs/MAPS.md +++ b/docs/MAPS.md @@ -120,34 +120,38 @@ Quand le relief orienté est affiché, une rose des vents donne la couleur de chaque orientation de pente (même formule CIELAB que le rendu) ; la clarté porte le relief local (clair = bosse, sombre = creux). -## Pile de couches +## Affichage : relief et précision -La carte ne sert que les couches de `PANEL_VIZ` (`index.py`) : aujourd'hui le -seul **relief orienté**, qui fusionne openness locale et orientation des pentes -en une image. Les autres visualisations présentes sur disque ne sont ni listées -ni servies en tuiles. +La carte ne sert que les couches de `PANEL_VIZ` (`index.py`) : le **relief +orienté** (couche d'affichage principal, qui fusionne openness locale et +orientation des pentes) et la **précision** (`densite_sol` : densité des points +sol retenus pour le MNT). Les autres visualisations présentes sur disque ne +sont ni listées ni servies en tuiles. -Chaque couche du panneau porte trois réglages, tous persistés et transportés -par le lien de partage : +Il n'y a plus de pile de couches. Le panneau propose trois modes, un clic +chacun, ou la touche **P** pour passer au suivant : -- **ordre** — poignée ⠿ en glisser-déposer, ou boutons ▲▼ (les seuls - utilisables au doigt). La liste est encadrée par « Haut de pile — devant » et - « Bas de pile — derrière ». Pendant un glisser, la destination est explicite : - la ligne tirée s'estompe et une **barre d'insertion** marque le point de - chute, au-dessus ou en dessous de la ligne survolée selon la moitié visée ; - à l'arrivée, la couche déplacée clignote brièvement et est ramenée dans le - champ de vision ; -- **opacité** — curseur par couche ; -- **mode de fusion** — `mix-blend-mode` CSS : Normal, Produit, Écran, - Superposé, Doux, Lumière crue, Différence, Luminosité. La superposition ne se - réduit donc pas à de la transparence : « Produit » pose un ombrage sur une - rampe de couleur sans la délaver, « Superposé » creuse le contraste local, - « Différence » fait ressortir les écarts entre deux couches. +- **Relief** — le relief orienté seul ; +- **Précision** — la densité seule, en 16 gris (échelle log fixe : niveau k à + partir de 0,25 × 2^(k/2) pts/m², noir ≤ 0,35 ou aucun point, blanc ≥ 45) ; + se lit comme une carte de fiabilité géométrique ; +- **Les deux** — la précision en « produit » (`mix-blend-mode: multiply`) sur + le relief, opacité réglable (60 % par défaut) : les zones où le relief est + interpolé s'assombrissent sans masquer le relief. + +La légende de la précision (16 paliers, info-bulle en pts/m² sur chaque +palier) s'affiche dès que la précision est visible. Les couches LiDAR vivent +dans un conteneur **isolé** (`isolation: isolate`) : le produit n'agit que sur +le relief, jamais sur le fond de carte. + +Le lien de partage transporte la couche principale, le mode et l'opacité : +`#z/lat/lng&M=relief_oriente&P=both:60&B=1:85:1`. Les anciens liens de la pile +(`&L=…`) s'ouvrent sans erreur sur l'état courant. ### Figer la configuration -Le bouton **★ Définir par défaut** enregistre l'état courant — ordre, couches -allumées, opacités, modes de fusion, fond de carte et fusion de la pile — dans +Le bouton **★ Définir par défaut** enregistre l'affichage courant — couche +principale, mode, opacité de la précision, fond de carte — dans `output/.map-defaults.json`. Tout navigateur sans réglage local part alors de cette configuration ; **↺ Réinitialiser** oublie l'état local et y revient. @@ -157,14 +161,11 @@ curl -X DELETE http://localhost:8975/api/map/defaults # retour au registre ``` Sans fichier enregistré, les défauts viennent du registre du pipeline -(`DEFAULT_LAYERS`, `DEFAULT_OPACITY`, `DEFAULT_BLEND` dans `index.py`). Les -valeurs reçues sont filtrées : couches inconnues écartées, opacités bornées à -0–1, modes de fusion validés. - -Les couches LiDAR vivent dans un conteneur **isolé** (`isolation: isolate`) : -les fusions agissent entre elles, jamais sur le fond de carte à travers les -zones sans donnée. La fusion de la **pile entière sur le fond** se règle à -part, en bas du panneau. +(`DEFAULT_VIZ`, `PRECISION_VIZ`, `DEFAULT_VIEW_MODE`, +`DEFAULT_PRECISION_OPACITY` dans `index.py`). Les valeurs reçues sont +filtrées : couche inconnue (ou la précision elle-même) refusée comme +principale, mode validé, opacités bornées à 0–1. Un fichier de l'ancienne pile +(`order`/`on`/`blend`) est ignoré, sauf le fond. > **Licence** — LiDAR HD est diffusé sous **Licence Ouverte 2.0** : > l'attribution IGN est obligatoire et doit rester visible chez le client. diff --git a/lidar_pipeline/dtm.py b/lidar_pipeline/dtm.py index 7407cb5..1e077e9 100644 --- a/lidar_pipeline/dtm.py +++ b/lidar_pipeline/dtm.py @@ -1446,6 +1446,200 @@ def _repair_laz_with_laspy(input_laz, output_las): return False +# ------------------------------------------------------------ +# Comblement des vides entre points (petits trous) +# ------------------------------------------------------------ +# À 0,2 m, ~80 % des pixels ne reçoivent aucun point : le MNT est comblé entre +# les points mesurés. Un comblement « à distance fixe » (tout pixel à moins de +# 1 m d'un point) faisait grossir chaque point isolé en pastille plate de 2 m +# et extrapolait une bande de 1 m dans chaque trou (plateaux recopiés, fausses +# pentes : pastilles roses et liserés colorés dans le relief orienté). Le +# comblement est désormais borné à l'ENVELOPPE des points : fermeture +# morphologique dont le rayon suit l'espacement local des points (court en +# zone dense, long sous couvert clairsemé), puis suppression des îlots isolés. +GAP_FILL_TAG = "LIDAR_GAP_FILL" # tag GeoTIFF : version du comblement +GAP_FILL_VERSION = 3 # 1 (absent) = distance fixe 1 m ; 3 = + densité +GAP_RADIUS_K = 1.5 # rayon = K × espacement local des points +GAP_RADII_M = (1.0, 1.5, 2.0, 3.0) # paliers de rayon (m) ; 1 m mini : trous < 2 m comblés +GAP_DENSITY_WINDOW_M = 5.0 # fenêtre de mesure de la densité locale +GAP_MIN_ISLAND_M2 = 1.0 # îlots de données plus petits : supprimés + + +def _morph_step(mask, square, erode): + """Un pas 3×3 (carré ou croix) de dilatation/érosion binaire par tranches + numpy (~10× plus rapide que scipy.ndimage). Au bord, le voisin manquant + compte comme le pixel lui-même.""" + op = np.logical_and if erode else np.logical_or + out = mask.copy() + op(out[1:], mask[:-1], out=out[1:]) + op(out[:-1], mask[1:], out=out[:-1]) + src = out.copy() if square else mask # carré : séparable (colonnes puis lignes) + op(out[:, 1:], src[:, :-1], out=out[:, 1:]) + op(out[:, :-1], src[:, 1:], out=out[:, :-1]) + return out + + +def _morph_disk(mask, radius_px, erode=False, first_step=1): + """Dilatation/érosion par un disque approché : octogone, pas carrés + (impairs) et croix (pairs) alternés. `first_step` poursuit une chaîne.""" + for step in range(first_step, first_step + radius_px): + mask = _morph_step(mask, step % 2 == 1, erode) + return mask + + +def _gap_fill_envelope(valid, resolution): + """Masque des pixels à combler ou conserver : enveloppe des points mesurés. + + Fermeture morphologique (dilatation puis érosion) de rayon variable : les + vides ENTRE points voisins sont comblés, rien n'est étendu vers l'extérieur + (un point isolé reste un pixel). Le rayon vaut GAP_RADIUS_K × l'espacement + local des points, mesuré sur GAP_DENSITY_WINDOW_M dans la zone couverte + (le bord d'un trou ne fait pas chuter la densité), arrondi au palier + GAP_RADII_M le plus proche. Les îlots de moins de GAP_MIN_ISLAND_M2 sont + retirés (retours épars sur l'eau, dans une cour…). + """ + from scipy import ndimage as nd + radii_px = sorted({max(1, int(round(r / resolution))) for r in GAP_RADII_M}) + r_max = radii_px[-1] + # Masque étendu en miroir : la bordure de la grille n'érode pas les + # données et ne comble pas les bandes vides. Une seule chaîne de + # dilatations sert tous les paliers (disque de rayon r = r premiers pas). + pad = 2 * r_max + dilated = {} + grown = np.pad(valid, pad, mode="symmetric") + for step in range(1, r_max + 1): + grown = _morph_disk(grown, 1, first_step=step) + if step in radii_px: + dilated[step] = grown + del grown + + # Densité locale : part de pixels mesurés dans la zone couverte (points à + # moins du plus grand rayon), convertie en espacement moyen des points. + # Calcul sur une grille de ~1 m (fenêtre de 5 m : largement suffisant). + covered = dilated[r_max][pad:-pad, pad:-pad] + h, w = valid.shape + f = max(1, int(round(1.0 / resolution))) + hc, wc = -(-h // f), -(-w // f) + + def _block_mean(a): + full = np.zeros((hc * f, wc * f), dtype=np.float32) + full[:h, :w] = a + return full.reshape(hc, f, wc, f).mean(axis=(1, 3)) + + win = max(3, int(round(GAP_DENSITY_WINDOW_M / (resolution * f))) | 1) + frac = nd.uniform_filter(_block_mean(valid), win, mode="nearest") + frac_cov = nd.uniform_filter(_block_mean(covered), win, mode="nearest") + with np.errstate(divide="ignore", invalid="ignore"): + # rayon voulu en pixels = K × espacement / résolution + wanted = GAP_RADIUS_K / np.sqrt(frac / np.maximum(frac_cov, 1e-6)) + # Rayon voulu agrandi en bilinéaire (sinon changements de palier en + # créneaux de 1 m sur le bord des trous), puis palier le plus proche. + wanted = np.nan_to_num(np.clip(wanted, 0, 2 * r_max), nan=2 * r_max) + if f > 1: + wanted = nd.zoom(wanted.astype(np.float32), f, order=1, grid_mode=True, + mode="nearest")[:h, :w] + mids = [(a + b) / 2 for a, b in zip(radii_px, radii_px[1:])] + level = np.digitize(wanted, mids).astype(np.uint8) + del covered, frac, frac_cov, wanted + + # Érosion par palier, restreinte à l'emprise de ses pixels (les grands + # rayons ne concernent que les zones clairsemées) + envelope = valid.copy() + for i, r in enumerate(radii_px): + sel = (level == i) & ~valid + rows = np.flatnonzero(sel.any(axis=1)) + if rows.size == 0: + continue + cols = np.flatnonzero(sel.any(axis=0)) + y0, y1 = rows[0], rows[-1] + 1 + x0, x1 = cols[0], cols[-1] + 1 + m = 2 * r # marge : l'érosion au bord du découpage reste hors emprise + crop = dilated[r][y0 + pad - m:y1 + pad + m, x0 + pad - m:x1 + pad + m] + closed = _morph_disk(crop, r, erode=True)[m:-m, m:-m] + envelope[y0:y1, x0:x1] |= sel[y0:y1, x0:x1] & closed + del dilated, level + + # Îlots trop petits : retirés (8-connexité, surface en m²) + labels, n = nd.label(envelope, structure=np.ones((3, 3), dtype=bool)) + if n: + area = np.bincount(labels.ravel()) * resolution * resolution + small = area < GAP_MIN_ISLAND_M2 + small[0] = False + if small.any(): + envelope &= ~small[labels] + return envelope + + +def _fill_small_gaps(dtm, resolution): + """Comble les vides entre points dans l'enveloppe des données. + + Returns: + Tuple (dtm, filled_count, removed_count) : pixels comblés et pixels + mesurés retirés avec les îlots isolés. + """ + valid = ~np.isnan(dtm) + if valid.all() or not valid.any(): + return dtm, 0, 0 + envelope = _gap_fill_envelope(valid, resolution) + from rasterio.fill import fillnodata + # Tout pixel de l'enveloppe est à moins du plus grand rayon d'un point + max_px = max(1, int(round(max(GAP_RADII_M) / resolution))) + # fillnodata écrit dans le tableau fourni : copie + filled = fillnodata(dtm.copy(), mask=valid, max_search_distance=max_px) + to_fill = envelope & ~valid & ~np.isnan(filled) + removed = valid & ~envelope + out = dtm.copy() + out[to_fill] = filled[to_fill] + out[removed] = np.nan + return out, int(to_fill.sum()), int(removed.sum()) + + +# Densité de points sol (couche « densite_sol ») : comptage des points retenus +# pour le MNT en mailles de DENSITY_CELL_M, moyenné sur DENSITY_SMOOTH × +# DENSITY_SMOOTH mailles (pts/m² sur 9 m² : moins de bruit de petits entiers). +# Écrite à côté du DTM (*_dtm*_density.tif) : les visualisations ne reçoivent +# que le MNT, plus les points. +DENSITY_CELL_M = 1.0 +DENSITY_SMOOTH = 3 + + +def density_path(dtm_path): + """Fichier annexe de densité d'un DTM.""" + dtm_path = Path(dtm_path) + return dtm_path.with_name(f"{dtm_path.stem}_density.tif") + + +def _write_density(xs, ys, bounds, dtm_path): + """Écrit la densité de points sol (pts/m², GeoTIFF float32 à 1 m).""" + from scipy import ndimage as nd + min_x, min_y, max_x, max_y = bounds + cell = DENSITY_CELL_M + nx = max(1, int(np.ceil((max_x - min_x) / cell - 1e-6))) + ny = max(1, int(np.ceil((max_y - min_y) / cell - 1e-6))) + top = max_y + counts, _, _ = np.histogram2d(top - ys, xs - min_x, bins=(ny, nx), + range=[[0, ny * cell], [0, nx * cell]]) + density = nd.uniform_filter(counts, DENSITY_SMOOTH, mode="nearest") / (cell * cell) + out = density_path(dtm_path) + with rasterio.open( + out, 'w', driver='GTiff', height=ny, width=nx, count=1, + dtype='float32', crs='EPSG:2154', + transform=from_bounds(min_x, top - ny * cell, min_x + nx * cell, top, nx, ny), + compress='deflate', + ) as dst: + dst.write(density.astype('float32'), 1) + return out + + +def read_dtm_gap_fill(dtm_path): + """Version du comblement enregistrée dans un DTM (1 si absente/illisible).""" + try: + with rasterio.open(dtm_path) as src: + return int(src.tags().get(GAP_FILL_TAG, 1) or 1) + except Exception: + return 1 + + def _interpolate_holes(dtm, downsample=8): """Fill remaining NaN holes with a terrain-aware surface interpolation. @@ -1834,6 +2028,13 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, ys = np.concatenate([ys, ny]) zs = np.concatenate([zs, nz]) + # Densité des points sol retenus (couche densite_sol), best-effort + try: + _write_density(xs, ys, (min_x, min_y, max_x, max_y), + dtm_dir / f"{basename}_dtm{output_suffix}.tif") + except Exception as e: + logger.warning(f" Densité de points non écrite ({e})") + t_raster = time.perf_counter() dtm = bin_mean_2d(xs, ys, zs, width, height, (min_x, max_x), (min_y, max_y)) @@ -1861,24 +2062,18 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, # dense ce retour est la végétation, qui imprimerait les arbres dans # le MNT. - # Fill small gaps (< 1 m from data) precisely — comme avant + # Vides entre points comblés dans l'enveloppe des données (rayon + # adapté à la densité locale), îlots isolés retirés nan_count = np.count_nonzero(np.isnan(dtm)) if nan_count > 0: total = dtm.size nan_pct = 100.0 * nan_count / total logger.info(f" {nan_count:,} pixels sans données ({nan_pct:.1f}%)") - - max_gap_pixels = max(1, int(1.0 / resolution)) t_fill = time.perf_counter() - from rasterio.fill import fillnodata - valid_mask = ~np.isnan(dtm) - dtm_filled = fillnodata(dtm, mask=valid_mask, max_search_distance=max_gap_pixels) - small_gap_mask = np.isnan(dtm) & ~np.isnan(dtm_filled) - filled_count = np.count_nonzero(small_gap_mask) - if filled_count > 0: - dtm = np.where(small_gap_mask, dtm_filled, dtm) - logger.info(f" {filled_count:,} petits trous comblés " - f"(< {max_gap_pixels}px, {time.perf_counter() - t_fill:.1f}s)") + dtm, filled_count, removed_count = _fill_small_gaps(dtm, resolution) + logger.info(f" {filled_count:,} pixels comblés entre points, " + f"{removed_count:,} pixels d'îlots isolés retirés " + f"({time.perf_counter() - t_fill:.1f}s)") # Save as GeoTIFF output_tif = dtm_dir / f"{basename}_dtm{output_suffix}.tif" @@ -1894,6 +2089,9 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, compress='lzw' ) as dst: dst.write(dtm.astype('float32'), 1) + # Version du comblement : un DTM d'une version antérieure est + # régénéré (pastilles et liserés de l'ancien comblement). + dst.update_tags(**{GAP_FILL_TAG: str(GAP_FILL_VERSION)}) if used_edge_buffer > 0: # Tampon de raccord inscrit dans le fichier : changement de # --edge-buffer ⇒ invalidation automatique du cache DTM. diff --git a/lidar_pipeline/index.py b/lidar_pipeline/index.py index d638637..69932fa 100644 --- a/lidar_pipeline/index.py +++ b/lidar_pipeline/index.py @@ -6,7 +6,7 @@ Ce module produit aujourd'hui les artefacts dont le pipeline et la carte ont besoin : - registres partagés : VIZ_LABELS/VIZ_LEGENDS (libellés, légendes), défauts - d'affichage (DEFAULT_LAYERS/OPACITY/BLEND), correspondances mot-clé ↔ + d'affichage (DEFAULT_VIZ, PRECISION_VIZ, VIEW_MODES), correspondances mot-clé ↔ étape du pipeline (KEYWORD_TO_STEP) ; - vignettes (index_thumbs/) et sous-tuiles (index_subtiles/) servant de paliers sources à la pyramide XYZ (tiles.py) ; @@ -46,6 +46,7 @@ VIZ_LABELS = { 'solar': 'Éclairage solaire', 'anomaly': 'Carte d\'anomalies', 'relief_oriente': 'Relief orienté', + 'densite_sol': 'Précision (densité de points sol)', 'ortho': 'Orthophoto IGN', 'topo': 'Carte topographique IGN', } @@ -185,6 +186,15 @@ VIZ_LEGENDS = { 'gradient': None, 'ticks': None, }, + 'densite_sol': { + 'title': 'Précision géométrique (densité de points sol)', + 'legend': 'Points sol retenus pour le MNT, par m² (moyenne sur 3 × 3 m)\n16 gris, échelle log fixe : 2 niveaux = densité doublée\nNoir = ≤ 0,35 pt/m² ou aucun point (relief interpolé ou absent)\nBlanc = ≥ 45 pts/m²', + 'description': 'Où le relief est mesuré (clair) et où il est interpolé (sombre) — en « produit » sur le relief, assombrit les zones peu fiables', + 'cmap': 'gray', + 'gradient': ('#000000', '#202020', '#404040', '#606060', '#808080', + '#a0a0a0', '#c0c0c0', '#e0e0e0', '#ffffff'), + 'ticks': ('≤ 0,35 pt/m²', '≥ 45 pts/m²'), + }, 'ortho': { 'title': 'Photographie Aérienne IGN', 'legend': 'Orthophotographie\nImage aérienne', @@ -203,21 +213,14 @@ VIZ_LEGENDS = { }, } -# Couches activées par défaut à l'ouverture de la carte : seules celles-ci -# sont allumées, toutes les autres sont désactivées (mais restent accessibles -# via le panneau). Ordre = pile basse → haute : une couche opaque (opacité 1) -# doit rester en bas pour ne pas masquer les autres. Le relief orienté -# fusionne déjà openness et aspect en une seule image (plus d'empilement). -DEFAULT_LAYERS = ('relief_oriente',) - -# Opacité par défaut des couches allumées (0–1) ; les absentes restent à 1. -DEFAULT_OPACITY = {} - -# Mode de fusion (mix-blend-mode CSS) par couche : les couches d'ombrage en -# niveaux de gris fondent en « multiply » — elles assombrissent le relief sans -# délaver la rampe de couleurs en dessous, contrairement à la transparence. -# Les couches absentes restent en « normal » (simple transparence). -DEFAULT_BLEND = {'hillshade_multi': 'multiply', 'solar': 'multiply', 'svf': 'multiply'} +# Affichage de la carte : UNE couche principale (DEFAULT_VIZ) et la couche +# « précision » (densité de points sol), montrée seule ou superposée en +# « produit » sur la principale (assombrit les zones interpolées sans masquer +# le relief). Plus de pile de couches réordonnable. +PRECISION_VIZ = 'densite_sol' +VIEW_MODES = ('relief', 'precision', 'both') # principale / précision / les deux +DEFAULT_VIEW_MODE = 'relief' +DEFAULT_PRECISION_OPACITY = 0.6 # opacité de la précision en mode « les deux » # Couche « principale » (mise en avant ; index_tiles.json + affichage en # direct de la webapp) : le relief orienté, qui fusionne openness et aspect. @@ -228,7 +231,7 @@ DEFAULT_VIZ = 'relief_oriente' # carte (panneau, tuiles XYZ, TileJSON, WMTS). Les autres visualisations # restent calculables avec --only mais ne sont plus proposées. # None = toutes les visualisations présentes sur disque. -PANEL_VIZ = ('relief_oriente',) +PANEL_VIZ = ('relief_oriente', 'densite_sol') # Correspondance mot-clé de fichier de sortie → nom d'étape --only du pipeline # (les trois visualisations dont le nom de sortie diffère du nom d'étape, @@ -255,14 +258,13 @@ def step_to_keyword(step): return STEP_TO_KEYWORD.get(step, step) -def default_layers_present(all_viz_keys): - """Couches activées par défaut (DEFAULT_LAYERS) présentes sur le disque. - - Ordre conservé (DEFAULT_LAYERS) ; sert de base à l'état initial du - panneau (toutes les autres couches sont désactivées). - """ - present = set(all_viz_keys) - return [v for v in DEFAULT_LAYERS if v in present] +def default_main_layer(all_viz_keys): + """Couche principale par défaut présente sur le disque (DEFAULT_VIZ, sinon + la première qui n'est pas la précision), ou 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,)): @@ -289,11 +291,15 @@ def cells_with_all_viz(vis_dir, viz_keys, resolutions=(0.5,)): # les viz exclues retombent en repli dalle entière. _CARTO_SUBTILED_VIZ = () -# Encodage AVIF des sous-tuiles (benchmark sur dalles d'aspect réelles : -# q55 ≈ −45 % vs WebP q82, visuellement propre sur rampes de couleur). -# speed 0 (lent/meilleur) → 10 (rapide) ; 6+8 quasi identique en taille, -# on garde 8 (2× plus rapide) pour 2500×2500 px. -_SUBTILE_AVIF_QUALITY = 55 +# Encodage AVIF des sous-tuiles : q75 en 4:2:0, encodé UNE fois depuis le +# raster d'origine (write_subtiles appelé par tif_to_crop). Mesuré sur le +# relief orienté (2 dalles réelles, galerie de compression) : l'ancienne +# chaîne dalle q60 → sous-tuile q55 donnait 18,1 dB / SSIM 0,84 pour +# 3,9 Mo par dalle ; q75 unique 19,5 dB / SSIM 0,93 pour 8,8 Mo. Au-delà, +# le 4:2:0 plafonne (la teinte du relief, pixel par pixel, est moyennée par +# 2 × 2) : seul le 4:4:4 irait plus loin (q75 : 28 dB, 14 Mo). +# speed 9 : encodage rapide (cf. rendering.AVIF_SPEED). +_SUBTILE_AVIF_QUALITY = 75 _SUBTILE_AVIF_SPEED = 9 # Vignette de sous-tuile (px) : à la vue d'ensemble chaque sous-tuile active @@ -303,11 +309,15 @@ _SUBTILE_AVIF_SPEED = 9 _SUBTILE_THUMB_PX = 160 # Couches à contenu photographique ou traits fins (orthophoto, carte topo) : -# q55, calé sur des rampes de couleur, bave textes et linéaires — qualité -# relevée pour ces couches uniquement. +# qualité propre (égale aujourd'hui à celle des rampes de couleur). _SUBTILE_AVIF_QUALITY_DETAIL = 75 _SUBTILE_DETAIL_VIZ = frozenset({'ortho', 'topo'}) +# Couches en aplats de niveaux codés (cf. rendering.LOSSLESS_GRAY_KEYWORDS) : +# sous-tuiles en WebP sans perte (niveaux de gris ; 3× plus léger que l'AVIF +# q100, seul AVIF exact) et vignette intermédiaire sans perte. +_SUBTILE_LOSSLESS_VIZ = frozenset({'densite_sol'}) + # Vignette intermédiaire (px) : pallier entre la vignette 256 px et l'image # pleine résolution, pour éviter d'étirer la vignette ou de décoder l'AVIF # complet dès qu'une tuile dépasse ~300 px à l'écran. @@ -759,6 +769,74 @@ def _fallback_full_dalle(entries, viz_key, info): entry['viz'][viz_key] = dict(info) +def _subtile_ext(viz_key): + """Extension des sous-tuiles pleine résolution d'une couche.""" + return '.webp' if viz_key in _SUBTILE_LOSSLESS_VIZ else '.avif' + + +def _save_subtile_thumb(quad, path): + """Écrit la vignette de sous-tuile (redimensionnée à _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'): + """Découpe une image de dalle (PIL, nord en haut) en k × k sous-tuiles : + pleine résolution, vignette intermédiaire et vignette. + + Appelée par build_index depuis l'image de dalle, et par le pipeline + (rendering.tif_to_crop) directement depuis le raster d'origine : la + sous-tuile ne subit alors qu'un seul encodage avec perte (réencoder la + dalle AVIF cumulait deux pertes). Encodage : AVIF 4:2:0 q75 + (_SUBTILE_AVIF_QUALITY, q75 aussi pour ortho/topo), WebP sans perte pour + les aplats de niveaux (_SUBTILE_LOSSLESS_VIZ). Lève en cas d'échec. + """ + 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 : ligne 0 = nord → la sous-tuile j (nord) est en haut + 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) + # Anciens fichiers d'une génération précédente (autre format, + # vignette 256 px sans taille dans le nom) + (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): """Découpe une dalle en sous-tuiles (crops AVIF) pour la carte interactive. @@ -796,42 +874,32 @@ def _build_subtiles(tile, offered_viz_keys, output_dir, sub_dir_name): } try: - resample = getattr(PILImage, 'LANCZOS', 1) thumb_suffix = f"_thumb{_SUBTILE_THUMB_PX}.webp" - def _save_thumb(quad, stem): - """Écrit la vignette de sous-tuile (redimensionnée à _SUBTILE_THUMB_PX).""" - 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))), resample) - out.save(str(output_dir / sub_dir_name / (stem + thumb_suffix)), - format='WEBP', quality=80) - for viz_key in offered_viz_keys: info = tile['viz'].get(viz_key) if not info: continue - avif_quality = (_SUBTILE_AVIF_QUALITY_DETAIL - if viz_key in _SUBTILE_DETAIL_VIZ - else _SUBTILE_AVIF_QUALITY) + ext = _subtile_ext(viz_key) stems = {key: f"{tile['dir_name']}_{viz_key}_{key[0]}_{key[1]}" for key in entries} # Régénère si au moins un fichier manque ou est périmé (dalle # source recalculée depuis — comme les vignettes). L'URL pleine # porte un suffixe ?v= d'invalidation : on l'ôte pour le chemin. + # Des sous-tuiles écrites par le pipeline depuis le raster + # d'origine (tif_to_crop) sont plus récentes que la dalle : elles + # sont gardées telles quelles (un seul encodage avec perte). 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) - # AVIF + vignette intermédiaire d'un côté, vignettes de l'autre : - # un changement de taille de vignette (dalle source inchangée) - # re-découpe depuis les AVIF de sous-tuiles déjà présents, sans - # ré-encoder ces derniers (des minutes de CPU par rebuild). - heavy_fresh = all(_fresh(stem + '.avif') and _fresh(stem + '_mid.webp') + # Pleine résolution + vignette intermédiaire d'un côté, vignettes + # de l'autre : un changement de taille de vignette (dalle source + # inchangée) re-découpe depuis les sous-tuiles déjà présentes, sans + # ré-encoder ces dernières (des minutes de CPU par 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: @@ -839,11 +907,11 @@ def _build_subtiles(tile, offered_viz_keys, output_dir, sub_dir_name): f"{tile['dir_name']}/{viz_key} ({len(stems)} découpages)") try: for stem in stems.values(): - with PILImage.open(str(output_dir / sub_dir_name / (stem + '.avif'))) as q: + 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_thumb(q, stem) + _save_subtile_thumb(q, output_dir / sub_dir_name / (stem + thumb_suffix)) # Vignette 256 px d'une génération précédente (output_dir / sub_dir_name / (stem + '_thumb.webp')).unlink(missing_ok=True) except Exception as e: @@ -871,38 +939,11 @@ def _build_subtiles(tile, offered_viz_keys, output_dir, sub_dir_name): _fallback_full_dalle(entries, viz_key, info) if img is None: continue - if img.mode not in ('RGB', 'L'): - img = img.convert('RGB') - W, H = img.size - cut_failed = False - for (i, j), stem in stems.items(): - # Image : ligne 0 = nord → la sous-tuile j (nord) est en haut - 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)) - try: - quad.save(str(output_dir / sub_dir_name / (stem + '.avif')), - format='AVIF', quality=avif_quality, - speed=_SUBTILE_AVIF_SPEED) - except Exception as e: - logger.warning(f"Encodage AVIF impossible ({stem}), " - f"repli dalle entière pour {viz_key} : {e}") - cut_failed = True - break - # Supprime les anciens crops d'une génération précédente - # (repli .webp, vignette 256 px sans taille dans le nom) - (output_dir / sub_dir_name / (stem + '.webp')).unlink(missing_ok=True) - (output_dir / sub_dir_name / (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))), resample) - mid_img.save(str(output_dir / sub_dir_name / (stem + '_mid.webp')), - format='WEBP', quality=82) - _save_thumb(quad, stem) - if cut_failed: + try: + write_subtiles(output_dir, tile['dir_name'], viz_key, img, k, sub_dir_name) + except Exception as e: + logger.warning(f"Encodage des sous-tuiles impossible ({tile['dir_name']}), " + f"repli dalle entière pour {viz_key} : {e}") _fallback_full_dalle(entries, viz_key, info) continue # Chaque URL est versionnée par la mtime de SON fichier (cf. @@ -916,7 +957,7 @@ def _build_subtiles(tile, offered_viz_keys, output_dir, sub_dir_name): 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}.avif{_file_v(stem + '.avif')}", + 'full': f"{sub_dir_name}/{stem}{ext}{_file_v(stem + ext)}", } except Exception as e: logger.warning(f"Sous-tuilage abandonné pour {tile['dir_name']}: {e}") diff --git a/lidar_pipeline/mapserve.py b/lidar_pipeline/mapserve.py index 3c095c9..5672d3e 100644 --- a/lidar_pipeline/mapserve.py +++ b/lidar_pipeline/mapserve.py @@ -878,11 +878,12 @@ def tilejson(layer: str, request: Request): # --------------------------------------------------------------------------- -# Configuration par défaut de la pile de couches +# Configuration d'affichage par défaut # --------------------------------------------------------------------------- -# L'état réglé depuis la carte (ordre, visibilité, opacité, fusion) peut être -# figé comme configuration servie à TOUT nouveau navigateur : le réglage utile -# ne se perd pas dans un localStorage et se partage sans lien. +# L'affichage réglé depuis la carte (couche principale, mode relief / +# précision / les deux, opacité de la précision, fond) peut être figé comme +# configuration servie à TOUT nouveau navigateur : le réglage utile ne se perd +# pas dans un localStorage et se partage sans lien. DEFAULTS_FILE = OUTPUT_DIR / ".map-defaults.json" _defaults_lock = threading.Lock() @@ -903,38 +904,26 @@ def _clamp01(value, fallback=1.0): return fallback -BLEND_MODES = ("normal", "multiply", "screen", "overlay", "soft-light", - "hard-light", "difference", "luminosity") - - def _sanitize_defaults(req): - """Retient ce qui est connu et borné : couches réelles, opacités 0–1, fusions valides.""" - known = set(tiles_mod.available_layers(OUTPUT_DIR)) - order = [k for k in (req.order or []) if k in known] - for key in sorted(known - set(order)): - order.append(key) - on = [k for k in (req.on or []) if k in known] - opacity = {k: _clamp01(v) for k, v in (req.opacity or {}).items() if k in known} - blend = {k: v for k, v in (req.blend or {}).items() - if k in known and v in BLEND_MODES} + """Retient ce qui est connu et borné : couche réelle, mode valide, opacités 0–1.""" + from .index import DEFAULT_PRECISION_OPACITY, PRECISION_VIZ, VIEW_MODES + known = set(tiles_mod.available_layers(OUTPUT_DIR)) - {PRECISION_VIZ} base_in = req.base if isinstance(req.base, dict) else {} base = {"on": bool(base_in.get("on", True)), "opacity": _clamp01(base_in.get("opacity"), 0.85), "dark": bool(base_in.get("dark", True))} - stack = req.stack_blend if req.stack_blend in BLEND_MODES else "normal" - return {"order": order, "on": on, "opacity": opacity, "blend": blend, - "base": base, "stack_blend": stack, "saved_at": time.time()} + return {"main": req.main if req.main in known else None, + "mode": req.mode if req.mode in VIEW_MODES else "relief", + "precision_opacity": _clamp01(req.precision_opacity, DEFAULT_PRECISION_OPACITY), + "base": base, "saved_at": time.time()} class DefaultsRequest(BaseModel): - order: list = Field(default_factory=list, - description="ordre de pile, du bas vers le haut") - on: list = Field(default_factory=list, description="couches allumées") - opacity: dict = Field(default_factory=dict, description="opacité par couche, 0–1") - blend: dict = Field(default_factory=dict, description="mode de fusion par couche") + main: Optional[str] = Field(None, description="couche d'affichage principal") + mode: str = Field("relief", description="relief | precision | both") + precision_opacity: float = Field(0.6, description="opacité de la précision en mode both, 0–1") base: dict = Field(default_factory=dict, description="fond de carte : on, opacity, dark") - stack_blend: str = Field("normal", description="fusion de la pile sur le fond") @app.get("/api/map/defaults") @@ -952,8 +941,8 @@ def set_defaults(req: DefaultsRequest): DEFAULTS_FILE.parent.mkdir(parents=True, exist_ok=True) tmp.write_text(json.dumps(data, ensure_ascii=False), encoding="utf-8") os.replace(tmp, DEFAULTS_FILE) - logger.info(f"Configuration de pile par défaut enregistrée : " - f"{len(data['on'])} couche(s) allumée(s)") + logger.info(f"Configuration d'affichage par défaut enregistrée : " + f"{data['main'] or 'principale du registre'}, mode {data['mode']}") return {"enregistré": True, "defaults": data} @@ -973,28 +962,26 @@ def clear_defaults(): @app.get("/api/map/meta") def map_meta(): """Couches, zooms, emprise et version — tout l'état initial de l'interface.""" - from .index import DEFAULT_BLEND, DEFAULT_OPACITY, default_layers_present + from .index import (DEFAULT_PRECISION_OPACITY, DEFAULT_VIEW_MODE, PRECISION_VIZ, + VIEW_MODES, default_main_layer) infos = _layer_infos() keys = [i["key"] for i in infos] - saved = _load_defaults() + saved = _load_defaults() or {} # Configuration figée depuis la carte, sinon réglages du registre (index.py). - # Couches figées toutes disparues (ex. passage au seul relief orienté) : - # registre, sinon un nouveau navigateur arrive sur une carte sans LiDAR. - saved_on = [k for k in saved["on"] if k in keys] if saved else [] - if saved and (saved_on or not saved["on"]): - default_layers = saved_on - else: - default_layers = default_layers_present(keys) - default_opacity = saved["opacity"] if saved else dict(DEFAULT_OPACITY) - default_blend = saved["blend"] if saved else dict(DEFAULT_BLEND) + # Couche figée disparue, ou fichier de l'ancienne pile (clés order/on) : + # couche principale du registre. + main = saved.get("main") + if main not in keys or main == PRECISION_VIZ: + main = default_main_layer(keys) + mode = saved.get("mode") if saved.get("mode") in VIEW_MODES else DEFAULT_VIEW_MODE return { "layers": infos, - "default_layers": default_layers, - "default_opacity": default_opacity, - "default_blend": default_blend, - "default_order": saved["order"] if saved else None, - "default_base": saved["base"] if saved else None, - "default_stack_blend": saved["stack_blend"] if saved else "normal", + "default_main": main, + "precision_layer": PRECISION_VIZ if PRECISION_VIZ in keys else None, + "default_mode": mode, + "default_precision_opacity": _clamp01(saved.get("precision_opacity"), + DEFAULT_PRECISION_OPACITY), + "default_base": saved.get("base"), "defaults_saved": bool(saved), # L'interface consomme des tuiles 512 px (2× moins de requêtes qu'en # 256 : décisif en HTTP/1.1) ; les clients OSM gardent le 256 canonique. diff --git a/lidar_pipeline/mapui.py b/lidar_pipeline/mapui.py index 6bec269..bf24031 100644 --- a/lidar_pipeline/mapui.py +++ b/lidar_pipeline/mapui.py @@ -67,19 +67,34 @@ _MAP_HTML = """ -
Haut de pile — devant
-
-
Bas de pile — derrière
-
- Fusion de la pile sur le fond - +
+ Affichage + + +
+
- +
@@ -220,49 +235,36 @@ button:hover { background: rgba(255, 255, 255, 0.12); border-color: var(--accent .icon-btn:hover { color: var(--text); background: none; } .link-btn { flex: 1; font-size: 11.5px; } -/* Une couche tient sur deux lignes : identité (poignée, œil, nom, ordre) puis - réglages (fusion + opacité). Les boutons ▲▼ doublent le glisser-déposer — - seuls utilisables au doigt. */ +/* Ligne du fond de carte : œil, nom, opacité. */ .layer-row { display: grid; grid-template-columns: 14px 24px minmax(0, 1fr) auto 34px; align-items: center; gap: 6px; padding: 5px 4px; border-radius: 8px; user-select: none; } .layer-ctl { grid-column: 1 / -1; display: grid; grid-template-columns: minmax(0, 1fr) 84px; gap: 6px; align-items: center; padding: 2px 0 2px 38px; } -.layer-move { display: inline-flex; gap: 2px; } -.layer-move button { padding: 0 5px; font-size: 10px; line-height: 16px; color: var(--muted); - background: rgba(255, 255, 255, 0.05); } -.layer-move button:hover:not(:disabled) { color: var(--text); } -.layer-move button:disabled { opacity: 0.25; cursor: default; } -.blend { background: rgba(255, 255, 255, 0.06); color: var(--muted); - border: 1px solid var(--border); border-radius: 6px; padding: 2px 4px; - font: inherit; font-size: 11px; min-width: 0; width: 100%; height: 22px; } -.blend:hover { color: var(--text); } -.blend option { background-color: #10141b; color: var(--text); } .layer-row:hover { background: rgba(255, 255, 255, 0.05); } -/* Réorganisation : la destination doit être lisible AVANT de lâcher. Une - barre d'insertion épaisse marque le point de chute, la ligne déplacée - s'estompe, et à l'arrivée la couche clignote pour que l'œil la retrouve. */ -.layer-row.dragging { opacity: 0.35; outline: 1px dashed var(--accent); } -.layer-row.drop-before { box-shadow: inset 0 3px 0 0 var(--accent); } -.layer-row.drop-after { box-shadow: inset 0 -3px 0 0 var(--accent); } -.layer-row.drop-before::before, .layer-row.drop-after::before { - content: '▸'; position: absolute; margin-left: -12px; color: var(--accent); - font-size: 11px; line-height: 1; -} -.layer-row.drop-before::before { margin-top: -14px; } -.layer-row.drop-after::before { margin-top: 14px; } -.layer-row { position: relative; } -.layer-row.just-moved { animation: moved 1s ease-out; } -@keyframes moved { - 0% { background: rgba(90, 169, 230, 0.38); } - 100% { background: transparent; } -} -body.reordering, body.reordering * { cursor: grabbing !important; } -.stack-edge { - display: flex; align-items: center; gap: 6px; - font-size: 10px; letter-spacing: 0.04em; text-transform: uppercase; - color: var(--muted); padding: 2px 6px; -} -.stack-edge::after { content: ''; flex: 1; height: 1px; background: var(--border); } +/* Affichage : une couche principale, la précision seule ou superposée. */ +.view-row { display: grid; grid-template-columns: auto minmax(0, 1fr); align-items: center; + gap: 8px; padding: 8px 4px 4px; border-top: 1px solid var(--border); margin-top: 4px; } +.view-label { font-size: 10.5px; font-weight: 700; letter-spacing: 0.06em; text-transform: uppercase; color: var(--muted); } +.view-row .layer-name { color: var(--text); } +.view-sel { background: rgba(255, 255, 255, 0.06); color: var(--text); border: 1px solid var(--border); + border-radius: 6px; padding: 2px 4px; font: inherit; font-size: 12px; min-width: 0; width: 100%; height: 24px; } +.view-sel option { background-color: #10141b; color: var(--text); } +.mode-seg { display: grid; grid-template-columns: repeat(3, 1fr); margin: 6px 4px 2px; + border: 1px solid var(--border); border-radius: 8px; overflow: hidden; } +.mode-seg button { border: none; border-radius: 0; background: rgba(255, 255, 255, 0.04); + color: var(--muted); padding: 6px 4px; font-size: 12px; } +.mode-seg button + button { border-left: 1px solid var(--border); } +.mode-seg button:hover { color: var(--text); background: rgba(255, 255, 255, 0.1); } +.mode-seg button[aria-pressed="true"] { background: var(--accent); color: #0b0e13; font-weight: 600; } +.mode-seg button:focus-visible { outline: 2px solid var(--accent); outline-offset: -2px; } +.prec-ctl { display: grid; grid-template-columns: auto minmax(0, 1fr) 34px; gap: 8px; + align-items: center; padding: 6px 4px 0; font-size: 11.5px; color: var(--muted); } +.prec-legend { padding: 8px 4px 2px; } +.prec-bar { display: grid; grid-template-columns: repeat(16, 1fr); height: 10px; + border: 1px solid var(--border); border-radius: 3px; overflow: hidden; } +.prec-ticks { display: flex; justify-content: space-between; font-size: 10.5px; color: var(--muted); + font-variant-numeric: tabular-nums; margin-top: 2px; } +.prec-cap { font-size: 10.5px; color: var(--muted); margin-top: 2px; } .grip { color: var(--muted); cursor: grab; text-align: center; font-size: 12px; } .grip::before { content: '⠿'; } #baseRow .grip::before { content: ''; } @@ -371,35 +373,29 @@ _MAP_JS = r"""'use strict'; d'image, aucune rotation CSS, aucune mosaïque maison. */ const UI_VERSION = "__UI_VERSION__"; -const LS_KEY = 'lidarMapLayers_v1'; +const LS_KEY = 'lidarMapView_v2'; const el = (id) => document.getElementById(id); const SVG_EYE = ''; const SVG_EYE_OFF = ''; -// Modes de fusion proposés (CSS mix-blend-mode) : la superposition de couches -// ne se réduit pas à de la transparence — « Produit » assombrit le relief sans -// délaver la rampe du dessous, « Écran » éclaircit, « Superposé »/« Doux » -// creusent le contraste local, « Différence » fait ressortir les écarts. -const BLEND_MODES = [ - ['normal', 'Normal'], - ['multiply', 'Produit'], - ['screen', 'Écran'], - ['overlay', 'Superposé'], - ['soft-light', 'Doux'], - ['hard-light', 'Lumière crue'], - ['difference', 'Différence'], - ['luminosity', 'Luminosité'], -]; -const BLEND_KEYS = BLEND_MODES.map(m => m[0]); -const validBlend = (v) => (BLEND_KEYS.includes(v) ? v : null); -const blendOptions = (current) => BLEND_MODES.map( - m => '').join(''); +// Affichage : une couche principale (relief orienté) et la couche +// « précision » (densité de points sol), montrée seule ou superposée en +// « produit » (mix-blend-mode multiply) sur la principale : les zones peu +// fiables s'assombrissent sans masquer le relief. Touche P : mode suivant. +const VIEW_MODES = [['relief', 'Relief'], ['precision', 'Précision'], ['both', 'Les deux']]; +const VIEW_KEYS = VIEW_MODES.map(m => m[0]); +const validMode = (v) => (VIEW_KEYS.includes(v) ? v : null); +const clamp01 = (v, d) => { + const n = Number(v); + return Number.isFinite(n) ? Math.min(1, Math.max(0, n)) : d; +}; +// Paliers de la précision : niveau k à partir de 0,25 × 2^(k/2) pts/m² +// (visualizations.density_levels), gris = colormap « gray » sur 0..15. +const PREC_LEVELS = 16; let META = null; let STATE = null; const LAYERS = new Map(); // clé → L.TileLayer -let dragKey = null; const map = L.map('map', { zoomControl: false, @@ -506,33 +502,28 @@ function toast(msg, ms) { } // --- état persisté -------------------------------------------------------- +function precisionKey() { + const p = META && META.precision_layer; + return p && META.layers.some(l => l.key === p) ? p : null; +} + +function mainKeys(meta) { + return meta.layers.map(l => l.key).filter(k => k !== meta.precision_layer); +} + +function validMain(meta, k) { + return mainKeys(meta).includes(k) ? k : null; +} + function defaultState(meta) { - // `default_order` n'existe que si une configuration a été figée depuis la - // carte (POST /api/map/defaults) ; sinon on part du registre du pipeline. - const order = (Array.isArray(meta.default_order) && meta.default_order.length) - ? meta.default_order.filter(k => meta.layers.some(l => l.key === k)) - : meta.layers.map(l => l.key); - for (const l of meta.layers) if (!order.includes(l.key)) order.push(l.key); - const on = {}, opacity = {}, blend = {}; - const wanted = (meta.default_layers || []).filter(k => order.includes(k)); - // Ordre figé depuis la carte : respecté tel quel. Sinon, les couches - // allumées par défaut passent en bas de pile (bas → haut). - const stacked = (Array.isArray(meta.default_order) && meta.default_order.length) - ? order - : wanted.concat(order.filter(k => !wanted.includes(k))); - for (const k of stacked) { - on[k] = wanted.includes(k); - const o = Number((meta.default_opacity || {})[k]); - opacity[k] = Number.isFinite(o) ? Math.min(1, Math.max(0, o)) : 1; - blend[k] = (meta.default_blend || {})[k] || 'normal'; - } const b = meta.default_base || {}; return { - order: stacked, on, opacity, blend, + main: validMain(meta, meta.default_main) || mainKeys(meta)[0] || null, + mode: validMode(meta.default_mode) || 'relief', + precOpacity: clamp01(meta.default_precision_opacity, 0.6), base: { on: b.on !== undefined ? !!b.on : true, opacity: b.opacity !== undefined ? Number(b.opacity) : 0.85, dark: b.dark !== undefined ? !!b.dark : true }, - stackBlend: validBlend(meta.default_stack_blend) || 'normal', }; } @@ -543,33 +534,19 @@ function loadState(meta) { const raw = localStorage.getItem(LS_KEY); if (raw) { const s = JSON.parse(raw); - const order = (Array.isArray(s.order) ? s.order : []).filter(k => def.order.includes(k)); - for (const k of def.order) if (!order.includes(k)) order.push(k); - const on = {}, opacity = {}, blend = {}; - for (const k of order) { - on[k] = !!(s.on && s.on[k]); - const o = Number(s.opacity ? s.opacity[k] : NaN); - opacity[k] = Number.isFinite(o) ? Math.min(1, Math.max(0, o)) : def.opacity[k]; - blend[k] = (s.blend && s.blend[k]) || def.blend[k] || 'normal'; - } - st = { order, on, opacity, blend, - base: Object.assign({}, def.base, s.base || {}), - stackBlend: validBlend(s.stackBlend) || 'normal' }; + st = { main: validMain(meta, s.main) || def.main, + mode: validMode(s.mode) || def.mode, + precOpacity: clamp01(s.precOpacity, def.precOpacity), + base: Object.assign({}, def.base, s.base || {}) }; } } catch (e) { st = def; } + // Lien partagé : prime sur l'état local (les anciens liens &L=… de la + // pile de couches sont ignorés, la vue s'ouvre alors sur l'état courant). const hash = parseHash(); - if (hash && hash.layers) { - const order = hash.layers.order.filter(k => def.order.includes(k)); - for (const k of def.order) if (!order.includes(k)) order.push(k); - st = { order, on: hash.layers.on, opacity: hash.layers.opacity, - blend: Object.assign({}, st.blend, hash.layers.blend || {}), - base: hash.base || st.base, - stackBlend: hash.stackBlend || st.stackBlend || 'normal' }; - for (const k of order) { - if (st.on[k] === undefined) st.on[k] = false; - if (st.opacity[k] === undefined) st.opacity[k] = 1; - if (!st.blend[k]) st.blend[k] = def.blend[k] || 'normal'; - } + if (hash) { + if (hash.main && validMain(meta, hash.main)) st.main = hash.main; + if (hash.view) { st.mode = hash.view.mode; st.precOpacity = hash.view.opacity; } + if (hash.base) st.base = hash.base; } return st; } @@ -579,33 +556,23 @@ function saveState() { } // --- hash de vue ---------------------------------------------------------- -// #z/lat/lng&L=clé:opacité:visible,...&B=on:opacité:sombre +// #z/lat/lng&M=principale&P=mode:opacité&B=on:opacité:sombre function parseHash() { const m = /^#([\d.]+)\/(-?[\d.]+)\/(-?[\d.]+)(.*)$/.exec(location.hash || ''); if (!m) return null; const out = { zoom: parseFloat(m[1]), lat: parseFloat(m[2]), lng: parseFloat(m[3]), - layers: null, base: null }; + main: null, view: null, base: null }; const seg = {}; for (const part of m[4].split('&')) { const i = part.indexOf('='); if (i > 0) seg[part.slice(0, i)] = part.slice(i + 1); } - if (seg.L) { - // clé:opacité:visible[:fusion] — le 4e champ est optionnel (liens anciens) - const order = [], on = {}, opacity = {}, blend = {}; - for (const item of seg.L.split(',')) { - const f = item.split(':'); - if (f.length < 3) continue; - order.push(f[0]); - on[f[0]] = f[2] === '1'; - const o = Number(f[1]); - opacity[f[0]] = Number.isFinite(o) ? Math.min(100, Math.max(0, o)) / 100 : 1; - const b = validBlend(f[3]); - if (b) blend[f[0]] = b; - } - if (order.length) out.layers = { order, on, opacity, blend }; + if (seg.M) out.main = seg.M; + if (seg.P) { + const f = seg.P.split(':'); + const mode = validMode(f[0]); + if (mode) out.view = { mode, opacity: clamp01(Number(f[1]) / 100, 0.6) }; } - if (seg.S) out.stackBlend = validBlend(seg.S) || 'normal'; if (seg.B) { const f = seg.B.split(':'); if (f.length === 3) { @@ -620,12 +587,10 @@ function viewHash(full) { const c = map.getCenter(); let h = '#' + map.getZoom().toFixed(1) + '/' + c.lat.toFixed(6) + '/' + c.lng.toFixed(6); if (full) { - h += '&L=' + STATE.order.map(k => - k + ':' + Math.round(STATE.opacity[k] * 100) + ':' + (STATE.on[k] ? '1' : '0') - + ':' + (STATE.blend[k] || 'normal')).join(','); + if (STATE.main) h += '&M=' + STATE.main; + h += '&P=' + STATE.mode + ':' + Math.round(STATE.precOpacity * 100); h += '&B=' + (STATE.base.on ? '1' : '0') + ':' + Math.round(STATE.base.opacity * 100) + ':' + (STATE.base.dark ? '1' : '0'); - h += '&S=' + (STATE.stackBlend || 'normal'); } return h; } @@ -662,11 +627,10 @@ const EvenLevelTileLayer = ZoomAnimTileLayer.extend({ }); function buildLayers() { - // Conteneur isolé : `isolation: isolate` cantonne les mix-blend-mode à la - // pile LiDAR. Sans lui, « Produit » assombrirait aussi le fond OSM à travers - // les zones sans donnée. La fusion de la pile ENTIÈRE sur le fond reste - // réglable à part (STATE.stackBlend). Les panneaux sont réutilisés si déjà - // créés (rafraîchissement après un run de génération). + // Conteneur isolé : `isolation: isolate` cantonne le « produit » de la + // précision à la couche principale. Sans lui, il assombrirait aussi le fond + // OSM. Les panneaux sont réutilisés si déjà créés (rafraîchissement après + // un run de génération). const stack = map.getPane('lidarStack') || map.createPane('lidarStack'); stack.style.zIndex = 400; stack.style.isolation = 'isolate'; @@ -703,20 +667,30 @@ function buildLayers() { applyLayers(); } +// Mode réellement appliqué : sans couche précision servie, relief seul. +function effectiveMode() { + return precisionKey() ? STATE.mode : 'relief'; +} + function applyLayers() { - const stack = map.getPane('lidarStack'); - if (stack) stack.style.mixBlendMode = STATE.stackBlend || 'normal'; - let z = 400; - for (const key of STATE.order) { - const layer = LAYERS.get(key); - if (!layer) continue; + const prec = precisionKey(); + const mode = effectiveMode(); + for (const [key, layer] of LAYERS) { const pane = map.getPane('lidar-' + key); - z += 1; - pane.style.zIndex = String(z); - pane.style.mixBlendMode = STATE.blend[key] || 'normal'; - if (STATE.on[key]) { + let show, opacity = 1, blend = 'normal'; + if (key === prec) { + show = mode !== 'relief'; + if (mode === 'both') { opacity = STATE.precOpacity; blend = 'multiply'; } + } else { + show = key === STATE.main && mode !== 'precision'; + } + if (pane) { + pane.style.zIndex = key === prec ? '402' : '401'; + pane.style.mixBlendMode = blend; + } + if (show) { if (!map.hasLayer(layer)) layer.addTo(map); - layer.setOpacity(STATE.opacity[key]); + layer.setOpacity(opacity); } else if (map.hasLayer(layer)) { map.removeLayer(layer); } @@ -729,156 +703,60 @@ function applyLayers() { updateRose(); } -function moveLayer(key, delta) { - // STATE.order est bas → haut ; le panneau affiche haut en premier. - const ord = STATE.order; - const i = ord.indexOf(key); - const j = i + delta; - if (i < 0 || j < 0 || j >= ord.length) return; - ord.splice(j, 0, ord.splice(i, 1)[0]); - saveState(); renderPanel(key); applyLayers(); +function setMode(mode, announce) { + if (!validMode(mode)) return; + STATE.mode = mode; + saveState(); renderPanel(); applyLayers(); + if (announce) toast('Affichage : ' + VIEW_MODES.find(m => m[0] === mode)[1], 1200); } -// Insère `src` avant ou après `target` dans l'ordre AFFICHÉ (haut → bas). -function reorderRelative(src, target, before) { - const top = STATE.order.slice().reverse(); - const from = top.indexOf(src); - if (from < 0) return; - top.splice(from, 1); - let at = top.indexOf(target); - if (at < 0) return; - if (!before) at += 1; - top.splice(at, 0, src); - STATE.order = top.reverse(); - saveState(); renderPanel(src); applyLayers(); -} - -function clearDropMarks() { - document.querySelectorAll('.layer-row.drop-before, .layer-row.drop-after') - .forEach(r => r.classList.remove('drop-before', 'drop-after')); -} - -function renderPanel(movedKey) { - const host = el('layers'); - host.innerHTML = ''; - const labels = new Map(META.layers.map(l => [l.key, l])); - // Haut de pile en premier (STATE.order est bas → haut). - const topFirst = STATE.order.slice().reverse(); - topFirst.forEach((key, idx) => { - const info = labels.get(key); - if (!info) return; - const on = !!STATE.on[key]; - const row = document.createElement('div'); - row.className = 'layer-row' + (on ? ' on' : '') - + (key === movedKey ? ' just-moved' : ''); - row.dataset.layer = key; - row.title = info.description || info.label; - row.innerHTML = - '' + - '' + - '' + info.label + '' + - '' + - '' + - '' + - '' + - '' + Math.round(STATE.opacity[key] * 100) + '%' + - (on ? '
' + - '' + - '' + - '
' : ''); - row.draggable = false; - row.addEventListener('mousedown', (e) => { - row.draggable = !!(e.target.closest && e.target.closest('.grip')); - }); - row.querySelector('.eye').addEventListener('click', () => { - STATE.on[key] = !STATE.on[key]; - if (STATE.on[key] && STATE.opacity[key] === 0) STATE.opacity[key] = 1; - saveState(); renderPanel(); applyLayers(); - }); - row.querySelector('.up').addEventListener('click', () => moveLayer(key, 1)); - row.querySelector('.down').addEventListener('click', () => moveLayer(key, -1)); - const blend = row.querySelector('.blend'); - if (blend) blend.addEventListener('change', () => { - STATE.blend[key] = blend.value; - saveState(); - const pane = map.getPane('lidar-' + key); - if (pane) pane.style.mixBlendMode = blend.value; - }); - const range = row.querySelector('input[type="range"]'); - const pct = row.querySelector('.pct'); - if (range) { - range.addEventListener('input', () => { - STATE.opacity[key] = Number(range.value) / 100; - pct.textContent = range.value + '%'; - const layer = LAYERS.get(key); - if (layer && map.hasLayer(layer)) layer.setOpacity(STATE.opacity[key]); - }); - range.addEventListener('change', saveState); - } - row.addEventListener('dragstart', (e) => { - dragKey = key; - row.classList.add('dragging'); - document.body.classList.add('reordering'); - if (e.dataTransfer) { - e.dataTransfer.effectAllowed = 'move'; - e.dataTransfer.setData('text/plain', key); // exigé par Firefox - } - }); - row.addEventListener('dragend', () => { - dragKey = null; - row.classList.remove('dragging'); - row.draggable = false; - document.body.classList.remove('reordering'); - clearDropMarks(); - }); - row.addEventListener('dragover', (e) => { - e.preventDefault(); - if (!dragKey || dragKey === key) return; - // Moitié haute survolée → insertion au-dessus, sinon en dessous. - const box = row.getBoundingClientRect(); - const before = (e.clientY - box.top) < box.height / 2; - clearDropMarks(); - row.classList.add(before ? 'drop-before' : 'drop-after'); - }); - row.addEventListener('dragleave', () => { - row.classList.remove('drop-before', 'drop-after'); - }); - row.addEventListener('drop', (e) => { - e.preventDefault(); - if (!dragKey || dragKey === key) { clearDropMarks(); return; } - const before = row.classList.contains('drop-before'); - const src = dragKey; - clearDropMarks(); - document.body.classList.remove('reordering'); - reorderRelative(src, key, before); - }); - host.appendChild(row); - }); - // La couche qui vient d'être déplacée est ramenée dans le champ de vision - // (le panneau défile dès qu'il y a beaucoup de couches). - if (movedKey) { - const moved = host.querySelector('.layer-row.just-moved'); - if (moved && moved.scrollIntoView) moved.scrollIntoView({ block: 'nearest' }); +function renderPrecLegend() { + const bar = el('precBar'); + if (bar.childElementCount) return; + for (let k = 0; k < PREC_LEVELS; k++) { + const g = Math.round(255 * k / (PREC_LEVELS - 1)); + const lo = 0.25 * Math.pow(2, k / 2); + const hi = 0.25 * Math.pow(2, (k + 1) / 2); + const cell = document.createElement('i'); + cell.style.background = 'rgb(' + g + ',' + g + ',' + g + ')'; + cell.title = k === 0 ? '≤ ' + hi.toFixed(2).replace('.', ',') + ' pt/m²' + : k === PREC_LEVELS - 1 ? '≥ ' + Math.round(lo) + ' pts/m²' + : lo.toFixed(lo < 10 ? 1 : 0).replace('.', ',') + '–' + + hi.toFixed(hi < 10 ? 1 : 0).replace('.', ',') + ' pts/m²'; + bar.appendChild(cell); } } +function renderPanel() { + const labels = new Map(META.layers.map(l => [l.key, l])); + const mains = mainKeys(META); + const info = labels.get(STATE.main); + const name = el('mainName'), sel = el('mainSel'); + name.textContent = info ? info.label : 'Aucune couche'; + name.title = info ? (info.description || info.label) : ''; + // Plusieurs couches principales possibles : menu ; sinon simple libellé. + const many = mains.length > 1; + name.hidden = many; + sel.hidden = !many; + if (many) { + sel.innerHTML = mains.map(k => '').join(''); + } + const prec = precisionKey(); + el('precBlock').hidden = !prec; + if (!prec) return; + renderPrecLegend(); + const mode = effectiveMode(); + el('modeSeg').querySelectorAll('button').forEach(b => + b.setAttribute('aria-pressed', String(b.dataset.mode === mode))); + el('precCtl').hidden = mode !== 'both'; + el('precOpacity').value = Math.round(STATE.precOpacity * 100); + el('precPct').textContent = Math.round(STATE.precOpacity * 100) + '%'; + el('precLegend').hidden = mode === 'relief'; +} + function renderBaseRow() { const b = STATE.base; - const sel = el('stackBlend'); - if (sel && !sel.options.length) { - sel.innerHTML = blendOptions(STATE.stackBlend || 'normal'); - sel.addEventListener('change', () => { - STATE.stackBlend = sel.value; - saveState(); applyLayers(); - }); - } else if (sel) { - sel.value = STATE.stackBlend || 'normal'; - } const toggle = el('baseToggle'); toggle.classList.toggle('on', b.on); toggle.innerHTML = b.on ? SVG_EYE : SVG_EYE_OFF; @@ -1062,7 +940,8 @@ function drawRose() { } function updateRose() { - const on = !!(STATE && STATE.on && STATE.on.relief_oriente && LAYERS.get('relief_oriente')); + const on = !!(STATE && STATE.main === 'relief_oriente' && effectiveMode() !== 'precision' + && LAYERS.get('relief_oriente')); el('rose').hidden = !on; if (on) drawRose(); } @@ -1429,7 +1308,7 @@ window.addEventListener('hashchange', () => { if (!h) return; applyingHash = true; map.setView([h.lat, h.lng], h.zoom); - if (h.layers) { + if (h.main || h.view || h.base) { STATE = loadState(META); saveState(); renderPanel(); renderBaseRow(); applyLayers(); } @@ -1463,21 +1342,37 @@ el('btnShare').addEventListener('click', () => { }); el('btnDefault').addEventListener('click', () => { const body = { - order: STATE.order, - on: STATE.order.filter(k => STATE.on[k]), - opacity: Object.fromEntries(STATE.order.map(k => [k, STATE.opacity[k]])), - blend: Object.fromEntries(STATE.order.map(k => [k, STATE.blend[k] || 'normal'])), + main: STATE.main, + mode: STATE.mode, + precision_opacity: STATE.precOpacity, base: STATE.base, - stack_blend: STATE.stackBlend || 'normal', }; fetch('api/map/defaults', { method: 'POST', headers: { 'Content-Type': 'application/json' }, body: JSON.stringify(body), }).then(r => r.ok ? r.json() : Promise.reject('HTTP ' + r.status)) - .then(() => { META.default_order = STATE.order.slice(); toast('Configuration enregistrée comme défaut'); }) + .then(() => { + META.default_main = STATE.main; META.default_mode = STATE.mode; + META.default_precision_opacity = STATE.precOpacity; META.default_base = STATE.base; + toast('Affichage enregistré comme défaut'); + }) .catch(msg => toast('Échec de l\'enregistrement : ' + msg, 6000)); }); +el('modeSeg').querySelectorAll('button').forEach(b => + b.addEventListener('click', () => setMode(b.dataset.mode))); +el('mainSel').addEventListener('change', (e) => { + STATE.main = e.target.value; + saveState(); renderPanel(); applyLayers(); +}); +el('precOpacity').addEventListener('input', (e) => { + STATE.precOpacity = Number(e.target.value) / 100; + el('precPct').textContent = e.target.value + '%'; + const layer = LAYERS.get(precisionKey()); + if (layer && map.hasLayer(layer)) layer.setOpacity(STATE.precOpacity); +}); +el('precOpacity').addEventListener('change', saveState); + el('btnReset').addEventListener('click', () => { // Oublie l'état local et reprend la configuration servie par le serveur. try { localStorage.removeItem(LS_KEY); } catch (e) { /* privé */ } @@ -1529,6 +1424,14 @@ el('gps').addEventListener('click', () => { }, { enableHighAccuracy: true, timeout: 15000, maximumAge: 60000 }); }); document.addEventListener('keydown', (e) => { + // P : relief → précision → les deux → relief (hors saisie de texte) + const typing = e.target && e.target.closest && e.target.closest('input, select, textarea'); + if ((e.key === 'p' || e.key === 'P') && !typing && !e.ctrlKey && !e.metaKey && !e.altKey + && META && STATE && precisionKey()) { + const i = VIEW_KEYS.indexOf(STATE.mode); + setMode(VIEW_KEYS[(i + 1) % VIEW_KEYS.length], true); + return; + } if (e.key === 'Escape') { el('tilecard').hidden = true; el('usecard').hidden = true; if (genDrawing) genSetDrawing(false); diff --git a/lidar_pipeline/pipeline.py b/lidar_pipeline/pipeline.py index eaf0326..cd50907 100644 --- a/lidar_pipeline/pipeline.py +++ b/lidar_pipeline/pipeline.py @@ -90,10 +90,11 @@ from .visualizations import ( generate_flow_accumulation, generate_anomaly_mask, generate_relief_oriente, + generate_densite_sol, ) from .gpu import gpu_cleanup, num_gpus, available_gpu_ids, restrict_gpus, safe_gpu_call, gpu_worker_slots from .ign import generate_ign_overlay -from .rendering import tif_to_crop +from .rendering import tif_to_crop, LOSSLESS_GRAY_KEYWORDS # Ordered list of visualization steps. @@ -114,6 +115,7 @@ VIZ_STEPS = [ ('solar', generate_solar), ('anomaly', generate_anomaly_mask), ('relief_oriente', generate_relief_oriente), + ('densite_sol', generate_densite_sol), ('ortho', lambda d, b, v, r: generate_ign_overlay( d, b, v, r, layer='ORTHOIMAGERY.ORTHOPHOTOS', @@ -262,6 +264,8 @@ class LidarArchaeoPipeline: def _expected_output_path(name, basename, file_vis_dir, output_format='avif'): """Return the expected output filename for a visualization step.""" ext = 'avif' if output_format == 'avif' else 'webp' + if name in LOSSLESS_GRAY_KEYWORDS: + ext = 'webp' # aplats de niveaux : WebP sans perte (rendering.tif_to_crop) if name == 'pos_open': return file_vis_dir / f"{basename}_positive_openness.{ext}" elif name == 'neg_open': @@ -370,7 +374,7 @@ class LidarArchaeoPipeline: logger.info(f" Conversion images {fmt_label}:") for name, tif_file in vis_results.items(): if tif_file and isinstance(tif_file, Path) and tif_file.suffix == '.tif' and tif_file.exists(): - img_file = tif_to_crop(tif_file, file_vis_dir, resolution, keep_tif=self.keep_tif, quality=self.quality, output_format=self.output_format) + img_file = tif_to_crop(tif_file, file_vis_dir, resolution, keep_tif=self.keep_tif, quality=self.quality, output_format=self.output_format, subtiles_dir=self.output_dir) if img_file: logger.info(f" ✓ {img_file.name}") @@ -468,6 +472,12 @@ class LidarArchaeoPipeline: from .dtm import read_dtm_edge_buffer return abs(read_dtm_edge_buffer(dtm_path) - self.edge_buffer) < 1e-6 + def _gap_fill_matches(self, dtm_path): + """True si le DTM a été comblé par la version courante (tag + LIDAR_GAP_FILL, dtm.py) ; absent = ancien comblement à distance fixe.""" + from .dtm import read_dtm_gap_fill, GAP_FILL_VERSION + return read_dtm_gap_fill(dtm_path) == GAP_FILL_VERSION + def _fetch_edge_neighbors(self, files): """Télécharge les dalles LAZ voisines manquantes (raccord des bords). @@ -556,6 +566,10 @@ class LidarArchaeoPipeline: if not self._strip_align_matches(basename, res_suffix): logger.info(f" DTM{res_suffix} sans calage de faisceaux conforme — régénération (offsets verticaux mesurés et appliqués)") dtm_path.unlink() + elif not self._gap_fill_matches(dtm_path): + logger.info(f" DTM{res_suffix} au comblement d'une version antérieure — régénération " + f"(vides comblés dans l'enveloppe des points, rayon selon la densité)") + dtm_path.unlink() elif not self._edge_buffer_matches(dtm_path): from .dtm import read_dtm_edge_buffer recorded = read_dtm_edge_buffer(dtm_path) diff --git a/lidar_pipeline/rendering.py b/lidar_pipeline/rendering.py index 6beb02d..0ca6d54 100644 --- a/lidar_pipeline/rendering.py +++ b/lidar_pipeline/rendering.py @@ -160,8 +160,20 @@ COLORMAPS = { 'vmin_mode': 'fixed', 'vmin_val': 0, 'vmax_mode': 'fixed', 'vmax_val': 1, }, + # Densité de points sol : niveau 0..15 (visualizations.density_levels) + 'densite_sol': { + 'cmap': 'gray', + 'vmin_mode': 'fixed', 'vmin_val': 0, + 'vmax_mode': 'fixed', 'vmax_val': 15, + }, } +# Couches en aplats de gris codés (niveaux discrets) : image niveaux de gris +# (mode L) encodée SANS PERTE (AVIF q100 en L, seul réglage exact) — l'AVIF q55-60 déplace ~20 % des pixels jusqu'à 4 +# niveaux sur ces aplats. Elles restent à leur résolution propre (1 m pour +# la densité : ~150–260 Ko par dalle). +LOSSLESS_GRAY_KEYWORDS = ('densite_sol',) + # RGB entries (ortho/topo) are handled specially RGB_LEGENDS = { 'ortho': {}, @@ -824,7 +836,24 @@ def tif_to_png(tif_file, vis_dir, resolution, keep_tif=False, source_info=None, return None -def tif_to_crop(tif_file, vis_dir, resolution, keep_tif=False, quality=60, output_format='avif'): +def _write_subtiles_from(img, tif_file, vis_dir, resolution, subtiles_dir): + """Sous-tuiles d'une dalle depuis son image d'origine (best-effort : + en cas d'échec, build_index les découpera depuis la dalle).""" + from .index import VIZ_LABELS, _subdivision_k, write_subtiles + try: + k = _subdivision_k(resolution) + stem = Path(tif_file).stem + viz_key = next((v for v in sorted(VIZ_LABELS, key=len, reverse=True) + if stem.endswith(f"_{v}")), None) + if k <= 1 or viz_key is None: + return + write_subtiles(subtiles_dir, Path(vis_dir).name, viz_key, img, k) + except Exception as e: + logger.warning(f" Sous-tuiles non écrites ({Path(tif_file).name}) : {e}") + + +def tif_to_crop(tif_file, vis_dir, resolution, keep_tif=False, quality=60, output_format='avif', + subtiles_dir=None): """Convert GeoTIFF to a cropped visualization image (no legend, no overlay). Applies colormap and saves the image as a pure 1×1 km square. @@ -837,6 +866,9 @@ def tif_to_crop(tif_file, vis_dir, resolution, keep_tif=False, quality=60, outpu keep_tif: If True, keep the source TIFF after conversion. quality: Image quality (1-100). Use 100 for lossless. output_format: Output format ('webp' or 'avif'). + subtiles_dir: dossier de sortie du pipeline : les sous-tuiles de la + carte (index_subtiles) y sont écrites depuis l'image d'origine, + un seul encodage avec perte (build_index les trouve à jour). Returns: Path to output image file, or None on failure. @@ -844,7 +876,8 @@ def tif_to_crop(tif_file, vis_dir, resolution, keep_tif=False, quality=60, outpu if not tif_file or not tif_file.exists(): return None - ext = 'avif' if output_format == 'avif' else 'webp' + lossless_gray = any(k in tif_file.stem for k in LOSSLESS_GRAY_KEYWORDS) + ext = 'webp' if lossless_gray else ('avif' if output_format == 'avif' else 'webp') output_file = vis_dir / f"{tif_file.stem}.{ext}" try: @@ -892,12 +925,22 @@ def tif_to_crop(tif_file, vis_dir, resolution, keep_tif=False, quality=60, outpu # Save as AVIF/WebP img = PILImage.fromarray(rgb_data) pil_format = 'AVIF' if output_format == 'avif' else 'WEBP' - if quality >= 100: + if lossless_gray: + # WebP sans perte en niveaux de gris : exact et 3× plus léger que + # l'AVIF q100 (Pillow ignore `lossless=True` en AVIF : q75 avec perte) + img = PILImage.fromarray(rgb_data[:, :, 0]) + img.save(str(output_file), format='WEBP', lossless=True) + elif quality >= 100: img.save(str(output_file), format=pil_format, lossless=True) else: img.save(str(output_file), format=pil_format, quality=quality, **({'speed': AVIF_SPEED} if pil_format == 'AVIF' else {})) + # Sous-tuiles de la carte depuis l'image d'origine (après la dalle : + # plus récentes qu'elle, build_index ne les ré-encode pas) + if subtiles_dir is not None: + _write_subtiles_from(img, tif_file, vis_dir, resolution, subtiles_dir) + # Delete source TIFF (unless --keep-tif) if not keep_tif: tif_file.unlink(missing_ok=True) diff --git a/lidar_pipeline/tests/test_dtm.py b/lidar_pipeline/tests/test_dtm.py index 8d45fd2..60583a4 100644 --- a/lidar_pipeline/tests/test_dtm.py +++ b/lidar_pipeline/tests/test_dtm.py @@ -140,6 +140,105 @@ class TestInterpolateHoles: assert np.isnan(filled).all() +class TestFillSmallGaps: + """Comblement borné à l'enveloppe des points (plus de pastilles ni de liseré).""" + + @staticmethod + def _grid(n=200, step=2, value=10.0): + """Semis régulier de points (1 pixel sur `step`) sur n × n pixels.""" + dtm = np.full((n, n), np.nan) + dtm[::step, ::step] = value + return dtm + + def test_isolated_point_is_removed_not_grown(self): + from lidar_pipeline.dtm import _fill_small_gaps + dtm = np.full((100, 100), np.nan) + dtm[50, 50] = 5.0 + out, filled, removed = _fill_small_gaps(dtm, 0.2) + assert np.isnan(out).all() # ni pastille, ni point seul + assert (filled, removed) == (0, 1) + + def test_gaps_between_points_filled_without_edge_band(self): + from lidar_pipeline.dtm import _fill_small_gaps + dtm = self._grid() + dtm[:, 100:] = np.nan # grand trou à l'est + out, filled, _ = _fill_small_gaps(dtm, 0.2) + assert filled > 0 + assert not np.isnan(out[10:190, 10:99]).any() # vides entre points comblés + assert np.isnan(out[:, 99:]).all() # rien d'extrapolé dans le trou + + def test_radius_follows_local_density(self): + """Semis clairsemé (1 point / 1,8 m, vides de 2,5 m en diagonale) + comblé ; trou de 3 m en zone dense conservé ; trou de 1,6 m comblé.""" + from lidar_pipeline.dtm import _fill_small_gaps + sparse = self._grid(step=9) + out, _, _ = _fill_small_gaps(sparse, 0.2) + assert not np.isnan(out[30:170, 30:170]).any() + dense = self._grid(step=1) + dense[90:105, 90:105] = np.nan # trou de 3 m dans un semis plein + dense[40:48, 40:48] = np.nan # trou de 1,6 m (voiture) + out, _, _ = _fill_small_gaps(dense, 0.2) + assert np.isnan(out[95:100, 95:100]).all() + assert not np.isnan(out[40:48, 40:48]).any() + + def test_morph_disk_matches_scipy(self): + from scipy import ndimage as nd + from lidar_pipeline.dtm import _morph_disk + rng = np.random.default_rng(0) + mask = rng.random((60, 70)) < 0.05 + square = np.ones((3, 3), dtype=bool) + cross = nd.generate_binary_structure(2, 1) + ref = mask + for i in range(5): + ref = nd.binary_dilation(ref, structure=square if i % 2 == 0 else cross) + assert (_morph_disk(mask, 5) == ref).all() + ref_e = ref + for i in range(5): + ref_e = nd.binary_erosion(ref_e, structure=square if i % 2 == 0 else cross, + border_value=1) + assert (_morph_disk(ref, 5, erode=True) == ref_e).all() + + def test_dtm_records_gap_fill_version(self, tmp_output_dir): + import laspy + import rasterio + from lidar_pipeline.dtm import (create_dtm_fast, read_dtm_gap_fill, + GAP_FILL_VERSION) + hdr = laspy.LasHeader(version='1.2', point_format=0) + las = laspy.LasData(hdr) + las.x, las.y, las.z = [0.5, 1.5, 0.5, 1.5], [0.5, 0.5, 1.5, 1.5], [10.0] * 4 + las.write(str(tmp_output_dir / "g.las")) + out = create_dtm_fast(tmp_output_dir / "g.las", "g", tmp_output_dir, 1.0, + force=True, strip_align=False) + assert read_dtm_gap_fill(out) == GAP_FILL_VERSION + with rasterio.open(str(out)) as src: + src.tags() # lisible + legacy = tmp_output_dir / "legacy.tif" + with rasterio.open(str(legacy), "w", driver="GTiff", width=2, height=2, + count=1, dtype="float32") as dst: + dst.write(np.zeros((1, 2, 2), dtype="float32")) + assert read_dtm_gap_fill(legacy) == 1 + + +class TestDensitySidecar: + def test_dtm_writes_ground_density(self, tmp_output_dir): + """4 points par m² sur 10 × 10 m : densité 4 au cœur, grille de 1 m.""" + import laspy + import rasterio + from lidar_pipeline.dtm import create_dtm_fast, density_path + g = (np.arange(20) + 0.25) / 2.0 # pas de 0,5 m + xx, yy = np.meshgrid(g, g) + hdr = laspy.LasHeader(version='1.2', point_format=0) + las = laspy.LasData(hdr) + las.x, las.y, las.z = xx.ravel(), yy.ravel(), np.full(xx.size, 10.0) + las.write(str(tmp_output_dir / "d.las")) + out = create_dtm_fast(tmp_output_dir / "d.las", "d", tmp_output_dir, 0.5, + force=True, strip_align=False) + with rasterio.open(density_path(out)) as src: + dens = src.read(1) + assert abs(src.transform.a - 1.0) < 1e-9 + assert np.allclose(dens[2:-2, 2:-2], 4.0) + + class TestDetectGroundMethod: def _make_mock_las(self, num_returns, z_values): """Create a mock laspy object with specified NumberOfReturns and z.""" diff --git a/lidar_pipeline/tests/test_index.py b/lidar_pipeline/tests/test_index.py index 8151f1d..9e05d61 100644 --- a/lidar_pipeline/tests/test_index.py +++ b/lidar_pipeline/tests/test_index.py @@ -287,17 +287,15 @@ def test_subtiles_cover_all_layers(tmp_path): def test_panel_restricted_to_kept_layers(): - """Le panneau est restreint aux couches conservées (PANEL_VIZ). - - Les couches allumées par défaut (DEFAULT_LAYERS) y figurent toutes — le - panneau peut en proposer davantage (ex. pente éteinte par défaut), et - chaque opacité par défaut vise une couche allumée. - """ - from lidar_pipeline.index import (PANEL_VIZ, DEFAULT_LAYERS, DEFAULT_OPACITY, - KEYWORD_TO_STEP) + """Le panneau est restreint aux couches conservées (PANEL_VIZ) : la couche + principale par défaut et la précision y figurent.""" + from lidar_pipeline.index import (PANEL_VIZ, DEFAULT_VIZ, PRECISION_VIZ, VIEW_MODES, + DEFAULT_VIEW_MODE, KEYWORD_TO_STEP, default_main_layer) from lidar_pipeline.pipeline import VIZ_STEPS - assert set(DEFAULT_LAYERS) <= set(PANEL_VIZ) - assert set(DEFAULT_OPACITY) <= set(DEFAULT_LAYERS) + assert DEFAULT_VIZ in PANEL_VIZ and PRECISION_VIZ in PANEL_VIZ + assert DEFAULT_VIEW_MODE in VIEW_MODES + assert default_main_layer(["densite_sol", "aspect"]) == "aspect" + assert default_main_layer(["densite_sol"]) is None # Chaque couche conservée correspond à une étape --only valide steps = {name for name, _ in VIZ_STEPS} for key in PANEL_VIZ: diff --git a/lidar_pipeline/tests/test_mapserve.py b/lidar_pipeline/tests/test_mapserve.py index b074a6e..841a0e1 100644 --- a/lidar_pipeline/tests/test_mapserve.py +++ b/lidar_pipeline/tests/test_mapserve.py @@ -210,7 +210,9 @@ def test_map_meta(tmp_path, monkeypatch): assert meta["tile_size"] == 512 and meta["zoom_offset"] == -1 assert meta["max_native_zoom"] == 19 assert meta["bounds"] and meta["stamp"] > 0 - assert meta["default_layers"] == ["relief_oriente"] + assert meta["default_main"] == "relief_oriente" + assert meta["precision_layer"] is None # densité absente de la fixture + assert meta["default_mode"] == "relief" def test_only_panel_layers_are_served(tmp_path, monkeypatch): @@ -256,31 +258,31 @@ def test_healthz_and_assets(tmp_path, monkeypatch): mapserve.assets("introuvable.js") -def test_ui_offers_reorder_and_blend_modes(): - """Pile de couches : réordonnancement au doigt ET modes de fusion. +def test_ui_main_layer_and_precision_modes(): + """Une couche d'affichage principal et la précision : relief seul, + précision seule, ou précision en « produit » sur le relief. - La superposition ne se réduit pas à de la transparence : chaque couche - porte un mix-blend-mode, cantonné à la pile LiDAR par un conteneur isolé - (sinon « Produit » assombrirait aussi le fond de carte). + Plus de pile réordonnable ni de fusion par couche. Le produit est + cantonné à la pile LiDAR par un conteneur isolé (sinon il assombrirait + aussi le fond de carte). Bascule rapide au clavier (touche P). """ from lidar_pipeline.mapui import _MAP_CSS, _MAP_HTML, _MAP_JS - # Modes proposés - for mode in ("multiply", "screen", "overlay", "soft-light", "difference", - "luminosity"): - assert f"'{mode}'" in _MAP_JS, mode - # Réordonnancement : glisser-déposer ET boutons (seuls utilisables au doigt) - assert "function moveLayer(" in _MAP_JS - assert "moveLayer(key, 1)" in _MAP_JS and "moveLayer(key, -1)" in _MAP_JS - assert "dragstart" in _MAP_JS and "layer-move" in _MAP_CSS - # Pile isolée + fusion par couche et fusion de la pile sur le fond - assert "createPane('lidarStack')" in _MAP_JS - assert "isolation = 'isolate'" in _MAP_JS - assert "map.createPane(paneName, stack)" in _MAP_JS - assert "pane.style.mixBlendMode" in _MAP_JS - assert "STATE.stackBlend" in _MAP_JS and 'id="stackBlend"' in _MAP_HTML - # Le lien partagé transporte la fusion (4e champ, rétrocompatible) - assert "(STATE.blend[k] || 'normal')" in _MAP_JS - assert "const b = validBlend(f[3]);" in _MAP_JS + for mode in ('data-mode="relief"', 'data-mode="precision"', 'data-mode="both"'): + assert mode in _MAP_HTML + assert 'id="precOpacity"' in _MAP_HTML and 'id="precLegend"' in _MAP_HTML + assert 'id="mainSel"' in _MAP_HTML + assert "blend = 'multiply'" in _MAP_JS + assert "createPane('lidarStack')" in _MAP_JS and "isolation = 'isolate'" in _MAP_JS + assert "e.key === 'p'" in _MAP_JS + # Légende en 16 paliers exacts (échelle log de la densité) + assert "const PREC_LEVELS = 16;" in _MAP_JS and "repeat(16, 1fr)" in _MAP_CSS + # Pile et fusions supprimées + for gone in ("function moveLayer(", "reorderRelative", "stackBlend", "dragstart", "blendOptions"): + assert gone not in _MAP_JS, gone + assert "Haut de pile" not in _MAP_HTML + # Lien partagé : principale, mode:opacité, fond + assert "'&M=' + STATE.main" in _MAP_JS + assert "'&P=' + STATE.mode" in _MAP_JS def test_map_ui_tile_load_indicator_and_selection_pane(): @@ -301,85 +303,73 @@ def test_map_ui_tile_load_indicator_and_selection_pane(): def test_defaults_roundtrip(tmp_path, monkeypatch): - """L'état réglé depuis la carte devient la configuration servie à tous.""" + """L'affichage réglé depuis la carte devient la configuration servie à tous.""" import lidar_pipeline.mapserve as mapserve - _setup(tmp_path, monkeypatch) + _setup(tmp_path, monkeypatch, layers=("relief_oriente", "aspect", "densite_sol")) monkeypatch.setattr(mapserve, "DEFAULTS_FILE", tmp_path / ".map-defaults.json") # Sans enregistrement : réglages du registre (index.py) meta = mapserve.map_meta() assert meta["defaults_saved"] is False - assert meta["default_order"] is None - assert meta["default_stack_blend"] == "normal" + assert meta["default_main"] == "relief_oriente" + assert meta["precision_layer"] == "densite_sol" + assert meta["default_mode"] == "relief" + assert meta["default_precision_opacity"] == 0.6 - req = mapserve.DefaultsRequest( - order=["slope", "aspect"], on=["aspect"], - opacity={"aspect": 0.6, "slope": 1.0}, - blend={"aspect": "multiply", "slope": "normal"}, - base={"on": True, "opacity": 0.4, "dark": False}, - stack_blend="soft-light") + req = mapserve.DefaultsRequest(main="aspect", mode="both", precision_opacity=0.4, + base={"on": True, "opacity": 0.4, "dark": False}) assert mapserve.set_defaults(req)["enregistré"] is True meta = mapserve.map_meta() assert meta["defaults_saved"] is True - assert meta["default_order"] == ["slope", "aspect"] - assert meta["default_layers"] == ["aspect"] - assert meta["default_opacity"]["aspect"] == 0.6 - assert meta["default_blend"]["aspect"] == "multiply" + assert meta["default_main"] == "aspect" + assert meta["default_mode"] == "both" + assert meta["default_precision_opacity"] == 0.4 assert meta["default_base"] == {"on": True, "opacity": 0.4, "dark": False} - assert meta["default_stack_blend"] == "soft-light" # Persisté sur disque : survit au redémarrage du conteneur assert (tmp_path / ".map-defaults.json").is_file() - assert mapserve.get_defaults()["defaults"]["on"] == ["aspect"] + assert mapserve.get_defaults()["defaults"]["mode"] == "both" # Retrait : retour aux réglages du registre assert mapserve.clear_defaults()["supprimé"] is True assert mapserve.map_meta()["defaults_saved"] is False -def test_defaults_obsolete_layers_fall_back(tmp_path, monkeypatch): - """Défauts figés sur des couches qui ne sont plus servies : registre. - - Sinon un nouveau navigateur arrive sur une carte sans aucune couche - LiDAR allumée (cas de la prod après le passage au seul relief orienté). - """ +def test_defaults_legacy_stack_file_falls_back(tmp_path, monkeypatch): + """Fichier de l'ancienne pile (order/on/blend) ou couche disparue : + couche principale du registre, mode relief ; le fond est conservé.""" import json import lidar_pipeline.mapserve as mapserve _setup(tmp_path, monkeypatch, layers=("relief_oriente",)) defaults = tmp_path / ".map-defaults.json" monkeypatch.setattr(mapserve, "DEFAULTS_FILE", defaults) defaults.write_text(json.dumps({ - "order": ["aspect", "positive_openness"], - "on": ["aspect", "positive_openness"], "opacity": {}, "blend": {}, - "base": {"on": True, "opacity": 0.85, "dark": True}, - "stack_blend": "normal"}), encoding="utf-8") + "order": ["aspect"], "on": ["aspect"], "opacity": {}, "blend": {}, + "base": {"on": True, "opacity": 0.5, "dark": True}, "stack_blend": "normal"}), + encoding="utf-8") meta = mapserve.map_meta() - assert meta["default_layers"] == ["relief_oriente"] - # Tout éteint volontairement (couches connues) : respecté - defaults.write_text(json.dumps({ - "order": ["relief_oriente"], "on": [], "opacity": {}, "blend": {}, - "base": {}, "stack_blend": "normal"}), encoding="utf-8") - assert mapserve.map_meta()["default_layers"] == [] + assert meta["default_main"] == "relief_oriente" + assert meta["default_mode"] == "relief" + assert meta["default_base"]["opacity"] == 0.5 + defaults.write_text(json.dumps({"main": "aspect", "mode": "vaudou"}), encoding="utf-8") + meta = mapserve.map_meta() + assert meta["default_main"] == "relief_oriente" and meta["default_mode"] == "relief" def test_defaults_sanitised(tmp_path, monkeypatch): - """Couches inconnues écartées, opacités bornées, fusions validées.""" + """Couche inconnue ou précision refusée comme principale, mode validé, opacités bornées.""" import lidar_pipeline.mapserve as mapserve - _setup(tmp_path, monkeypatch) + _setup(tmp_path, monkeypatch, layers=("relief_oriente", "densite_sol")) monkeypatch.setattr(mapserve, "DEFAULTS_FILE", tmp_path / ".map-defaults.json") data = mapserve.set_defaults(mapserve.DefaultsRequest( - order=["aspect", "inexistante"], on=["aspect", "inexistante"], - opacity={"aspect": 5, "slope": "x", "inexistante": 0.5}, - blend={"aspect": "vaudou", "slope": "screen"}, - base={"opacity": -3}, stack_blend="vaudou"))["defaults"] - assert "inexistante" not in data["order"] and "inexistante" not in data["on"] - # Les couches présentes mais non citées complètent l'ordre - assert set(data["order"]) == {"aspect", "slope"} - assert data["opacity"]["aspect"] == 1.0 # borné à 1 - assert data["opacity"]["slope"] == 1.0 # valeur illisible → opaque - assert data["blend"] == {"slope": "screen"} # mode inconnu écarté + main="inexistante", mode="vaudou", precision_opacity=5, + base={"opacity": -3}))["defaults"] + assert data["main"] is None + assert data["mode"] == "relief" + assert data["precision_opacity"] == 1.0 # borné à 1 assert data["base"]["opacity"] == 0.0 # borné à 0 - assert data["stack_blend"] == "normal" # repli + data = mapserve.set_defaults(mapserve.DefaultsRequest(main="densite_sol"))["defaults"] + assert data["main"] is None # la précision n'est pas une principale def test_ui_applies_server_defaults(): @@ -387,38 +377,13 @@ def test_ui_applies_server_defaults(): from lidar_pipeline.mapui import _MAP_HTML, _MAP_JS assert 'id="btnDefault"' in _MAP_HTML and 'id="btnReset"' in _MAP_HTML assert "api/map/defaults" in _MAP_JS - assert "meta.default_order" in _MAP_JS + assert "meta.default_main" in _MAP_JS and "meta.default_mode" in _MAP_JS assert "meta.default_base" in _MAP_JS - assert "validBlend(meta.default_stack_blend)" in _MAP_JS + assert "meta.default_precision_opacity" in _MAP_JS # Réinitialiser oublie l'état local avant de reprendre celui du serveur assert "localStorage.removeItem(LS_KEY)" in _MAP_JS -def test_ui_reorder_shows_destination(): - """Réorganiser doit être lisible : repère de chute, fantôme, arrivée signalée. - - Sans repère, on lâche une couche sans savoir où elle atterrit — la pile - n'est pas un détail, elle décide de ce qu'on voit. - """ - from lidar_pipeline.mapui import _MAP_CSS, _MAP_HTML, _MAP_JS - # Barre d'insertion selon la moitié survolée, effacée entre deux survols - assert ".layer-row.drop-before" in _MAP_CSS and ".layer-row.drop-after" in _MAP_CSS - assert "clearDropMarks" in _MAP_JS - assert "box.height / 2" in _MAP_JS - assert "row.classList.add(before ? 'drop-before' : 'drop-after')" in _MAP_JS - # Ligne déplacée estompée pendant le glisser - assert ".layer-row.dragging" in _MAP_CSS and "classList.add('dragging')" in _MAP_JS - # Insertion relative (avant/après la ligne visée), pas un simple échange - assert "function reorderRelative(" in _MAP_JS - # Arrivée mise en évidence puis ramenée dans le champ de vision - assert "just-moved" in _MAP_CSS and "@keyframes moved" in _MAP_CSS - assert "scrollIntoView" in _MAP_JS - assert "renderPanel(key)" in _MAP_JS and "renderPanel(src)" in _MAP_JS - # Extrémités de la pile nommées - assert 'id="edgeTop"' in _MAP_HTML and 'id="edgeBottom"' in _MAP_HTML - assert "Haut de pile" in _MAP_HTML and "Bas de pile" in _MAP_HTML - - def test_ui_uses_standard_tilelayer(): """L'interface s'appuie sur le LOD natif de Leaflet, pas sur un palier maison.""" from lidar_pipeline.mapui import _MAP_JS, render_html, ui_version diff --git a/lidar_pipeline/tests/test_pipeline.py b/lidar_pipeline/tests/test_pipeline.py index d0ca77b..21e9d25 100644 --- a/lidar_pipeline/tests/test_pipeline.py +++ b/lidar_pipeline/tests/test_pipeline.py @@ -20,16 +20,16 @@ class TestVizSteps: assert len(names) == len(set(names)), "VIZ_STEPS has duplicate names" def test_expected_visualization_count(self): - """Should have 16 visualizations (14 terrain + ortho + topo).""" + """17 visualisations : 14 terrain + densité de points + ortho + topo.""" from lidar_pipeline.pipeline import VIZ_STEPS - assert len(VIZ_STEPS) == 16 + assert len(VIZ_STEPS) == 17 - def test_default_run_produces_only_relief(self, tmp_path): - """Sans --only : seule la couche affichée (relief orienté) est produite ; - --only reste libre pour les autres visualisations.""" + def test_default_run_produces_only_panel_layers(self, tmp_path): + """Sans --only : seules les couches affichées (relief orienté, densité + de points) sont produites ; --only reste libre pour les autres.""" from lidar_pipeline.pipeline import LidarArchaeoPipeline p = LidarArchaeoPipeline(tmp_path, tmp_path / "out") - assert [n for n, _ in p.viz_steps] == ["relief_oriente"] + assert [n for n, _ in p.viz_steps] == ["relief_oriente", "densite_sol"] p = LidarArchaeoPipeline(tmp_path, tmp_path / "out2", only_viz=["slope"]) assert [n for n, _ in p.viz_steps] == ["slope"] diff --git a/lidar_pipeline/tests/test_rendering.py b/lidar_pipeline/tests/test_rendering.py index 0993fcc..41b0886 100644 --- a/lidar_pipeline/tests/test_rendering.py +++ b/lidar_pipeline/tests/test_rendering.py @@ -217,3 +217,52 @@ class TestCoreTileWindow: (637900, 6626900, 639100, 6628100), 1200) with rasterio.open(tif) as src: assert _core_tile_window(tif, src) is None + + +class TestDensiteSolCrop: + """Densité de points : WebP sans perte en niveaux de gris, sous-tuiles + écrites depuis l'image d'origine (un seul encodage).""" + + def test_lossless_gray_levels_and_subtiles(self, tmp_path): + from PIL import Image as PILImage + from lidar_pipeline.rendering import tif_to_crop + base = "LHD_FXX_0660_6701_PTS_LAMB93_IGN69" + vis = tmp_path / "visualisations" / f"{base}_r0p2" + vis.mkdir(parents=True) + levels = np.tile(np.arange(16, dtype=np.float32).repeat(4), (64, 1)) # 64 × 64 + tif = TestTifToCrop._write_named_tif(vis, f"{base}_densite_sol.tif", levels) + out = tif_to_crop(tif, vis, 0.2, output_format='avif', subtiles_dir=tmp_path) + assert out is not None and out.suffix == ".webp" + # WebP n'a pas de mode gris : relu en RGB à canaux égaux + rgb = np.asarray(PILImage.open(str(out)).convert("RGB")).astype(int) + assert (rgb[..., 0] == rgb[..., 1]).all() and (rgb[..., 1] == rgb[..., 2]).all() + img = PILImage.fromarray(rgb[..., 0].astype(np.uint8)) + row = np.asarray(img)[0, ::4].astype(int) + assert len(set(row.tolist())) == 16 and np.all(np.diff(row) > 0) + # Sous-tuiles 2 × 2 (0,2 m/px) en WebP sans perte, identiques à la dalle + q = tmp_path / "index_subtiles" / f"{base}_r0p2_densite_sol_0_1.webp" + assert q.exists() + assert np.array_equal(np.asarray(PILImage.open(str(q)).convert("L")), np.asarray(img)[:32, :32]) + assert (tmp_path / "index_subtiles" / f"{base}_r0p2_densite_sol_0_1_mid.webp").exists() + + def test_relief_subtiles_encoded_from_source(self, tmp_path): + """Relief : sous-tuiles AVIF écrites par tif_to_crop, plus récentes que la dalle.""" + import pytest + from lidar_pipeline.rendering import tif_to_crop + base = "LHD_FXX_0660_6701_PTS_LAMB93_IGN69" + vis = tmp_path / "visualisations" / f"{base}_r0p2" + vis.mkdir(parents=True) + rgb = np.random.default_rng(1).integers(0, 255, (3, 64, 64)).astype('uint8') + tif = vis / f"{base}_relief_oriente.tif" + with rasterio.open(tif, 'w', driver='GTiff', height=64, width=64, count=3, + dtype='uint8', crs='EPSG:2154', + transform=from_bounds(660000, 6700000, 661000, 6701000, 64, 64)) as dst: + dst.write(rgb) + try: + out = tif_to_crop(tif, vis, 0.2, output_format='avif', subtiles_dir=tmp_path) + except Exception: + pytest.skip("encodeur AVIF indisponible") + sub = tmp_path / "index_subtiles" + quads = sorted(sub.glob(f"{base}_r0p2_relief_oriente_?_?.avif")) + assert len(quads) == 4 + assert all(q.stat().st_mtime_ns >= out.stat().st_mtime_ns for q in quads) diff --git a/lidar_pipeline/tests/test_tiles.py b/lidar_pipeline/tests/test_tiles.py index afaf62f..69675b0 100644 --- a/lidar_pipeline/tests/test_tiles.py +++ b/lidar_pipeline/tests/test_tiles.py @@ -687,3 +687,28 @@ def test_existing_cache_without_registry_is_refreshed_once(tmp_path): tiles._seen_cache.clear() tiles.source_index(tmp_path, force=True) assert tiles.cached_tile(tmp_path, "slope", z, x, y)[1] == "fresh" + + +def test_webp_subtiles_indexed_and_rendered_nearest(tmp_path): + """Couche densité : quadrants .webp sans perte reconnus comme palier fin, + et rendus au plus proche voisin (16 gris exacts même agrandis).""" + from PIL import Image + from lidar_pipeline import tiles + base = _make_dalle(tmp_path, 1054, 6882, ["densite_sol"]) + sub = tmp_path / "index_subtiles" + sub.mkdir(parents=True, exist_ok=True) + for i in range(2): + for j in range(2): + im = Image.new("L", (8, 8), 0) + im.paste(255, (0, 0, 4, 8)) # moitié blanche, moitié noire + im.save(str(sub / f"{base}_densite_sol_{i}_{j}.webp"), format="WEBP", lossless=True) + tiers = tiles.source_index(tmp_path, force=True)["densite_sol"][(1054, 6882)] + quads = next(t for t in tiers if len(t) == 4) + assert all(q.path.suffix == ".webp" for q in quads) + assert "densite_sol" in tiles.NEAREST_LAYERS + # z17 : résolution visée (~0,8 m) atteinte par les quadrants (0,5 m ici) + x, y = _tile_of_cell(1054, 6882, 17) + img = tiles.render_tile(tmp_path, "densite_sol", 17, x, y) + assert img is not None + grays = {p[0] for p in img.getdata() if p[3] == 255} + assert grays <= {0, 255} # aucun gris intermédiaire inventé diff --git a/lidar_pipeline/tests/test_visualizations.py b/lidar_pipeline/tests/test_visualizations.py index 0cc4be4..e598d2b 100644 --- a/lidar_pipeline/tests/test_visualizations.py +++ b/lidar_pipeline/tests/test_visualizations.py @@ -553,3 +553,30 @@ class TestPriorityFlood: result = _priority_flood(dem, nodata) assert result[2, 2] == 999.0 + + +class TestDensiteSol: + def test_density_levels_log_scale(self): + from lidar_pipeline.visualizations import density_levels + d = np.array([0.0, np.nan, 0.25, 0.36, 0.5, 1.0, 45.0, 45.3, 1000.0]) + assert density_levels(d).tolist() == [0, 0, 0, 1, 2, 4, 14, 15, 15] + + def test_generate_reads_density_sidecar(self, tmp_path): + import rasterio + from rasterio.transform import from_bounds + from lidar_pipeline.dtm import density_path + from lidar_pipeline.visualizations import generate_densite_sol + dem = tmp_path / "T_dtm_r0p2.tif" + dem.touch() + with rasterio.open(density_path(dem), 'w', driver='GTiff', width=4, height=1, + count=1, dtype='float32', crs='EPSG:2154', + transform=from_bounds(0, 0, 4, 1, 4, 1)) as dst: + dst.write(np.array([[0.0, 1.0, 4.0, 64.0]], dtype='float32'), 1) + out = generate_densite_sol(dem, "T", tmp_path, 0.2) + with rasterio.open(out) as src: + assert src.read(1).tolist() == [[0, 4, 8, 15]] + assert src.width == 4 # grille de 1 m conservée + + def test_missing_sidecar_returns_none(self, tmp_path): + from lidar_pipeline.visualizations import generate_densite_sol + assert generate_densite_sol(tmp_path / "X_dtm.tif", "X", tmp_path, 0.2) is None diff --git a/lidar_pipeline/tiles.py b/lidar_pipeline/tiles.py index 272daef..5c9c261 100644 --- a/lidar_pipeline/tiles.py +++ b/lidar_pipeline/tiles.py @@ -48,6 +48,9 @@ AVIF_SPEED = 9 # encodage rapide (cf. rendering.AVIF_SPEED : ×7, +3 % # l'active quand la bande passante compte (consultation mobile). PNG_PALETTE = os.environ.get("LIDAR_TILE_PNG_PALETTE", "") == "1" +# Couches en aplats de niveaux codés, rééchantillonnées au plus proche voisin +NEAREST_LAYERS = frozenset({'densite_sol'}) + # Demi-circonférence équatoriale : emprise du Web Mercator (EPSG:3857). ORIGIN = 20037508.342789244 @@ -452,8 +455,10 @@ def _build_index(output_dir): quads = {} for f in sub_dir.iterdir(): name = f.name + # « .webp » plein en dernier : les vignettes finissent aussi en .webp for suffix, tier in ((".avif", "full"), ("_mid.webp", "mid"), - (f"_thumb{_SUBTILE_THUMB_PX}.webp", "thumb")): + (f"_thumb{_SUBTILE_THUMB_PX}.webp", "thumb"), + (".webp", "full")): if not name.endswith(suffix): continue m = _SUBTILE_RE.search(name[:-len(suffix)] + ".avif") @@ -810,7 +815,9 @@ def render_tile(output_dir, layer, z, x, y, scale=1): 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 + # Aplats de niveaux codés (densité, 1 m/px) : plus proche voisin, sinon + # l'agrandissement aux zooms 18–19 invente des gris intermédiaires. + 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) diff --git a/lidar_pipeline/visualizations.py b/lidar_pipeline/visualizations.py index 5dfb823..d69161e 100644 --- a/lidar_pipeline/visualizations.py +++ b/lidar_pipeline/visualizations.py @@ -1129,6 +1129,48 @@ def generate_relief_oriente(dem_file, basename, vis_dir, resolution, shared=None return None +# Densité de points sol : 16 niveaux sur une échelle log fixe (même gris = +# même densité partout, mosaïque jointive). Le niveau k commence à +# DENSITY_BASE × 2^(k/2) pts/m² : 0 = ≤ 0,35 (noir, y compris sans aucun +# point), 15 = ≥ 45 (blanc) ; deux niveaux = densité doublée. +DENSITY_LEVELS = 16 +DENSITY_BASE = 0.25 + + +def density_levels(density): + """Niveau 0..15 d'une densité (pts/m²).""" + with np.errstate(divide="ignore", invalid="ignore"): + k = np.floor(2.0 * np.log2(np.asarray(density, dtype=np.float64) / DENSITY_BASE)) + return np.clip(np.nan_to_num(k, nan=0.0, neginf=0.0), 0, DENSITY_LEVELS - 1).astype(np.uint8) + + +def generate_densite_sol(dem_file, basename, vis_dir, resolution, shared=None): + """Densité des points sol retenus pour le MNT, en 16 niveaux de gris. + + Lue dans le fichier annexe écrit à la rastérisation (dtm.density_path) : + grille de 1 m, gardée telle quelle (la densité n'a pas plus de détail — + image 25× plus légère qu'à 0,2 m). Sortie : niveau 0..15 (float32). + """ + logger.info(" → Densité de points sol...") + t0 = time.time() + output = vis_dir / f"{basename}_densite_sol.tif" + try: + from .dtm import density_path + src_path = density_path(dem_file) + if not src_path.exists(): + logger.error(f" ✗ Densité absente ({src_path.name}) : régénérer le DTM") + return None + with rasterio.open(src_path) as src: + density = src.read(1) + transform, crs = src.transform, src.crs + _save_tif(output, density_levels(density).astype(np.float32), transform, crs) + logger.info(f" ✓ Densité de points terminée ({time.time()-t0:.1f}s)") + return output + except Exception as e: + logger.error(f" ✗ Erreur densité de points: {e}", exc_info=True) + return None + + def generate_mslrm(dem_file, basename, vis_dir, resolution, shared=None): """Multi-Scale Relief Model (MSRM) - LRM at adaptive scales combined (GPU if available).