"""IGN (French National Geographic Institute) tile download and overlay generation. Downloads WMTS tiles from the IGN Geoportail and composites them into GeoTIFF overlays matching the LiDAR DTM extent. """ import logging import math import time from pathlib import Path import numpy as np import rasterio try: from rasterio.warp import transform as warp_transform HAS_WARP = True except ImportError: HAS_WARP = False logger = logging.getLogger("lidar") def _lat_lon_to_tile(lat, lon, zoom): """Convert lat/lon to Web Mercator tile coordinates.""" n = 2 ** zoom col = int((lon + 180) / 360 * n) lat_rad = math.radians(lat) row = int((1 - math.log(math.tan(lat_rad) + 1 / math.cos(lat_rad)) / math.pi) / 2 * n) return col, row def _lat_lon_to_px(lat, lon, zoom, tile_size=256): """Convert lat/lon to Web Mercator pixel coordinates.""" n = 2 ** zoom px_x = (lon + 180) / 360 * n * tile_size lat_rad = math.radians(lat) px_y = (1 - math.log(math.tan(lat_rad) + 1 / math.cos(lat_rad)) / math.pi) / 2 * n * tile_size return px_x, px_y def download_ign_tiles(min_x, max_x, min_y, max_y, layer, zoom_level=15): """Download IGN WMTS tiles for the given bounds using Web Mercator (PM). Args: min_x, max_x, min_y, max_y: Bounds in Lambert 93. layer: IGN WMTS layer name. zoom_level: WMTS zoom level (default 15). Returns: numpy array (H, W, 3) uint8, or None on failure. """ import urllib.request import io from PIL import Image as PILImage if not HAS_WARP: logger.error(" ✗ rasterio.warp non disponible pour conversion de coordonnées") return None try: l93_xs = [min_x, max_x] l93_ys = [max_y, min_y] lons, lats = warp_transform('EPSG:2154', 'EPSG:4326', l93_xs, l93_ys) nw_lat, nw_lon = lats[0], lons[0] se_lat, se_lon = lats[1], lons[1] except Exception: logger.error(" ✗ Impossible de convertir les coordonnées") return None wmts_url = "https://data.geopf.fr/wmts" tile_matrix_set = "PM" tile_size = 256 col_min, row_min = _lat_lon_to_tile(nw_lat, nw_lon, zoom_level) col_max, row_max = _lat_lon_to_tile(se_lat, se_lon, zoom_level) nw_px_x, nw_px_y = _lat_lon_to_px(nw_lat, nw_lon, zoom_level) se_px_x, se_px_y = _lat_lon_to_px(se_lat, se_lon, zoom_level) out_width = int(se_px_x - nw_px_x) out_height = int(se_px_y - nw_px_y) if out_width <= 0 or out_height <= 0 or out_width > 5000 or out_height > 5000: return None composite = np.full((out_height, out_width, 3), 255, dtype=np.uint8) tiles_downloaded = 0 fmt = "image/png" if 'PLAN' in layer else "image/jpeg" for col in range(col_min, col_max + 1): for row in range(row_min, row_max + 1): url = ( f"{wmts_url}?SERVICE=WMTS&VERSION=1.0.0&REQUEST=GetTile" f"&LAYER={layer}&STYLE=normal" f"&TILEMATRIXSET={tile_matrix_set}" f"&TILEMATRIX={zoom_level}&TILECOL={col}&TILEROW={row}" f"&FORMAT={fmt}" ) try: req = urllib.request.Request(url, headers={'User-Agent': 'Mozilla/5.0'}) with urllib.request.urlopen(req, timeout=10) as response: tile_data = response.read() tile_img = PILImage.open(io.BytesIO(tile_data)).convert('RGB') tile_arr = np.array(tile_img) tile_origin_x = col * tile_size tile_origin_y = row * tile_size px_x = int(tile_origin_x - nw_px_x) px_y = int(tile_origin_y - nw_px_y) x_off = max(0, -px_x) y_off = max(0, -px_y) dst_x_start = max(0, px_x) dst_y_start = max(0, px_y) dst_x_end = min(out_width, px_x + tile_size) dst_y_end = min(out_height, px_y + tile_size) src_x = x_off src_y = y_off src_w = dst_x_end - dst_x_start src_h = dst_y_end - dst_y_start if src_w > 0 and src_h > 0 and tile_arr.shape[0] >= src_y + src_h and tile_arr.shape[1] >= src_x + src_w: composite[dst_y_start:dst_y_end, dst_x_start:dst_x_end] = \ tile_arr[src_y:src_y+src_h, src_x:src_x+src_w] tiles_downloaded += 1 except Exception as e: if tiles_downloaded == 0 and col == col_min and row == row_min: logger.error(f" ✗ Erreur tuile IGN: {e}") continue logger.info(f" → {tiles_downloaded} tuiles IGN téléchargées ({layer})") if tiles_downloaded == 0: return None return composite def generate_ign_overlay(dem_file, basename, vis_dir, resolution, layer, title, legend_label, description, out_suffix): """Generate an IGN basemap overlay visualization. Args: dem_file: Path to DTM GeoTIFF. basename: Base name for output file. vis_dir: Output directory for visualizations. resolution: Grid resolution in m/px. layer: IGN WMTS layer name. title: Title for the overlay. legend_label: Legend text. description: Description text. out_suffix: Suffix for output filename. Returns: Path to output GeoTIFF, or None on failure. """ logger.info(f" → {title}...") t0 = time.time() output = vis_dir / f"{basename}_{out_suffix}.tif" try: with rasterio.open(dem_file) as src: dem = src.read(1) transform = src.transform crs = src.crs height, width = dem.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 zoom = 15 result = download_ign_tiles(min_x, max_x, min_y, max_y, layer, zoom) if result is None: logger.error(" ✗ Impossible de télécharger les tuiles IGN") return None from PIL import Image as PILImage ign_pil = PILImage.fromarray(result) ign_resized = ign_pil.resize((width, height), PILImage.LANCZOS) ign_arr = np.array(ign_resized) with rasterio.open( output, 'w', driver='GTiff', height=height, width=width, count=3, dtype='uint8', crs=crs, transform=transform, ) as dst: for i in range(3): dst.write(ign_arr[:, :, i], i + 1) logger.info(f" ✓ {title} terminé ({time.time()-t0:.1f}s)") return output except Exception as e: logger.error(f" ✗ Erreur IGN overlay: {e}", exc_info=True) return None