From bc92689accac7dfa6c9345d856e8bc98c600b0d3 Mon Sep 17 00:00:00 2001 From: Jacquin Antoine Date: Sun, 10 May 2026 03:22:08 +0200 Subject: [PATCH] =?UTF-8?q?Orthophoto=20IGN:=20zoom=20adaptatif=20pour=20r?= =?UTF-8?q?=C3=A9solution=20optimale?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Zoom WMTS calculé automatiquement à partir de la résolution DTM (0.5m/px → zoom 19, 0.2m/px → zoom 20) au lieu de zoom 15 fixe - Orthophoto: zoom max 20 (~0.15m/px), carte topo: zoom max 19 - Limite image 10000x10000px (au lieu de 5000x5000) - Log du nombre de tuiles et de la taille avant téléchargement Co-Authored-By: Claude Opus 4.6 --- lidar_pipeline/ign.py | 51 +++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 49 insertions(+), 2 deletions(-) diff --git a/lidar_pipeline/ign.py b/lidar_pipeline/ign.py index f401cc0..dd4b128 100644 --- a/lidar_pipeline/ign.py +++ b/lidar_pipeline/ign.py @@ -21,6 +21,40 @@ except ImportError: logger = logging.getLogger("lidar") +def _optimal_zoom_level(resolution, lat, layer): + """Compute optimal WMTS zoom level for a given ground resolution. + + Chooses the zoom level where the WMTS ground resolution is at least + as fine as the requested resolution, capped to avoid downloading + too many tiles. + + Args: + resolution: Desired ground resolution in m/px (DTM resolution). + lat: Latitude in degrees (for Web Mercator resolution calculation). + layer: IGN WMTS layer name (different max zoom per layer). + + Returns: + Integer zoom level. + """ + # Web Mercator ground resolution: m/px = (2π * R * cos(lat)) / (256 * 2^z) + # Solve for z: z = log2(2π * R * cos(lat) / (256 * resolution)) + R = 6378137.0 # WGS84 equatorial radius + lat_rad = math.radians(lat) + numerator = 2 * math.pi * R * math.cos(lat_rad) + z_exact = math.log2(numerator / (256 * resolution)) + + zoom = max(15, min(20, int(math.ceil(z_exact)))) + + # Cap zoom per layer: orthophoto supports up to 20, topographic map up to 19 + if 'ORTHO' in layer: + zoom = min(zoom, 20) + elif 'PLAN' in layer: + zoom = min(zoom, 19) + + logger.debug(f" Zoom WMTS: {zoom} (résolution cible: {resolution}m/px)") + return zoom + + def _lat_lon_to_tile(lat, lon, zoom): """Convert lat/lon to Web Mercator tile coordinates.""" n = 2 ** zoom @@ -81,9 +115,13 @@ def download_ign_tiles(min_x, max_x, min_y, max_y, layer, zoom_level=15): 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: + if out_width <= 0 or out_height <= 0 or out_width > 10000 or out_height > 10000: + logger.warning(f" Image IGN trop grande: {out_width}x{out_height}px — abandon") return None + total_tiles = (col_max - col_min + 1) * (row_max - row_min + 1) + logger.info(f" Zoom {zoom_level}: {total_tiles} tuiles à télécharger ({out_width}x{out_height}px)") + composite = np.full((out_height, out_width, 3), 255, dtype=np.uint8) tiles_downloaded = 0 @@ -179,7 +217,16 @@ def generate_ign_overlay(dem_file, basename, vis_dir, resolution, layer, title, max_y = top_left_y min_y = top_left_y - height * pixel_size_y - zoom = 15 + # Compute center latitude for zoom calculation + center_lat = (max_y + min_y) / 2 + center_lon = (min_x + max_x) / 2 + try: + clons, clats = warp_transform('EPSG:2154', 'EPSG:4326', [center_lon], [center_lat]) + center_lat_wgs84 = clats[0] + except Exception: + center_lat_wgs84 = 46.0 # fallback: center of France + + zoom = _optimal_zoom_level(resolution, center_lat_wgs84, layer) result = download_ign_tiles(min_x, max_x, min_y, max_y, layer, zoom)