np.percentile crash avec IndexError quand valid_data est vide (TIF entièrement NaN). Ajout d'un guard qui retourne une colormap par défaut si aucune donnée valide n'existe.
605 lines
27 KiB
Python
605 lines
27 KiB
Python
"""Rendering module: colormap registry, GeoTIFF-to-WebP conversion, and PDF report generation.
|
|
|
|
Contains:
|
|
- COLORMAPS: registry mapping filename keywords to (cmap, title, legend, description)
|
|
- tif_to_png(): convert a GeoTIFF to a WebP visualization with legend, scale bar, north arrow
|
|
- generate_pdf_report(): generate an A3 PDF report with all visualizations
|
|
"""
|
|
|
|
import logging
|
|
import time
|
|
from pathlib import Path
|
|
|
|
import numpy as np
|
|
import rasterio
|
|
from PIL import Image as PILImage
|
|
|
|
try:
|
|
from rasterio.warp import transform as warp_transform
|
|
HAS_WARP = True
|
|
except ImportError:
|
|
HAS_WARP = False
|
|
|
|
import matplotlib
|
|
matplotlib.use('Agg')
|
|
import matplotlib.pyplot as plt
|
|
from matplotlib import rcParams
|
|
from matplotlib.gridspec import GridSpec
|
|
from matplotlib.patches import Polygon as MplPolygon, Rectangle as RectPatch
|
|
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
|
|
|
|
rcParams['figure.dpi'] = 150
|
|
rcParams['savefig.dpi'] = 300
|
|
rcParams['font.size'] = 10
|
|
|
|
logger = logging.getLogger("lidar")
|
|
|
|
|
|
# ============================================================
|
|
# Colormap registry
|
|
# ============================================================
|
|
# Each entry: keyword → (cmap, title, legend_label, description, vmin_mode, vmax_mode)
|
|
# vmin_mode/vmax_mode: 'percentile_X_Y' or '0_max_X' or 'symmetric_X_Y' or 'fixed_0_1'
|
|
# For RGB images (ortho/topo), special handling is done in tif_to_png.
|
|
|
|
COLORMAPS = {
|
|
'hillshade': {
|
|
'cmap': 'gray',
|
|
'title': 'Hillshade Multidirectionnel',
|
|
'legend': 'Illumination combinée de 6 directions\nBlanc = Face éclairée | Noir = Zone d\'ombre',
|
|
'description': 'Ombres portées révélant micro-relief (murs, fossés, terrasses)',
|
|
'vmin_mode': 'fixed', 'vmin_val': 0,
|
|
'vmax_mode': 'fixed', 'vmax_val': 1,
|
|
},
|
|
'slope': {
|
|
'cmap': 'inferno',
|
|
'title': 'Pente (Inclinaison du terrain)',
|
|
'legend': 'Inclinaison en degrés\nMin: {vmin:.1f}° | Max: {vmax:.1f}°\nClair = Forte pente | Sombre = Terrain plat',
|
|
'description': 'Murs, talus et bords ressortent en clair — terrain plat en sombre',
|
|
'vmin_mode': 'fixed', 'vmin_val': 0,
|
|
'vmax_mode': 'percentile', 'vmax_pct': 95,
|
|
},
|
|
'aspect': {
|
|
'cmap': 'hsv',
|
|
'title': 'Aspect (Direction des pentes)',
|
|
'legend': 'Direction vers laquelle le terrain descend\nRouge=Nord | Vert=Est | Cyan=Sud | Bleu=Ouest',
|
|
'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),
|
|
},
|
|
'svf': {
|
|
'cmap': 'viridis',
|
|
'title': 'Sky-View Factor (Ray-tracing 16 directions)',
|
|
'legend': 'Fraction de ciel visible depuis chaque point\nSombre = Encaissé (fossés, vallées, rues)\nClair = Dégagé (sommets, plateformes, plateaux)',
|
|
'description': 'Ray-tracing sur 16 azimuts — élimine l\'éclairage, détecte structures linéaires et enclos',
|
|
'vmin_mode': 'percentile', 'vmin_pct': 5,
|
|
'vmax_mode': 'percentile', 'vmax_pct': 95,
|
|
},
|
|
'mslrm': {
|
|
'cmap': 'RdBu_r',
|
|
'title': 'MSRM - Multi-Scale Relief Model (5 échelles)',
|
|
'legend': 'Relief combiné à 5 échelles (5m à 100m)\nRouge = Surélévation (mur, tumulus, levée)\nBleu = Dépression (fossé, douve, fossé)\n\nDifférence avec LRM:\nLRM = 1 échelle (15m)\nMSRM = 5 échelles combiné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 (2 é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 échelle fine 5m + large 100m',
|
|
'description': 'Identifie la position topographique — utile pour repérer crêtes vs vallées à grande échelle',
|
|
'vmin_mode': 'symmetric', 'sym_pct': (2, 98),
|
|
},
|
|
'depressions': {
|
|
'cmap': 'YlOrRd',
|
|
'title': 'Dépressions (Remplissage hydrologique)',
|
|
'legend': 'Profondeur des cuvettes détectées (m)\nTransparent = Pas de dépression\nJaune = Dépression légère | Rouge = Dépression profonde\nMax: {vmax:.2f}m',
|
|
'description': 'Simule le remplissage d\'eau — détecte dolines, sinkholes, cuvettes et zones inondables',
|
|
'vmin_mode': 'fixed', 'vmin_val': 0,
|
|
'vmax_mode': 'percentile', 'vmax_pct': 99,
|
|
},
|
|
'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é de Surface (Écart-type local 5m)',
|
|
'legend': 'Irrégularité du terrain dans un voisinage de 5m\nSombre = Surface lisse (route, mur, sol plat)\nClair = Surface rugueuse (végétation, ruines, pierres)\nMax: {vmax:.2f}m',
|
|
'description': 'Mesure la variabilité locale — surfaces anthropiques lisses vs naturelles rugueuses',
|
|
'vmin_mode': 'fixed', 'vmin_val': 0,
|
|
'vmax_mode': 'percentile', 'vmax_pct': 97,
|
|
},
|
|
'anomalies': {
|
|
'cmap': 'coolwarm',
|
|
'title': 'Anomalies Statistiques (Z-score x Moran\'s I)',
|
|
'legend': 'Anomalies topographiques significatives\nRouge vif = Surélévation anormale (mur, tumulus)\nBleu vif = Dépression anormale (fossé, doline)\nBlanc/gris = Normal\n\nCombine z-score (intensité) et\nMoran\'s I (regroupement spatial)',
|
|
'description': 'Détecte uniquement les anomalies statistiquement significatives — filtre le bruit de fond',
|
|
'vmin_mode': 'symmetric', 'sym_pct': (5, 95),
|
|
},
|
|
'wavelet': {
|
|
'cmap': 'cividis',
|
|
'title': 'Ondelette Mexican Hat (CWT multi-échelle)',
|
|
'legend': 'Réponse de la transformée en ondelette à 5 échelles\nÉchelles: 2m, 5m, 10m, 20m, 50m\n\nClair = Structure détectée à cette échelle\nSombre = Pas de structure\n\nOptimisé pour formes circulaires:\ntumulus, enclos, fossés annulaires',
|
|
'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
|
|
RGB_LEGENDS = {
|
|
'ortho': {
|
|
'title': 'Photographie Aérienne IGN',
|
|
'legend': 'Orthophotographie\nImage aérienne',
|
|
'description': 'Photographie aérienne IGN (Orthophoto)',
|
|
},
|
|
'topo': {
|
|
'title': 'Carte Topographique IGN',
|
|
'legend': 'Carte IGN\nPlan topographique',
|
|
'description': 'Carte topographique IGN (Plan IGN)',
|
|
},
|
|
}
|
|
|
|
|
|
def _apply_colormap(data, tif_file):
|
|
"""Apply the registered colormap normalization to data based on filename.
|
|
|
|
Returns (data, cmap, title, legend_label, description, is_rgb).
|
|
"""
|
|
name = str(tif_file).lower()
|
|
|
|
# Check for RGB first
|
|
for key in RGB_LEGENDS:
|
|
if key in name:
|
|
info = RGB_LEGENDS[key]
|
|
return data, None, info['title'], info['legend'], info['description'], True
|
|
|
|
# Find matching colormap
|
|
for key, info in COLORMAPS.items():
|
|
if key in name:
|
|
valid_data = np.asarray(data.compressed() if hasattr(data, 'compressed') else data.flatten())
|
|
valid_data = valid_data[~np.isnan(valid_data)]
|
|
|
|
if len(valid_data) == 0:
|
|
logger.warning(f" Aucune donnée valide pour {Path(tif_file).name} — colormap ignorée")
|
|
return data, 'terrain', Path(tif_file).stem.replace('_', ' ').title(), '', '', False
|
|
|
|
vmin = vmax = None
|
|
|
|
# Compute vmin/vmax based on mode
|
|
if info['vmin_mode'] == 'fixed':
|
|
vmin = info['vmin_val']
|
|
elif info['vmin_mode'] == 'percentile':
|
|
vmin = np.percentile(valid_data, info['vmin_pct'])
|
|
elif info['vmin_mode'] == 'symmetric':
|
|
vmax_abs = max(abs(np.percentile(valid_data, info['sym_pct'][0])),
|
|
abs(np.percentile(valid_data, info['sym_pct'][1])), 0.001)
|
|
vmin = -vmax_abs
|
|
vmax = vmax_abs # symmetric mode sets both vmin and vmax
|
|
|
|
if vmax is None:
|
|
# Only compute vmax if not already set by symmetric mode
|
|
if info.get('vmax_mode') == 'fixed':
|
|
vmax = info['vmax_val']
|
|
elif info.get('vmax_mode') == 'percentile':
|
|
vmax = np.percentile(valid_data, info['vmax_pct'])
|
|
elif info.get('vmax_mode') == 'symmetric':
|
|
vmax_abs = max(abs(np.percentile(valid_data, info['sym_pct'][0])),
|
|
abs(np.percentile(valid_data, info['sym_pct'][1])), 0.001)
|
|
vmax = vmax_abs
|
|
|
|
# Apply normalization
|
|
if vmin is not None and vmax is not None:
|
|
data = np.clip((data - vmin) / max(vmax - vmin, 0.001), 0, 1)
|
|
|
|
legend = info['legend'].format(vmin=vmin or 0, vmax=vmax or 0)
|
|
return data, info['cmap'], info['title'], legend, info['description'], False
|
|
|
|
# Default: terrain colormap with percentile stretch
|
|
valid_data = np.asarray(data.compressed() if hasattr(data, 'compressed') else data.flatten())
|
|
valid_data = valid_data[~np.isnan(valid_data)]
|
|
if len(valid_data) == 0:
|
|
return data, 'terrain', Path(tif_file).stem.replace('_', ' ').title(), '', '', False
|
|
p2, p98 = np.percentile(valid_data, (2, 98))
|
|
data = np.clip((data - p2) / (p98 - p2), 0, 1)
|
|
title = Path(tif_file).stem.replace('_', ' ').title()
|
|
return data, 'terrain', title, 'Altitude normalisée', '', False
|
|
|
|
|
|
def tif_to_png(tif_file, vis_dir, resolution):
|
|
"""Convert GeoTIFF to visualization WebP with GPS coordinates, legend, and scale bar.
|
|
|
|
Args:
|
|
tif_file: Path to input GeoTIFF.
|
|
vis_dir: Output directory for the WebP file.
|
|
resolution: Grid resolution in m/px.
|
|
|
|
Returns:
|
|
Path to output WebP file, or None on failure.
|
|
"""
|
|
if not tif_file or not tif_file.exists():
|
|
return None
|
|
|
|
webp_file = vis_dir / f"{tif_file.stem}.webp"
|
|
|
|
try:
|
|
with rasterio.open(tif_file) as src:
|
|
is_rgb = src.count >= 3 and ('ortho' in str(tif_file) or 'topo' in str(tif_file))
|
|
|
|
if is_rgb:
|
|
data = src.read([1, 2, 3])
|
|
data = np.moveaxis(data, 0, -1)
|
|
else:
|
|
data = src.read(1)
|
|
|
|
nodata = src.nodata
|
|
transform = src.transform
|
|
crs = src.crs
|
|
|
|
if is_rgb:
|
|
height, width, _ = data.shape
|
|
else:
|
|
height, width = data.shape
|
|
|
|
top_left_x = transform.c
|
|
top_left_y = transform.f
|
|
pixel_size_x = transform.a
|
|
pixel_size_y = abs(transform.e)
|
|
|
|
min_x = top_left_x
|
|
max_x = top_left_x + width * pixel_size_x
|
|
max_y = top_left_y
|
|
min_y = top_left_y - height * pixel_size_y
|
|
|
|
# GPS coordinates
|
|
gps_coords = {}
|
|
if HAS_WARP and crs is not None:
|
|
try:
|
|
l93_xs = [min_x, max_x, min_x, max_x]
|
|
l93_ys = [max_y, max_y, min_y, min_y]
|
|
lons, lats = warp_transform(crs, 'EPSG:4326', l93_xs, l93_ys)
|
|
gps_coords = {
|
|
'NW': (lats[0], lons[0]),
|
|
'NE': (lats[1], lons[1]),
|
|
'SW': (lats[2], lons[2]),
|
|
'SE': (lats[3], lons[3]),
|
|
}
|
|
n_ticks = 5
|
|
tick_l93_x = np.linspace(min_x, max_x, n_ticks)
|
|
tick_l93_y_bottom = np.full(n_ticks, min_y)
|
|
tick_lons, tick_lats = warp_transform(crs, 'EPSG:4326', tick_l93_x, tick_l93_y_bottom)
|
|
gps_coords['x_ticks'] = list(zip(tick_lons, tick_lats))
|
|
tick_l93_y = np.linspace(min_y, max_y, n_ticks)
|
|
tick_l93_x_left = np.full(n_ticks, min_x)
|
|
tick_lons_y, tick_lats_y = warp_transform(crs, 'EPSG:4326', tick_l93_x_left, tick_l93_y)
|
|
gps_coords['y_ticks'] = list(zip(tick_lons_y, tick_lats_y))
|
|
except Exception:
|
|
gps_coords = {}
|
|
|
|
if nodata is not None and not is_rgb:
|
|
data = np.ma.masked_where((data == nodata) | np.isnan(data), data)
|
|
|
|
if not is_rgb:
|
|
valid_data = np.asarray(data.compressed() if hasattr(data, 'compressed') else data.flatten())
|
|
valid_data = valid_data[~np.isnan(valid_data)]
|
|
|
|
# Convert MaskedArray to plain ndarray (NaN-filled) to avoid numpy warnings
|
|
if isinstance(data, np.ma.MaskedArray):
|
|
data = np.ma.filled(data, np.nan)
|
|
|
|
# Apply colormap
|
|
data, cmap, title, legend_label, description, is_rgb_result = _apply_colormap(data, tif_file)
|
|
|
|
# Create figure — adapt width to data resolution for sharp rendering
|
|
# At high res (5000+px wide), we need a larger figure to avoid downsampling artifacts
|
|
fig_width = max(20, width / 150)
|
|
fig_width = min(fig_width, 40) # cap at 40 inches
|
|
map_aspect = height / width
|
|
fig = plt.figure(figsize=(fig_width, fig_width * map_aspect * 0.7 + 2.5),
|
|
facecolor='white')
|
|
gs = GridSpec(2, 1, height_ratios=[1.0, 0.06],
|
|
hspace=0.04, figure=fig,
|
|
left=0.06, right=0.88, top=0.93, bottom=0.08)
|
|
|
|
ax = fig.add_subplot(gs[0])
|
|
if is_rgb:
|
|
im = ax.imshow(data, aspect='equal', origin='upper',
|
|
interpolation='bilinear')
|
|
else:
|
|
im = ax.imshow(data, cmap=cmap, aspect='equal', origin='upper',
|
|
interpolation='bilinear')
|
|
|
|
ax.set_title(f"{title}\n{description}", fontsize=15, fontweight='bold', pad=10)
|
|
|
|
if not is_rgb:
|
|
cbar = plt.colorbar(im, ax=ax, pad=0.02, shrink=0.85, aspect=30)
|
|
cbar.ax.tick_params(labelsize=9, width=1.5)
|
|
cbar.outline.set_linewidth(1.5)
|
|
cbar.set_label(legend_label, fontsize=10, fontweight='bold')
|
|
else:
|
|
ax.text(1.02, 0.5, legend_label, transform=ax.transAxes,
|
|
fontsize=10, fontweight='bold', rotation=90,
|
|
verticalalignment='center', horizontalalignment='left')
|
|
|
|
# GPS coordinate ticks
|
|
if gps_coords and 'x_ticks' in gps_coords:
|
|
x_pixel_positions = np.linspace(0, width - 1, len(gps_coords['x_ticks']))
|
|
x_labels = [f"{lon:.5f}E" for lon, lat in gps_coords['x_ticks']]
|
|
ax.set_xticks(x_pixel_positions)
|
|
ax.set_xticklabels(x_labels, fontsize=7, rotation=30)
|
|
ax.set_xlabel('Longitude', fontsize=9, fontweight='bold')
|
|
|
|
y_pixel_positions = np.linspace(0, height - 1, len(gps_coords['y_ticks']))
|
|
y_labels = [f"{lat:.5f}N" for lon, lat in gps_coords['y_ticks']]
|
|
ax.set_yticks(y_pixel_positions)
|
|
ax.set_yticklabels(y_labels, fontsize=7)
|
|
ax.set_ylabel('Latitude', fontsize=9, fontweight='bold')
|
|
else:
|
|
x_ticks_count = 5
|
|
x_positions = np.linspace(0, width - 1, x_ticks_count)
|
|
x_labels = [f"{(min_x + xp * pixel_size_x)/1000:.1f}" for xp in x_positions]
|
|
ax.set_xticks(x_positions)
|
|
ax.set_xticklabels(x_labels, fontsize=8)
|
|
ax.set_xlabel('Est (km) - Lambert 93', fontsize=9, fontweight='bold')
|
|
|
|
y_ticks_count = 5
|
|
y_positions = np.linspace(0, height - 1, y_ticks_count)
|
|
y_labels = [f"{(max_y - yp * pixel_size_y)/1000:.1f}" for yp in y_positions]
|
|
ax.set_yticks(y_positions)
|
|
ax.set_yticklabels(y_labels, fontsize=8)
|
|
ax.set_ylabel('Nord (km) - Lambert 93', fontsize=9, fontweight='bold')
|
|
|
|
ax.tick_params(axis='both', which='both', direction='out', length=3,
|
|
width=0.8, colors='black')
|
|
for spine in ax.spines.values():
|
|
spine.set_visible(True)
|
|
spine.set_color('black')
|
|
spine.set_linewidth(0.8)
|
|
|
|
# North arrow
|
|
north_ax = inset_axes(ax, width="4%", height="7%", loc='upper right',
|
|
bbox_to_anchor=(-0.03, 0.12, 1, 1), bbox_transform=ax.transAxes)
|
|
north_ax.set_xlim(0, 1)
|
|
north_ax.set_ylim(0, 1)
|
|
north_ax.axis('off')
|
|
north_ax.plot([0.5, 0.5], [0.1, 0.65], color='black', linewidth=2.5, zorder=10)
|
|
north_ax.add_patch(MplPolygon([[0.5, 0.2], [0.35, 0.4], [0.5, 0.65], [0.65, 0.4]],
|
|
closed=True, facecolor='black', edgecolor='black', zorder=9))
|
|
north_ax.text(0.5, 0.95, 'N', ha='center', va='top',
|
|
fontsize=13, fontweight='bold', color='black', zorder=11)
|
|
|
|
# Bottom info bar
|
|
info_ax = fig.add_subplot(gs[1])
|
|
info_ax.axis('off')
|
|
|
|
extent_km_x = (max_x - min_x) / 1000
|
|
extent_km_y = (max_y - min_y) / 1000
|
|
|
|
if is_rgb:
|
|
alt_min = alt_max = 0
|
|
else:
|
|
alt_min = float(np.nanmin(valid_data)) if len(valid_data) > 0 else 0
|
|
alt_max = float(np.nanmax(valid_data)) if len(valid_data) > 0 else 0
|
|
|
|
if gps_coords:
|
|
nw_lat, nw_lon = gps_coords['NW']
|
|
se_lat, se_lon = gps_coords['SE']
|
|
info_text = (
|
|
f"GPS: {nw_lat:.5f}N {nw_lon:.5f}E - {se_lat:.5f}N {se_lon:.5f}E | "
|
|
f"EPSG:2154 | Res: {resolution}m/px | "
|
|
f"Emprise: {extent_km_x:.1f}x{extent_km_y:.1f}km"
|
|
)
|
|
else:
|
|
info_text = (
|
|
f"EPSG:2154 | X: {min_x:.0f}-{max_x:.0f} Y: {min_y:.0f}-{max_y:.0f} | "
|
|
f"Res: {resolution}m/px | Emprise: {extent_km_x:.1f}x{extent_km_y:.1f}km"
|
|
)
|
|
|
|
info_ax.text(0.01, 0.5, info_text,
|
|
transform=info_ax.transAxes, fontsize=8.5,
|
|
verticalalignment='center', family='monospace',
|
|
bbox=dict(boxstyle='round,pad=0.3', facecolor='#f0f0f0',
|
|
edgecolor='#aaaaaa', alpha=0.95))
|
|
|
|
# Scale bar
|
|
scale_m = 100
|
|
pixels_per_meter = 1.0 / pixel_size_x
|
|
scale_px = int(scale_m * pixels_per_meter)
|
|
scale_bar_frac = scale_px / width
|
|
bar_left = 0.80
|
|
bar_bottom = 0.15
|
|
bar_width_frac = min(scale_bar_frac, 0.15)
|
|
bar_height = 0.35
|
|
|
|
info_ax.add_patch(RectPatch((bar_left, bar_bottom), bar_width_frac, bar_height,
|
|
facecolor='black', edgecolor='black', linewidth=0.5,
|
|
transform=info_ax.transAxes, clip_on=False))
|
|
info_ax.text(bar_left + bar_width_frac / 2, bar_bottom + bar_height + 0.12,
|
|
f"{scale_m} m", ha='center', va='bottom', fontsize=9, fontweight='bold',
|
|
transform=info_ax.transAxes)
|
|
|
|
fig.patch.set_facecolor('white')
|
|
|
|
# Save as PNG then convert to WebP — use higher DPI for large data
|
|
save_dpi = 200 if width > 3000 else 150
|
|
png_temp = vis_dir / f"{tif_file.stem}_temp.png"
|
|
plt.savefig(png_temp, dpi=save_dpi, bbox_inches='tight', pad_inches=0.15,
|
|
facecolor='white', format='png')
|
|
plt.close()
|
|
|
|
img = PILImage.open(str(png_temp))
|
|
img.save(str(webp_file), format='WEBP', lossless=True)
|
|
png_temp.unlink()
|
|
|
|
# Delete source TIFF
|
|
tif_file.unlink()
|
|
|
|
return webp_file
|
|
|
|
except Exception as e:
|
|
logger.error(f" Erreur conversion WebP: {e}", exc_info=True)
|
|
return None
|
|
|
|
|
|
def generate_pdf_report(basename, vis_dir, pdf_dir, resolution):
|
|
"""Generate A3 PDF report for a LiDAR file with all visualizations.
|
|
|
|
Page 1: Mise en situation (ortho + topo IGN side by side)
|
|
Pages 2+: Other visualizations (2 per page)
|
|
|
|
Args:
|
|
basename: Base name for the report file.
|
|
vis_dir: Directory containing WebP visualization files.
|
|
pdf_dir: Directory for output PDF.
|
|
resolution: Grid resolution (used in info text).
|
|
|
|
Returns:
|
|
Path to PDF file, or None on failure.
|
|
"""
|
|
from matplotlib.backends.backend_pdf import PdfPages
|
|
|
|
pdf_file = pdf_dir / f"{basename}_rapport.pdf"
|
|
logger.info(f" → Génération rapport PDF A3: {pdf_file.name}")
|
|
t0 = time.time()
|
|
|
|
# Look for WebPs in per-file subdirectory first, then fallback to main dir
|
|
file_vis_dir = vis_dir / basename
|
|
if file_vis_dir.exists():
|
|
png_files = sorted(file_vis_dir.glob("*.webp"))
|
|
else:
|
|
png_files = sorted(vis_dir.glob(f"{basename}_*.webp"))
|
|
if not png_files:
|
|
logger.warning(f" ✗ Aucune image trouvée pour {basename}")
|
|
return None
|
|
|
|
# Categorize
|
|
situ_files = []
|
|
analysis_files = []
|
|
|
|
for f in png_files:
|
|
name = f.stem.lower()
|
|
if 'ortho' in name:
|
|
situ_files.insert(0, f)
|
|
elif 'topo' in name:
|
|
situ_files.append(f)
|
|
else:
|
|
analysis_files.append(f)
|
|
|
|
# Sort analysis files by archaeological priority
|
|
order = ['mslrm', 'svf', 'negative_openness',
|
|
'positive_openness', 'sailore', 'depressions', 'hillshade_multi',
|
|
'lrm', 'tpi', 'slope', 'curvature', 'aspect',
|
|
'roughness', 'anomalies', 'wavelet', 'flow']
|
|
|
|
def sort_key(f):
|
|
name = f.stem.lower()
|
|
for i, key in enumerate(order):
|
|
if key in name:
|
|
return i
|
|
return len(order)
|
|
|
|
analysis_files.sort(key=sort_key)
|
|
|
|
a3_w, a3_h = 16.54, 11.69
|
|
|
|
try:
|
|
with PdfPages(str(pdf_file)) as pdf:
|
|
# Page 1: Mise en situation
|
|
if situ_files:
|
|
fig = plt.figure(figsize=(a3_w, a3_h), facecolor='white')
|
|
n_situ = len(situ_files)
|
|
if n_situ == 2:
|
|
gs = fig.add_gridspec(1, 2, wspace=0.05, left=0.03, right=0.97,
|
|
top=0.92, bottom=0.06)
|
|
else:
|
|
gs = fig.add_gridspec(1, max(n_situ, 1), wspace=0.05,
|
|
left=0.03, right=0.97, top=0.92, bottom=0.06)
|
|
|
|
fig.text(0.5, 0.97, f"Mise en situation - {basename}",
|
|
fontsize=20, fontweight='bold', ha='center', va='top')
|
|
|
|
for i, f in enumerate(situ_files):
|
|
ax = fig.add_subplot(gs[0, i])
|
|
img = np.array(PILImage.open(str(f)).convert('RGB'))
|
|
ax.imshow(img)
|
|
ax.axis('off')
|
|
title = f.stem.replace(basename + '_', '').replace('_', ' ').title()
|
|
ax.set_title(title, fontsize=12, fontweight='bold', pad=5)
|
|
|
|
pdf.savefig(fig, dpi=150)
|
|
plt.close(fig)
|
|
|
|
# Pages 2+: Analysis maps (2 per page)
|
|
for page_start in range(0, len(analysis_files), 2):
|
|
page_files = analysis_files[page_start:page_start + 2]
|
|
|
|
fig = plt.figure(figsize=(a3_w, a3_h), facecolor='white')
|
|
|
|
if len(page_files) == 2:
|
|
gs = fig.add_gridspec(1, 2, wspace=0.08, left=0.03, right=0.97,
|
|
top=0.93, bottom=0.05)
|
|
else:
|
|
gs = fig.add_gridspec(1, 1, left=0.05, right=0.95,
|
|
top=0.93, bottom=0.05)
|
|
|
|
for i, f in enumerate(page_files):
|
|
ax = fig.add_subplot(gs[0, i])
|
|
img = np.array(PILImage.open(str(f)).convert('RGB'))
|
|
ax.imshow(img)
|
|
ax.axis('off')
|
|
title = f.stem.replace(basename + '_', '').replace('_', ' ').title()
|
|
ax.set_title(title, fontsize=11, fontweight='bold', pad=3)
|
|
|
|
page_num = (page_start // 2) + 2
|
|
fig.text(0.99, 0.01, f"Page {page_num}", fontsize=8,
|
|
ha='right', va='bottom', color='gray')
|
|
|
|
pdf.savefig(fig, dpi=150)
|
|
plt.close(fig)
|
|
|
|
logger.info(f" ✓ Rapport PDF terminé ({time.time()-t0:.1f}s)")
|
|
return pdf_file
|
|
|
|
except Exception as e:
|
|
logger.error(f" ✗ Erreur PDF: {e}", exc_info=True)
|
|
return None |