diff --git a/Dockerfile b/Dockerfile index f46d1ba..3490119 100644 --- a/Dockerfile +++ b/Dockerfile @@ -42,7 +42,8 @@ RUN pip3 install --no-cache-dir \ fastapi \ uvicorn \ piexif \ - pillow-avif-plugin + pillow-avif-plugin \ + cmcrameri # Install CuPy for GPU acceleration (optional - will fallback to numpy if not available) RUN pip3 install --no-cache-dir cupy-cuda12x || echo "CuPy not available - GPU acceleration disabled" diff --git a/lidar_pipeline/ign.py b/lidar_pipeline/ign.py index 6d1d985..56138ab 100644 --- a/lidar_pipeline/ign.py +++ b/lidar_pipeline/ign.py @@ -73,7 +73,7 @@ def _lat_lon_to_px(lat, lon, zoom, tile_size=256): return px_x, px_y -def download_ign_tiles(min_x, max_x, min_y, max_y, layer, zoom_level=15): +def download_ign_tiles(min_x, max_x, min_y, max_y, layer, zoom_level=15, min_zoom=15): """Download IGN WMTS tiles for the given bounds using Web Mercator (PM). If the first tile returns 404, automatically retries at lower zoom levels. @@ -82,6 +82,7 @@ def download_ign_tiles(min_x, max_x, min_y, max_y, layer, zoom_level=15): min_x, max_x, min_y, max_y: Bounds in Lambert 93. layer: IGN WMTS layer name. zoom_level: WMTS zoom level (default 15). + min_zoom: Lowest zoom level to try when falling back (default 15). Returns: numpy array (H, W, 3) uint8, or None on failure. @@ -109,7 +110,6 @@ def download_ign_tiles(min_x, max_x, min_y, max_y, layer, zoom_level=15): tile_size = 256 # Try downloading at the requested zoom level; fall back to lower zooms on 404 - min_zoom = 15 for zoom in range(zoom_level, min_zoom - 1, -1): col_min, row_min = _lat_lon_to_tile(nw_lat, nw_lon, zoom) col_max, row_max = _lat_lon_to_tile(se_lat, se_lon, zoom) diff --git a/lidar_pipeline/rendering.py b/lidar_pipeline/rendering.py index 7c96aeb..54d6bdf 100644 --- a/lidar_pipeline/rendering.py +++ b/lidar_pipeline/rendering.py @@ -29,6 +29,13 @@ matplotlib.use('Agg') import matplotlib.pyplot as plt from matplotlib import rcParams from matplotlib.patches import Polygon as MplPolygon, Rectangle as RectPatch, FancyBboxPatch +from matplotlib.colors import ListedColormap + +try: + from cmcrameri import cm as cmc + HAS_CMCRAmeri = True +except ImportError: + HAS_CMCRAmeri = False rcParams['figure.dpi'] = 150 rcParams['savefig.dpi'] = 300 @@ -71,6 +78,77 @@ _FRANCE_OUTLINE_L93 = np.array([ # For RGB images (ortho/topo), special handling is done in tif_to_png. COLORMAPS = { + # === Famille RELIEF : rouge=surélévation, bleu=dépression === + # Roma (Crameri): perceptually uniform, CVD-friendly, dark center → near-zero values visible + # Falls back to RdBu_r if cmcrameri unavailable + 'curvature': { + 'cmap': 'roma' if HAS_CMCRAmeri else 'RdBu_r', + 'title': 'Courbure (Convexité/Concavité du terrain)', + 'legend': 'Changement de pente (1/m)\nRouge = Convexe (sommet de mur, levée)\nBleu = Concave (fond de fossé, dépression)', + 'description': 'Détecte les ruptures de pente — utile pour bords de terrasses et levées', + 'vmin_mode': 'symmetric', 'sym_pct': (5, 95), + }, + 'mslrm': { + 'cmap': 'roma' if HAS_CMCRAmeri else 'RdBu_r', + 'title': 'MSRM - Multi-Scale Relief Model (échelles adaptatives)', + 'legend': 'Relief combiné multi-échelles\nRouge = Surélévation (mur, tumulus, levée)\nBleu = Dépression (fossé, douve)\n\nLRM = 1 échelle (15m)\nMSRM = échelles combinées pondérées\nDétecte du micro au macro', + 'description': 'Combine LRM à 5 échelles — détecte structures de 5m à 100m simultanément', + 'vmin_mode': 'symmetric', 'sym_pct': (2, 98), + }, + 'lrm': { + 'cmap': 'roma' if HAS_CMCRAmeri else 'RdBu_r', + 'title': 'LRM - Local Relief Model (échelle unique 15m)', + 'legend': 'Écart local par rapport au terrain moyen (m)\nRouge = Surélévation (+{vmax:.2f}m)\nBleu = Dépression ({vmin:.2f}m)\nNoyau gaussien unique de 15m', + 'description': 'Micro-relief à 15m seulement — voir MSRM pour toutes les échelles', + 'vmin_mode': 'symmetric', 'sym_pct': (2, 98), + }, + 'tpi': { + 'cmap': 'roma' if HAS_CMCRAmeri else 'RdBu_r', + 'title': 'TPI - Topographic Position Index (4 échelles)', + 'legend': 'Position dans le paysage\nRouge = Plus haut que le voisinage (crête, plateau)\nBleu = Plus bas que le voisinage (fossé, vallée)\nCombine 4 échelles : 3m, 15m, 50m, 200m', + 'description': 'Identifie la position topographique — utile pour repérer crêtes vs vallées à grande échelle', + 'vmin_mode': 'symmetric', 'sym_pct': (2, 98), + }, + 'sailore': { + 'cmap': 'roma' if HAS_CMCRAmeri else 'RdBu_r', + 'title': 'SAILORE - LRM Auto-Adaptatif', + 'legend': 'Relief local adaptatif\nRouge = Surélévation | Bleu = Dépression\n\nLRM = noyau fixe 15m\nMSRM = 5 noyaux fixes\nSAILORE = noyau adapté à la pente\nPlat=grand noyau | Pente=petit noyau', + 'description': 'Noyau qui s\'adapte à la pente locale — terrain plat=grand noyau, pente=petit noyau', + 'vmin_mode': 'symmetric', 'sym_pct': (2, 98), + }, + 'aniso_open': { + 'cmap': 'roma' if HAS_CMCRAmeri else 'RdBu_r', + 'title': 'Openness Anisotropique (pondération directionnelle)', + 'legend': 'Openness positive - négative pondérée (degrés)\nRouge = Surélévation dominante (mur, levée)\nBleu = Dépression dominante (fossé, doline)\nPondère les directions NW-SE et NE-SW davantage', + 'description': 'Openness avec pondération anisotropique — détecte mieux les structures alignées NW-SE et NE-SW', + 'vmin_mode': 'symmetric', 'sym_pct': (2, 98), + }, + # === Famille OUVERTURE : clair=ouvert, sombre=fermé === + 'positive_openness': { + 'cmap': 'YlOrBr', + 'title': 'Openness Positive (ouverture vers le haut)', + 'legend': 'Angle d\'ouverture vers le haut (deg)\nClair = Vue dégagée vers le ciel (sommets, plateaux)\nSombre = Vue bloquée (vallées encaissées)', + 'description': 'Ray-tracing 8 directions — complémentaire de la négative pour détecter crêtes', + 'vmin_mode': 'percentile', 'vmin_pct': 10, + 'vmax_mode': 'percentile', 'vmax_pct': 98, + }, + 'negative_openness': { + 'cmap': 'PuBu', + 'title': 'Openness Negative (ouverture vers le bas)', + 'legend': 'Angle d\'ouverture vers le bas (deg)\nClair = Surplomb (bords de fossé, grottes)\nSombre = Terrain plat (fonds de vallée)\nMeilleur détecteur de cavités et dolines', + 'description': 'Ray-tracing 8 directions — détecte fossés, dolines, souterrains', + 'vmin_mode': 'percentile', 'vmin_pct': 10, + 'vmax_mode': 'percentile', 'vmax_pct': 98, + }, + 'svf': { + 'cmap': 'bone_r', + 'title': 'Sky-View Factor (fraction de ciel visible)', + 'legend': 'Proportion de ciel visible depuis chaque point\nClair = Ciel dégagé (sommet, plateau, levée)\nSombre = Ciel masqué (vallée, fossé, tranchée)\nContraste adapté aux valeurs réelles (percentiles 2-98)', + 'description': 'Détection de micro-relief — fossés sombres, levées claires, complémentaire de l\'openness', + 'vmin_mode': 'percentile', 'vmin_pct': 2, + 'vmax_mode': 'percentile', 'vmax_pct': 98, + }, + # === Famille SCALAIRE : propriétés non divergentes === 'hillshade': { 'cmap': 'gray', 'title': 'Hillshade Multidirectionnel', @@ -88,64 +166,13 @@ COLORMAPS = { 'vmax_mode': 'percentile', 'vmax_pct': 95, }, 'aspect': { - 'cmap': 'hsv', + 'cmap': 'twilight', 'title': 'Aspect (Direction des pentes)', - 'legend': 'Direction vers laquelle le terrain descend\nRouge=Nord | Vert=Est | Cyan=Sud | Bleu=Ouest', + 'legend': 'Direction vers laquelle le terrain descend\nCycle continu : Nord→Est→Sud→Ouest→Nord\nCouleurs perceptuellement uniformes (pas de saut de teinte)', 'description': 'Orientation des pentes — utile pour distinguer structures des formes naturelles', 'vmin_mode': 'fixed', 'vmin_val': 0, 'vmax_mode': 'fixed', 'vmax_val': 360, }, - 'curvature': { - 'cmap': 'RdYlBu_r', - 'title': 'Courbure (Convexité/Concavité du terrain)', - 'legend': 'Changement de pente (1/m)\nRouge = Convexe (sommet de mur, levée)\nBleu = Concave (fond de fossé, dépression)', - 'description': 'Détecte les ruptures de pente — utile pour bords de terrasses et levées', - 'vmin_mode': 'symmetric', 'sym_pct': (5, 95), - }, - 'mslrm': { - 'cmap': 'RdBu_r', - 'title': 'MSRM - Multi-Scale Relief Model (échelles adaptatives)', - 'legend': 'Relief combiné multi-échelles (adapté à la résolution)\nRouge = Surélévation (mur, tumulus, levée)\nBleu = Dépression (fossé, douve)\n\nDifférence avec LRM:\nLRM = 1 échelle (15m)\nMSRM = échelles combinées pondérées\nMSRM détecte du micro au macro', - 'description': 'Combine LRM à 5 échelles — détecte structures de 5m à 100m simultanément', - 'vmin_mode': 'symmetric', 'sym_pct': (2, 98), - }, - 'lrm': { - 'cmap': 'RdBu_r', - 'title': 'LRM - Local Relief Model (échelle unique 15m)', - 'legend': 'Écart local par rapport au terrain moyen (m)\nRouge = Surélévation (+{vmax:.2f}m)\nBleu = Dépression ({vmin:.2f}m)\nNoyau gaussien unique de 15m', - 'description': 'Micro-relief à 15m seulement — voir MSRM pour toutes les échelles', - 'vmin_mode': 'symmetric', 'sym_pct': (2, 98), - }, - 'positive_openness': { - 'cmap': 'YlOrBr', - 'title': 'Openness Positive (ouverture vers le haut)', - 'legend': 'Angle d\'ouverture vers le haut (deg)\nClair = Vue dégagée vers le ciel (sommets, plateaux)\nSombre = Vue bloquée (vallées encaissées)', - 'description': 'Ray-tracing 8 directions — complémentaire de la négative pour détecter crêtes', - 'vmin_mode': 'percentile', 'vmin_pct': 10, - 'vmax_mode': 'percentile', 'vmax_pct': 98, - }, - 'negative_openness': { - 'cmap': 'GnBu_r', - 'title': 'Openness Negative (ouverture vers le bas)', - 'legend': 'Angle d\'ouverture vers le bas (deg)\nClair = Surplomb (bords de fossé, grottes)\nSombre = Terrain plat (fonds de vallée)\nMeilleur détecteur de cavités et dolines', - 'description': 'Ray-tracing 8 directions — détecte fossés, dolines, souterrains', - 'vmin_mode': 'percentile', 'vmin_pct': 10, - 'vmax_mode': 'percentile', 'vmax_pct': 98, - }, - 'tpi': { - 'cmap': 'BrBG', - 'title': 'TPI - Topographic Position Index (4 échelles)', - 'legend': 'Position dans le paysage\nBrun/Sombre = Plus bas que le voisinage (fossé, vallée)\nVert/Clair = Plus haut que le voisinage (crête, plateau)\nCombine 4 échelles : 3m, 15m, 50m, 200m', - 'description': 'Identifie la position topographique — utile pour repérer crêtes vs vallées à grande échelle', - 'vmin_mode': 'symmetric', 'sym_pct': (2, 98), - }, - 'sailore': { - 'cmap': 'seismic', - 'title': 'SAILORE - LRM Auto-Adaptatif', - 'legend': 'Relief local adaptatif (m)\nRouge = Surélévation | Bleu = Dépression\n\nDifférence avec LRM/MSRM:\nLRM = noyau fixe 15m\nMSRM = 5 noyaux fixes\nSAILORE = noyau adapté à la pente\nPlat=grand noyau | Pente=petit noyau', - 'description': 'Noyau qui s\'adapte à la pente locale — terrain plat=grand noyau, pente=petit noyau', - 'vmin_mode': 'symmetric', 'sym_pct': (2, 98), - }, 'roughness': { 'cmap': 'magma', 'title': 'Rugosité Multi-Échelle (3m + 15m)', @@ -154,21 +181,6 @@ COLORMAPS = { 'vmin_mode': 'fixed', 'vmin_val': 0, 'vmax_mode': 'percentile', 'vmax_pct': 97, }, - 'svf': { - 'cmap': 'gray_r', - 'title': 'Sky-View Factor (fraction de ciel visible)', - 'legend': 'Proportion de ciel visible depuis chaque point\nBlanc = Ciel dégagé (sommet, plateau, levée)\nNoir = Ciel masqué (vallée, fossé, tranchée)\nMoyenne de cos²(angle horizon) sur 16 directions', - 'description': 'Détection de micro-relief — fossés sombres, levées claires, complémentaire de l\'openness', - 'vmin_mode': 'fixed', 'vmin_val': 0, - 'vmax_mode': 'fixed', 'vmax_val': 1, - }, - 'aniso_open': { - 'cmap': 'RdBu_r', - 'title': 'Openness Anisotropique (pondération directionnelle)', - 'legend': 'Openness positive - négative pondérée (degrés)\nRouge = Surélévation dominante (mur, levée)\nBleu = Dépression dominante (fossé, doline)\nPondère les directions NW-SE et NE-SW davantage', - 'description': 'Openness avec pondération anisotropique — détecte mieux les structures alignées NW-SE et NE-SW', - 'vmin_mode': 'symmetric', 'sym_pct': (2, 98), - }, 'wavelet': { 'cmap': 'cividis', 'title': 'Ondelette Mexican Hat (CWT multi-échelle)', @@ -176,14 +188,6 @@ COLORMAPS = { 'description': 'Transformée en ondelette 2D — excellente pour détecter structures circulaires', 'vmin_mode': 'symmetric', 'sym_pct': (2, 98), }, - 'flow': { - 'cmap': 'Blues', - 'title': 'Accumulation de Flux Hydrologique (D8)', - 'legend': 'Logarithme de l\'accumulation d\'eau\nBlanc = Pas de collecte (sommet, crête)\nBleu foncé = Collecte maximale (thalweg, fossé)\n\nSimule l\'écoulement de l\'eau de pluie\nDétecte fossés d\'enceinte, routes et drainage', - 'description': 'Algorithme D8 — simule le cheminement de l\'eau pour détecter fossés et routes antiques', - 'vmin_mode': 'fixed', 'vmin_val': 0, - 'vmax_mode': 'percentile', 'vmax_pct': 98, - }, } # RGB entries (ortho/topo) are handled specially @@ -300,22 +304,22 @@ def _download_location_map(min_x, max_x, min_y, max_y): center_lat = clats[0] center_lon = clons[0] - # Use a much lower zoom level for context (wider view) - # Zoom 10 gives ~150km per 256px tile — perfect for a small location map + # Use zoom 10 for context (wider view: ~150m/px) context_zoom = 10 - # Expand bounds by 3x in each direction for wider context - extent_x = max_x - min_x - extent_y = max_y - min_y - context_min_x = center_x - extent_x * 2 - context_max_x = center_x + extent_x * 2 - context_min_y = center_y - extent_y * 2 - context_max_y = center_y + extent_y * 2 + # Expand bounds to ~80km for regional context (shows nearby cities/rivers) + # At zoom 10, this gives ~500px wide image — ideal for a small inset + context_half = 40000 # 40km each side = 80km total + context_min_x = center_x - context_half + context_max_x = center_x + context_half + context_min_y = center_y - context_half + context_max_y = center_y + context_half result = download_ign_tiles( context_min_x, context_max_x, context_min_y, context_max_y, layer='GEOGRAPHICALGRIDSYSTEMS.PLANIGNV2', - zoom_level=context_zoom + zoom_level=context_zoom, + min_zoom=8 ) if result is not None: @@ -333,15 +337,7 @@ def _nice_scale(extent_m): Returns (scale_m, label) where label is like '100 m' or '500 m' or '1 km'. """ nice_scales = [50, 100, 200, 500, 1000, 2000, 5000, 10000] - # Pick the largest scale that is <= 15% of the extent - for s in nice_scales: - if s <= extent_m * 0.15: - continue - # s is the first one > 15% — take the one below - break - else: - s = nice_scales[-1] - # Actually pick the largest scale <= 20% of extent + # Pick the largest scale <= 20% of extent chosen = nice_scales[0] for s in nice_scales: if s <= extent_m * 0.20: @@ -676,7 +672,7 @@ def tif_to_png(tif_file, vis_dir, resolution, keep_tif=False, source_info=None, location_map = _download_location_map(min_x, max_x, min_y, max_y) if location_map is not None: # Draw IGN topo map as background - map_ax.imshow(location_map, aspect='auto', extent=[ + map_ax.imshow(location_map, aspect='equal', extent=[ min_x - (max_x - min_x) * 2, max_x + (max_x - min_x) * 2, min_y - (max_y - min_y) * 2, max_y + (max_y - min_y) * 2 ]) @@ -793,9 +789,9 @@ def generate_pdf_report(basename, vis_dir, pdf_dir, resolution): # Sort analysis files by archaeological priority order = ['mslrm', 'svf', 'negative_openness', - 'positive_openness', 'sailore', 'hillshade_multi', + 'positive_openness', 'aniso_open', 'sailore', 'hillshade_multi', 'lrm', 'tpi', 'slope', 'curvature', 'aspect', - 'roughness', 'anomalies', 'wavelet', 'flow'] + 'roughness', 'wavelet'] def sort_key(f): name = f.stem.lower() diff --git a/lidar_pipeline/visualizations.py b/lidar_pipeline/visualizations.py index 07991dc..f694cd6 100644 --- a/lidar_pipeline/visualizations.py +++ b/lidar_pipeline/visualizations.py @@ -1287,7 +1287,7 @@ def generate_flow(dem_file, basename, vis_dir, resolution, shared=None): def generate_svf(dem_file, basename, vis_dir, resolution, shared=None): """Sky-View Factor - fraction of sky visible from each point (GPU if available). - SVF = average of cos²(horizon_angle) across 8 directions. + SVF = average of cos²(horizon_angle) across 16 directions. High SVF (near 1) = open sky (ridgetop, plateau) Low SVF (near 0) = enclosed sky (valley, deep trench) diff --git a/requirements.txt b/requirements.txt index d9b122d..4dfe653 100644 --- a/requirements.txt +++ b/requirements.txt @@ -7,3 +7,4 @@ laspy[lazrs]>=2.5 scikit-image>=0.21 tqdm>=4.65 scikit-learn>=1.3 +cmcrameri>=2.0