Files
lidar_rendu/lidar_pipeline/rendering.py

977 lines
43 KiB
Python
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

"""Rendering module: colormap registry, GeoTIFF-to-image 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/AVIF visualization with legend, scale bar, north arrow
- generate_pdf_report(): generate an A3 PDF report with all visualizations
"""
import logging
import time
from datetime import datetime
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
# Cache for IGN location map tiles (avoid re-downloading for each visualization)
_location_map_cache = {}
import matplotlib
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
from matplotlib.ticker import ScalarFormatter
try:
from cmcrameri import cm as cmc
HAS_CMCRAmeri = True
except ImportError:
HAS_CMCRAmeri = False
rcParams['figure.dpi'] = 150
rcParams['savefig.dpi'] = 300
rcParams['font.size'] = 10
logger = logging.getLogger("lidar")
# ============================================================
# Simplified France outline in Lambert 93 (EPSG:2154)
# Used for location inset map on each visualization
# ============================================================
_FRANCE_OUTLINE_L93 = np.array([
[109000, 6385000], [134000, 6410000], [153000, 6430000], [173000, 6445000],
[200000, 6460000], [250000, 6475000], [300000, 6490000], [350000, 6500000],
[400000, 6505000], [450000, 6510000], [500000, 6510000], [550000, 6510000],
[600000, 6505000], [650000, 6500000], [700000, 6495000], [750000, 6485000],
[800000, 6470000], [840000, 6460000], [880000, 6450000], [920000, 6435000],
[950000, 6425000], [980000, 6415000], [1010000, 6405000], [1040000, 6395000],
[1060000, 6385000], [1080000, 6370000], [1100000, 6355000], [1120000, 6340000],
[1140000, 6320000], [1160000, 6300000], [1175000, 6280000], [1185000, 6260000],
[1190000, 6240000], [1195000, 6220000], [1198000, 6200000], [1196000, 6180000],
[1192000, 6160000], [1185000, 6140000], [1175000, 6120000], [1160000, 6100000],
[1140000, 6085000], [1120000, 6070000], [1095000, 6060000], [1070000, 6050000],
[1040000, 6040000], [1000000, 6035000], [950000, 6035000], [900000, 6035000],
[850000, 6040000], [800000, 6045000], [750000, 6050000], [700000, 6055000],
[650000, 6060000], [600000, 6065000], [550000, 6070000], [500000, 6075000],
[450000, 6080000], [400000, 6085000], [350000, 6095000], [300000, 6110000],
[250000, 6125000], [200000, 6145000], [160000, 6170000], [130000, 6200000],
[110000, 6230000], [100000, 6260000], [95000, 6290000], [100000, 6310000],
[105000, 6340000], [109000, 6385000],
])
# ============================================================
# Colormap registry
# ============================================================
# Each entry: keyword → (cmap, 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.
# Les textes (title/legend/description) viennent de VIZ_LEGENDS (index.py,
# source unique partagée avec l'export multi-dalles) — fusion ci-dessous.
COLORMAPS = {
# === Famille RELIEF : rouge=surélévation, bleu=dépression ===
# Diverging: rouge vif=positif, bleu vif=négatif, blanc=plat
'mslrm': {
'cmap': 'seismic',
'vmin_mode': 'fixed', 'vmin_val': -3,
'vmax_mode': 'fixed', 'vmax_val': 3,
},
'sailore': {
'cmap': 'seismic',
'vmin_mode': 'fixed', 'vmin_val': -3,
'vmax_mode': 'fixed', 'vmax_val': 3,
},
# === Famille OUVERTURE : séquentiel, valeurs normalisées ===
'positive_openness': {
'cmap': 'YlOrBr',
'vmin_mode': 'fixed', 'vmin_val': -3,
'vmax_mode': 'fixed', 'vmax_val': 3,
},
'negative_openness': {
'cmap': 'PuBu',
'vmin_mode': 'fixed', 'vmin_val': -3,
'vmax_mode': 'fixed', 'vmax_val': 3,
},
'svf': {
'cmap': 'hot_r',
'vmin_mode': 'fixed', 'vmin_val': 0,
'vmax_mode': 'fixed', 'vmax_val': 1,
},
# === Famille SCALAIRE : propriétés non divergentes ===
# Clé = mot-clé exact du nom de fichier de sortie (hillshade_multi).
'hillshade_multi': {
'cmap': 'gray',
'vmin_mode': 'fixed', 'vmin_val': 0,
'vmax_mode': 'fixed', 'vmax_val': 1,
},
'slope': {
'cmap': 'inferno',
'vmin_mode': 'fixed', 'vmin_val': 0,
'vmax_mode': 'fixed', 'vmax_val': 30,
},
'aspect': {
'cmap': 'twilight',
'vmin_mode': 'fixed', 'vmin_val': 0,
'vmax_mode': 'fixed', 'vmax_val': 360,
},
'roughness': {
'cmap': 'plasma',
# vmax fixe (p98 médian mesuré sur 20 dalles réelles) : un vmax au
# percentile par tuile rendait l'échelle non jointive entre tuiles.
'vmin_mode': 'fixed', 'vmin_val': 0,
'vmax_mode': 'fixed', 'vmax_val': 3.8,
},
'wavelet': {
'cmap': 'inferno',
# Nœuds (percentile → valeur) mesurés sur 20 tuiles réelles par
# résolution avec l'algorithme actuel (détendage 35 m, échelles
# 1-50 m) : médiane inter-tuiles des percentiles par tuile, robuste
# aux tuiles très structurées. Distribution plus large qu'avant
# détendage : le fond macro-relief supprimé, les petites structures
# percent bien plus au-dessus du bruit (p98 ≈ 9 vs 1,65 avant).
'knots': {
0.5: ([0.089, 0.213, 0.296, 0.445, 0.602, 0.783, 1.0, 1.274,
1.661, 2.32, 3.842, 5.661, 8.936, 11.93, 15.06],
[0.01, 0.05, 0.10, 0.20, 0.30, 0.40, 0.50, 0.60,
0.70, 0.80, 0.90, 0.95, 0.98, 0.99, 0.995]),
0.2: ([0.088, 0.213, 0.296, 0.445, 0.602, 0.783, 1.0, 1.287,
1.694, 2.358, 3.846, 5.687, 8.971, 11.928, 14.971],
[0.01, 0.05, 0.10, 0.20, 0.30, 0.40, 0.50, 0.60,
0.70, 0.80, 0.90, 0.95, 0.98, 0.99, 0.995]),
},
},
'flow_acc': {
'cmap': 'YlGn',
'vmin_mode': 'fixed', 'vmin_val': 0,
'vmax_mode': 'percentile', 'vmax_pct': 98,
},
'anomaly': {
'cmap': 'YlOrRd',
'vmin_mode': 'fixed', 'vmin_val': 0,
'vmax_mode': 'fixed', 'vmax_val': 1,
},
'solar': {
'cmap': 'gray',
'vmin_mode': 'fixed', 'vmin_val': 0,
'vmax_mode': 'fixed', 'vmax_val': 1,
},
}
# RGB entries (ortho/topo) are handled specially
RGB_LEGENDS = {
'ortho': {},
'topo': {},
}
# Fusion des textes de légende (titre / lecture du rendu / méthode de calcul)
# depuis la source unique VIZ_LEGENDS (index.py, sans dépendance lourde).
from .index import VIZ_LEGENDS
for _key, _info in COLORMAPS.items():
_info.update(VIZ_LEGENDS[_key])
for _key, _info in RGB_LEGENDS.items():
_info.update(VIZ_LEGENDS[_key])
def _apply_colormap(data, tif_file, resolution=None):
"""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 — sort by key length descending so 'mslrm' matches before 'lrm'
for key in sorted(COLORMAPS.keys(), key=len, reverse=True):
info = COLORMAPS[key]
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
knots = info.get('knots')
if knots is not None:
# Étalonnage quantile figé (appariement d'histogramme global,
# cf. normalisation radiométrique des mosaïques) : fonction de
# transfert en nœuds mesurés une fois sur un échantillon de
# tuiles — même valeur → même couleur sur toutes les tuiles,
# toute la palette utilisée, insensible aux queues locales.
if isinstance(knots, dict):
key = min(knots, key=lambda r: abs(float(r) - float(resolution or 0.5)))
knots = knots[key]
kv, kt = knots
data = np.interp(np.asarray(data, dtype=float), kv, kt, left=0.0, right=1.0)
vmin, vmax = kv[0], kv[-1]
else:
# 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))
# Garde contre les tuiles quasi constantes : plage nulle → division par zéro
span = max(p98 - p2, 1e-6)
data = np.clip((data - p2) / span, 0, 1)
title = Path(tif_file).stem.replace('_', ' ').title()
return data, 'terrain', title, 'Altitude normalisée', '', False
def _download_location_map(min_x, max_x, min_y, max_y):
"""Download a wide-area IGN topographic map for location context.
Downloads a zoomed-out IGN PLANIGNV2 tile covering ~30km around the
processed zone, giving regional context. Results are cached to
avoid re-downloading for each visualization in the same tile.
Args:
min_x, max_x, min_y, max_y: DTM bounds in Lambert 93.
Returns:
Tuple (image_array, bounds_dict) where bounds_dict has keys
'min_x', 'max_x', 'min_y', 'max_y' in Lambert 93, or None on failure.
"""
import hashlib
# Cache key based on rounded coordinates (1km grid)
cache_key = (round(min_x, -3), round(max_x, -3), round(min_y, -3), round(max_y, -3))
if cache_key in _location_map_cache:
return _location_map_cache[cache_key]
from .ign import download_ign_tiles, _optimal_zoom_level
if not HAS_WARP:
return None
try:
# Compute center coordinates for zoom calculation
center_x = (min_x + max_x) / 2
center_y = (min_y + max_y) / 2
clons, clats = warp_transform('EPSG:2154', 'EPSG:4326', [center_x], [center_y])
center_lat = clats[0]
center_lon = clons[0]
# Use zoom 10 for context (~150m/px — wide view, fast download)
context_zoom = 10
# Expand bounds to ~30km for regional context
# Large enough to see surrounding towns/rivers, small enough that
# the processed zone (typically 1km) is clearly visible as a red rectangle
context_half = 15000 # 15km each side = 30km 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,
min_zoom=8
)
if result is not None:
bounds = {
'min_x': context_min_x, 'max_x': context_max_x,
'min_y': context_min_y, 'max_y': context_max_y,
}
cached = (result, bounds)
_location_map_cache[cache_key] = cached
return cached
return result
except Exception as e:
logger.debug(f" Carte de localisation IGN non disponible: {e}")
return None
def _nice_scale(extent_m):
"""Choose a nice round scale distance that fits well in the image.
Returns (scale_m, label) where label is like '100 m' or '500 m' or '1 km'.
"""
nice_scales = [10, 20, 50, 100, 200, 500, 1000, 2000, 5000, 10000]
# Pick the largest scale that fits within 20% of extent
max_scale = extent_m * 0.20
chosen = nice_scales[0]
for s in nice_scales:
if s <= max_scale:
chosen = s
else:
break
if chosen >= 1000:
return chosen, f"{chosen // 1000} km"
return chosen, f"{chosen} m"
def tif_to_png(tif_file, vis_dir, resolution, keep_tif=False, source_info=None, quality=98, output_format='avif'):
"""Convert GeoTIFF to visualization image (WebP or AVIF) with GPS coordinates, legend, and scale bar.
Args:
tif_file: Path to input GeoTIFF.
vis_dir: Output directory for the image file.
resolution: Grid resolution in m/px.
keep_tif: If True, keep the source TIFF after conversion.
source_info: Dict with method/date/basename for metadata.
quality: Image quality (1-100). Use 100 for lossless. Default 95.
output_format: Output format ('webp' or 'avif'). Default 'webp'.
Returns:
Path to output image file, or None on failure.
"""
if not tif_file or not tif_file.exists():
return None
ext = 'avif' if output_format == 'avif' else 'webp'
output_file = vis_dir / f"{tif_file.stem}.{ext}"
try:
with rasterio.open(tif_file) as src:
is_rgb = src.count >= 3 and any(k in str(tif_file) for k in ('ortho', 'topo'))
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)]
# Track NaN mask before converting to plain ndarray
nan_mask = None
if not is_rgb:
if isinstance(data, np.ma.MaskedArray):
nan_mask = data.mask.copy()
data = np.ma.filled(data, np.nan)
elif np.any(np.isnan(data)):
nan_mask = np.isnan(data)
# For rendering: replace NaN with neutral value to avoid interpolation halos
if nan_mask is not None and np.any(nan_mask) and len(valid_data) > 0:
fill_value = float(np.median(valid_data))
data[nan_mask] = fill_value
nan_mask = nan_mask # keep for later
# Apply colormap
data, cmap, title, legend_label, description, is_rgb_result = _apply_colormap(data, tif_file, resolution=resolution)
# Apply NaN mask: make zones without data transparent
has_nan_mask = nan_mask is not None and not is_rgb_result
if has_nan_mask:
# data is normalized 0-1 from _apply_colormap; apply cmap to get RGBA
# Save the colormap for colorbar before converting to RGBA
saved_cmap = plt.get_cmap(cmap) if isinstance(cmap, str) else cmap
saved_vmin = float(np.nanmin(data)) if not nan_mask.all() else 0
saved_vmax = float(np.nanmax(data)) if not nan_mask.all() else 1
rgba = saved_cmap(data) # (H, W, 4) float RGBA
rgba[nan_mask, 3] = 0.0 # transparent where no data
data = rgba
is_rgba = True
else:
is_rgba = False
saved_cmap = None
saved_vmin = None
saved_vmax = None
# Create figure with FIXED layout for consistent data area position
# All visualizations use the same axes positions so they can be overlaid
fig_width = max(20, width / 150)
fig_width = min(fig_width, 40)
fig_height = fig_width * 0.7 + 2.0 # Fixed header + footer space
fig = plt.figure(figsize=(fig_width, fig_height), facecolor='white')
# Fixed data area position — identical for ALL visualization types
# This ensures overlay/superposition works across all WebP images
data_left = 0.08
data_bottom = 0.19
data_width_frac = 0.74
data_height_frac = 0.71
ax = fig.add_axes([data_left, data_bottom, data_width_frac, data_height_frac])
if is_rgba or 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}", fontsize=14, fontweight='bold', pad=10)
if description:
ax.text(0.5, 1.04, description, transform=ax.transAxes,
fontsize=10, fontstyle='italic', color='#555555',
ha='center', va='bottom')
# Colorbar/legend area — full height alongside data
cbar_left = data_left + data_width_frac + 0.02
cbar_width = 0.04
cbar_bottom = data_bottom
cbar_height = data_height_frac
if is_rgb:
# RGB: descriptive text label instead of gradient colorbar
cbar_ax = fig.add_axes([cbar_left, cbar_bottom, cbar_width, cbar_height])
cbar_ax.set_xticks([])
cbar_ax.set_yticks([])
cbar_ax.text(0.5, 0.5, legend_label, transform=cbar_ax.transAxes,
fontsize=9, fontweight='bold', rotation=90,
verticalalignment='center', horizontalalignment='center',
wrap=True)
cbar_ax.set_frame_on(False)
elif is_rgba and saved_cmap is not None:
cbar_ax = fig.add_axes([cbar_left, cbar_bottom, cbar_width, cbar_height])
sm = plt.cm.ScalarMappable(cmap=saved_cmap,
norm=plt.Normalize(vmin=saved_vmin, vmax=saved_vmax))
sm.set_array([])
cbar = plt.colorbar(sm, cax=cbar_ax)
cbar.ax.tick_params(labelsize=9, width=1.5)
cbar.ax.yaxis.set_major_formatter(ScalarFormatter(useOffset=False))
cbar.outline.set_linewidth(1.5)
cbar.set_label(legend_label, fontsize=10, fontweight='bold')
else:
cbar_ax = fig.add_axes([cbar_left, cbar_bottom, cbar_width, cbar_height])
cbar = plt.colorbar(im, cax=cbar_ax)
cbar.ax.tick_params(labelsize=9, width=1.5)
cbar.ax.yaxis.set_major_formatter(ScalarFormatter(useOffset=False))
cbar.outline.set_linewidth(1.5)
cbar.set_label(legend_label, fontsize=10, fontweight='bold')
# 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 — compass rose in bottom-right corner of data area
# Semi-transparent background for readability over any data
north_ax = fig.add_axes([data_left + data_width_frac - 0.07,
data_bottom + 0.01,
0.06, 0.14],
facecolor='none')
north_ax.set_xlim(-1.5, 1.5)
north_ax.set_ylim(-1.5, 1.5)
north_ax.axis('off')
north_ax.set_aspect('equal')
# Compass rose centered at (0, 0) — all 4 cardinals equidistant from center
# Semi-transparent white background circle
circle_bg = plt.Circle((0, 0), 1.0, facecolor='white', edgecolor='#888888',
linewidth=0.5, alpha=0.7, zorder=1)
north_ax.add_patch(circle_bg)
# N arrow (pointing up = North)
north_ax.annotate('N', xy=(0, 1.35), fontsize=9, fontweight='bold',
ha='center', va='bottom', color='#b22222', zorder=10)
north_ax.plot([0, 0], [-0.5, 1.0], color='#b22222', linewidth=2.0, zorder=10)
north_ax.add_patch(MplPolygon([[0, 0.5], [-0.2, 0.7], [0, 1.0], [0.2, 0.7]],
closed=True, facecolor='#b22222', edgecolor='#b22222', zorder=9))
# Cardinal ticks — all centered at (0, 0)
for angle, label in [(90, 'N'), (0, 'E'), (180, 'O'), (270, 'S')]:
rad = np.radians(angle)
north_ax.plot([1.0*np.cos(rad), 1.2*np.cos(rad)],
[1.0*np.sin(rad), 1.2*np.sin(rad)],
color='#555555', linewidth=0.8, zorder=5)
if label:
north_ax.text(1.35*np.cos(rad), 1.35*np.sin(rad), label,
fontsize=6, ha='center', va='center', color='#555555', zorder=5)
# Bottom info bar — enriched with source, method, date
info_ax = fig.add_axes([data_left, 0.015, data_width_frac + cbar_width + 0.02, 0.09])
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
# Build info lines
line1_parts = []
if gps_coords:
nw_lat, nw_lon = gps_coords['NW']
se_lat, se_lon = gps_coords['SE']
line1_parts.append(f"GPS: {nw_lat:.5f}°N {nw_lon:.5f}°E — {se_lat:.5f}°N {se_lon:.5f}°E")
else:
line1_parts.append(f"X: {min_x:.0f}–{max_x:.0f} Y: {min_y:.0f}–{max_y:.0f}")
line1_parts.append(f"EPSG:2154")
# Round resolution to avoid ugly decimals like 0.499999959
res_display = round(resolution, 2) if resolution < 1 else round(resolution, 1)
line1_parts.append(f"Res: {res_display}m/px")
line1_parts.append(f"Emprise: {extent_km_x:.1f}×{extent_km_y:.1f}km")
if not is_rgb:
line1_parts.append(f"Alt: {alt_min:.1f}–{alt_max:.1f}m")
line2_parts = []
line2_parts.append("Source: LiDAR HD IGN")
if source_info:
if source_info.get('method'):
line2_parts.append(f"Classif.: {source_info['method'].upper()}")
if source_info.get('date'):
line2_parts.append(f"Date: {source_info['date']}")
else:
line2_parts.append(datetime.now().strftime("Date: %Y-%m-%d"))
info_text_line1 = " | ".join(line1_parts)
info_text_line2 = " | ".join(line2_parts)
info_ax.text(0.01, 0.7, info_text_line1,
transform=info_ax.transAxes, fontsize=8,
verticalalignment='center', family='monospace',
bbox=dict(boxstyle='round,pad=0.2', facecolor='#f0f0f0',
edgecolor='#aaaaaa', alpha=0.95))
info_ax.text(0.01, 0.2, info_text_line2,
transform=info_ax.transAxes, fontsize=7.5,
verticalalignment='center', family='monospace',
color='#444444',
bbox=dict(boxstyle='round,pad=0.2', facecolor='#f8f8f8',
edgecolor='#cccccc', alpha=0.9))
# Scale bar — adaptive with alternating black/white segments
# Position: left of the location map to avoid overlap
extent_m_x = max_x - min_x
scale_m, scale_label = _nice_scale(extent_m_x)
pixels_per_meter = 1.0 / pixel_size_x
scale_px = int(scale_m * pixels_per_meter)
n_segments = 5
segment_px = scale_px / n_segments
bar_bottom_y = 0.55
bar_top_y = 0.85
bar_height = bar_top_y - bar_bottom_y
# Place scale bar so it ends before the location map (map starts at x=0.82 in fig coords)
# map_ax occupies [0.82, 0.02, 0.16, 0.13] in figure coords
# info_ax occupies [data_left, 0.015, width, 0.09]
# Scale bar end in fig coords = info_ax.left + (scale_start_x + scale_px/width) * info_ax.width
# We need: info_ax.left + (scale_start_x + scale_px/width) * info_ax.width < 0.80
scale_end_frac = scale_px / width # fraction of info_ax width
info_ax_width = data_width_frac + cbar_width + 0.02
# Calculate scale_start_x so scale bar ends at fig_x = 0.78 (leaving gap before map at 0.82)
max_scale_end_fig = 0.78
scale_end_in_info = (max_scale_end_fig - data_left) / info_ax_width
scale_start_x = max(0.05, scale_end_in_info - scale_end_frac)
for seg_i in range(n_segments):
color = 'black' if seg_i % 2 == 0 else 'white'
seg_left = scale_start_x + seg_i * segment_px / width
seg_width_frac = segment_px / width
info_ax.add_patch(RectPatch((seg_left, bar_bottom_y), seg_width_frac, bar_height,
facecolor=color, edgecolor='black', linewidth=0.5,
transform=info_ax.transAxes, clip_on=False))
info_ax.text(scale_start_x + scale_px / (2 * width), bar_top_y + 0.12,
f"{scale_label}", ha='center', va='bottom', fontsize=8, fontweight='bold',
transform=info_ax.transAxes)
# Scale end ticks
info_ax.plot([scale_start_x, scale_start_x], [bar_bottom_y - 0.05, bar_top_y + 0.05],
color='black', linewidth=1, transform=info_ax.transAxes, clip_on=False)
info_ax.plot([scale_start_x + scale_px / width, scale_start_x + scale_px / width],
[bar_bottom_y - 0.05, bar_top_y + 0.05],
color='black', linewidth=1, transform=info_ax.transAxes, clip_on=False)
# Location inset map — IGN topographic background with processed zone marker
# Positioned in lower-right corner, above the info bar
map_ax = fig.add_axes([0.82, 0.02, 0.16, 0.13])
# Try to download a wide-area IGN topo map for location context
location_result = _download_location_map(min_x, max_x, min_y, max_y)
if location_result is not None:
location_map, loc_bounds = location_result
# Draw IGN topo map as background with correct bounds
map_ax.imshow(location_map, aspect='equal', extent=[
loc_bounds['min_x'], loc_bounds['max_x'],
loc_bounds['min_y'], loc_bounds['max_y']
])
# Mark the processed zone with a red rectangle
rect_x1, rect_x2 = min_x, max_x
rect_y1, rect_y2 = min_y, max_y
map_ax.add_patch(RectPatch((rect_x1, rect_y1),
rect_x2 - rect_x1, rect_y2 - rect_y1,
facecolor='#ff3333', edgecolor='#cc0000',
linewidth=1.5, alpha=0.6, zorder=5))
else:
# Fallback: simplified France outline
map_ax.set_facecolor('#e8e8e8')
france = _FRANCE_OUTLINE_L93
map_ax.fill(france[:, 0] / 1000, france[:, 1] / 1000,
facecolor='#f5f0e6', edgecolor='#888888', linewidth=0.8)
rect_x1, rect_x2 = min_x / 1000, max_x / 1000
rect_y1, rect_y2 = min_y / 1000, max_y / 1000
map_ax.add_patch(RectPatch((rect_x1, rect_y1),
rect_x2 - rect_x1, rect_y2 - rect_y1,
facecolor='#ff3333', edgecolor='#cc0000',
linewidth=1.2, alpha=0.7, zorder=5))
map_ax.set_xlim(france[:, 0].min() / 1000 - 50, france[:, 0].max() / 1000 + 50)
map_ax.set_ylim(france[:, 1].min() / 1000 - 50, france[:, 1].max() / 1000 + 50)
map_ax.set_aspect('equal')
map_ax.tick_params(left=False, bottom=False, labelleft=False, labelbottom=False)
for spine in map_ax.spines.values():
spine.set_edgecolor('#aaaaaa')
spine.set_linewidth(0.5)
# Label with coordinates
if gps_coords:
nw_lat, nw_lon = gps_coords['NW']
se_lat, se_lon = gps_coords['SE']
map_ax.set_title(f"{nw_lat:.2f}°N {nw_lon:.2f}°E",
fontsize=6, pad=1, color='#333333')
else:
map_ax.set_title(f"X:{min_x/1000:.0f} Y:{min_y/1000:.0f} km L93",
fontsize=6, pad=1, color='#333333')
fig.patch.set_facecolor('white')
# Save figure to in-memory buffer (avoids disk I/O of temp PNG)
save_dpi = 200 if width > 3000 else 150
from io import BytesIO
buf = BytesIO()
try:
plt.savefig(buf, dpi=save_dpi, facecolor='white', format='png')
finally:
plt.close()
buf.seek(0)
img = PILImage.open(buf)
pil_format = 'AVIF' if output_format == 'avif' else 'WEBP'
if quality >= 100:
img.save(str(output_file), format=pil_format, lossless=True)
else:
img.save(str(output_file), format=pil_format, quality=quality)
# Delete source TIFF (unless --keep-tif)
if not keep_tif:
tif_file.unlink(missing_ok=True)
return output_file
except Exception as e:
logger.error(f" Erreur conversion {ext.upper()}: {e}", exc_info=True)
return None
def tif_to_crop(tif_file, vis_dir, resolution, keep_tif=False, quality=98, output_format='avif'):
"""Convert GeoTIFF to a cropped visualization image (no legend, no overlay).
Applies colormap and saves the image as a pure 1×1 km square.
Used for grid/map display where images must tile seamlessly.
Args:
tif_file: Path to input GeoTIFF.
vis_dir: Output directory for the image file.
resolution: Grid resolution in m/px.
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').
Returns:
Path to output image file, or None on failure.
"""
if not tif_file or not tif_file.exists():
return None
ext = 'avif' if output_format == 'avif' else 'webp'
output_file = vis_dir / f"{tif_file.stem}.{ext}"
try:
with rasterio.open(tif_file) as src:
is_rgb = src.count >= 3 and any(k in str(tif_file) for k in ('ortho', 'topo'))
if is_rgb:
data = src.read([1, 2, 3])
data = np.moveaxis(data, 0, -1)
else:
data = src.read(1)
# Apply colormap normalization
data, cmap_name, title, legend_label, description, is_rgb_result = _apply_colormap(data, tif_file, resolution=resolution)
if not is_rgb_result:
data = np.where(np.isnan(data), 0.0, data)
# Convert to RGB using colormap
if is_rgb_result:
# RGB images are already in RGB (uint8 depuis le TIF IGN, ou float 0-1)
if data.dtype == np.uint8:
rgb_data = data
else:
rgb_data = (np.clip(data, 0, 1) * 255).astype(np.uint8)
else:
# Normalize data to 0-1 range for colormap
cmap = plt.get_cmap(cmap_name)
rgb_float = cmap(data.clip(0, 1))
rgb_data = (rgb_float[:, :, :3] * 255).astype(np.uint8)
# Save as AVIF/WebP
img = PILImage.fromarray(rgb_data)
pil_format = 'AVIF' if output_format == 'avif' else 'WEBP'
if quality >= 100:
img.save(str(output_file), format=pil_format, lossless=True)
else:
img.save(str(output_file), format=pil_format, quality=quality)
# Delete source TIFF (unless --keep-tif)
if not keep_tif:
tif_file.unlink(missing_ok=True)
return output_file
except Exception as e:
logger.error(f" Erreur conversion crop {ext.upper()}: {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 images in per-file subdirectory first, then fallback to main dir
file_vis_dir = vis_dir / basename
png_files = []
if file_vis_dir.exists():
png_files = sorted(file_vis_dir.glob("*.avif")) + sorted(file_vis_dir.glob("*.webp"))
else:
png_files = sorted(vis_dir.glob(f"{basename}_*.avif")) + sorted(vis_dir.glob(f"{basename}_*.webp"))
# Deduplicate in case both formats exist
seen = set()
unique_files = []
for f in png_files:
if f not in seen:
seen.add(f)
unique_files.append(f)
if not unique_files:
logger.warning(f" ✗ Aucune image trouvée pour {basename}")
return None
png_files = unique_files
# 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', 'hillshade_multi',
'flow_acc', 'solar', 'slope', 'roughness', 'wavelet']
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