Orthophoto IGN: zoom adaptatif pour résolution optimale

- 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 <noreply@anthropic.com>

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)