Upgrade colormaps to Crameri scientific palette and fix IGN location map
- Replace RdBu_r with Crameri roma for all 6 relief-family colormaps (curvature, mslrm, lrm, tpi, sailore, aniso_open): perceptually uniform, CVD-friendly, dark center makes near-zero values visible - Replace hsv with twilight for aspect: perceptually uniform cyclic colormap with no hue discontinuity - Add cmcrameri dependency (falls back to RdBu_r if unavailable) - Fix IGN location map: pass min_zoom=8 to download_ign_tiles so zoom 10 context maps work (was hardcoded to min_zoom=15, blocking the loop entirely for zoom_level=10) - Expand location map context from 3x extent to 80km fixed radius for better regional context at zoom 10 - Remove dead code: flow colormap entry, _nice_scale dead loop, aspect='auto' in location map imshow, outdated PDF report order list - Fix SVF docstring: 8 directions → 16 directions
This commit is contained in:
@ -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)
|
||||
|
||||
@ -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()
|
||||
|
||||
@ -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)
|
||||
|
||||
|
||||
Reference in New Issue
Block a user