[EJP] add sources

Author: FaugereE

bluemath_tk/core/plotting/WMS_USGSTopo.py

@@ -163,58 +163,57 @@ def get_cartopy_scale(zoom: int) -> str:
         return "110m"
 
 
-def webmercator_to_lonlat(x: float, y: float) -> Tuple[float, float]:
-    """
-    Converts Web Mercator projection coordinates (meters) to lon/lat.
-
-    Parameters
-    ----------
-    x : float
-        X coordinate in meters (Web Mercator).
-    y : float
-        Y coordinate in meters (Web Mercator).
-
-    Returns
-    -------
-    lon : float
-        Longitude in degrees.
-    lat : float
-        Latitude in degrees.
-    """
-    lon = math.degrees(x / EARTH_RADIUS_M)
-    lat = math.degrees(2 * math.atan(math.exp(y / EARTH_RADIUS_M)) - math.pi / 2)
-    return lon, lat
-
-
 def plot_usgs_raster_map(
     lat_min: float,
     lat_max: float,
     lon_min: float,
     lon_max: float,
     zoom: int = None,
+    source: str = "arcgis",
     verbose: bool = True,
     mask_ocean: bool = False,
     add_features: bool = True,
     display_width_px: int = 1024,
-    figsize: Tuple[int, int] = (12, 8),
-    fig: plt.Figure = None,
-    ax: plt.Axes = None,
 ) -> Tuple[plt.Figure, plt.Axes]:
     """
     Downloads and displays a USGS raster map for the given bounding box.
 
     Parameters
     ----------
     lat_min : float
-        Minimum latitude of the region.
+        Minimum latitude of the bounding box.
     lat_max : float
-        Maximum latitude of the region.
+        Maximum latitude of the bounding box.
     lon_min : float
-        Minimum longitude of the region.
+        Minimum longitude of the bounding box.
     lon_max : float
-        Maximum longitude of the region.
+        Maximum longitude of the bounding box.
+    zoom : int, optional
+        Zoom level (default is None, which auto-calculates based on bounding box).
+    source : str, optional
+        Source of the map tiles (default is "arcgis"). Options include:
+        - "arcgis" (Esri World Imagery)
+        - "google" (Google Maps Satellite)
+        - "eox" (EOX Sentinel-2 Cloudless)
+        - "osm" (OpenStreetMap Standard Tiles)
+        - "amazon" (Amazon Terrarium Elevation Tiles)
+        - "esri_world" (Esri World Imagery)
+        - "geoportail_fr" (Geoportail France Orthophotos)
+        - "ign_spain_pnoa" (IGN Spain PNOA Orthophotos)
+    verbose : bool, optional
+        If True, prints additional information about the source and zoom level.
+    mask_ocean : bool, optional
+        If True, masks ocean areas in the map.
+    add_features : bool, optional
+        If True, adds borders, coastlines, and states to the map.
     display_width_px : int, optional
-        Approximate pixel width for display (default is 1024).
+        Width of the display in pixels (default is 1024).
+    Returns
+    -------
+    fig : plt.Figure
+        The matplotlib figure containing the map.
+    ax : plt.Axes
+        The matplotlib axes with the map projection.
     """
 
     tile_size = 256
@@ -229,11 +228,117 @@ def plot_usgs_raster_map(
     height = y_end - y_start + 1
 
     map_img = Image.new("RGB", (width * tile_size, height * tile_size))
-    tile_url = "https://basemap.nationalmap.gov/arcgis/rest/services/USGSImageryOnly/MapServer/tile/{z}/{y}/{x}"
+    if source == "arcgis":
+        tile_url = "https://services.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}"
+        z_max = 19
+        if verbose:
+            print(
+                "Using Esri World Imagery (ArcGIS):\n"
+                "- Free for public and non-commercial use\n"
+                "- Commercial use allowed under Esri terms of service\n"
+                "- Attribution required: 'Tiles © Esri — Sources: Esri, Earthstar Geographics, CNES/Airbus DS, "
+                "USDA, USGS, AeroGRID, IGN, and the GIS User Community'\n"
+                "- Max zoom ~19"
+            )
+
+    elif source == "google":
+        tile_url = "https://mt1.google.com/vt/lyrs=s&x={x}&y={y}&z={z}"
+        z_max = 21
+        if verbose:
+            print(
+                "Using Google Maps Satellite:\n"
+                "- NOT license-free\n"
+                "- Usage outside official Google Maps SDKs/APIs is prohibited\n"
+                "- Commercial use requires a paid license from Google\n"
+                "- May be blocked without an API key\n"
+                "- Max zoom ~21"
+            )
+
+    elif source == "eox":
+        tile_url = "https://tiles.maps.eox.at/wmts/1.0.0/s2cloudless-2024_3857/default/g/{z}/{y}/{x}.jpg"
+        z_max = 16
+        if verbose:
+            print(
+                "Using EOX Sentinel-2 Cloudless:\n"
+                "- Based on Copernicus Sentinel-2 data processed by EOX\n"
+                "- Licensed under Creative Commons BY-NC-SA 4.0 (non-commercial use only)\n"
+                "- Attribution required: 'Sentinel-2 cloudless – © EOX, based on modified Copernicus Sentinel data 2016–2024'\n"
+                "- Versions from 2016–2017 are under CC BY 4.0 (commercial use allowed)\n"
+                "- Max zoom ~16"
+            )
+
+    elif source == "osm":
+        tile_url = "https://tile.openstreetmap.org/{z}/{x}/{y}.png"
+        z_max = 19
+        if verbose:
+            print(
+                "Using OpenStreetMap Standard Tiles:\n"
+                "- Free and open-source (ODbL license)\n"
+                "- Commercial use allowed with attribution\n"
+                "- Attribution required: '© OpenStreetMap contributors'\n"
+                "- Max zoom ~19"
+            )
+    
+    elif source == "amazon":
+        tile_url = "https://s3.amazonaws.com/elevation-tiles-prod/terrarium/{z}/{x}/{y}.png"
+        z_max = 15
+        if verbose:
+            print(
+                "Using Amazon Terrarium Elevation Tiles:\n"
+                "- Free for public use\n"
+                "- Attribution recommended: 'Amazon Terrarium'\n"
+                "- Max zoom ~15"
+            )
+
+    elif source == "esri_world":
+        tile_url = "https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}"
+        z_max = 19
+        if verbose:
+            print(
+                "Using Esri World Imagery:\n"
+                "- High-resolution global imagery\n"
+                "- Free with attribution under Esri terms\n"
+                "- Max zoom ~19"
+            )
+
+    elif source == "geoportail_fr":
+        tile_url = ("https://data.geopf.fr/wmts?"
+                    "REQUEST=GetTile&SERVICE=WMTS&VERSION=1.0.0"
+                    "&STYLE=normal&TILEMATRIXSET=PM&FORMAT=image/jpeg"
+                    "&LAYER=ORTHOIMAGERY.ORTHOPHOTOS&TILEMATRIX={z}"
+                    "&TILEROW={y}&TILECOL={x}")
+        z_max = 19
+        if verbose:
+            print(
+                "Using Geoportail France (Orthophotos):\n"
+                "- Aerial orthophotos from the French National Institute (IGN)\n"
+                "- Free for public use under Etalab license\n"
+                "- Attribution: Geoportail France / IGN\n"
+                "- Max zoom ~19"
+            )
+    elif source == "ign_spain_pnoa":
+        tile_url = ("https://tms-pnoa-ma.idee.es/1.0.0/pnoa-ma/{z}/{x}/{y_inv}.jpeg")
+        z_max = 19
+        if verbose:
+            print(
+                "Using IGN Spain PNOA Orthophotos:\n"
+                "- High-resolution aerial imagery (PNOA program)\n"
+                "- Provided by IGN/CNIG (Government of Spain)\n"
+                "- Free to use under Creative Commons BY 4.0 license\n"
+                "- Attribution required: 'Ortofotografía PNOA – © IGN / CNIG (Gobierno de España) – CC BY 4.0'\n"
+                "- Max zoom ~19"
+            )
+    
+    if zoom > z_max:
+        zoom = z_max
 
     for x in range(x_start, x_end + 1):
         for y in range(y_start, y_end + 1):
-            url = tile_url.format(z=zoom, x=x, y=y)
+            if source == "ign_spain_pnoa":
+                y_inv = (2**zoom - 1) - y
+                url = tile_url.format(z=zoom, x=x, y_inv=y_inv)
+            else:
+                url = tile_url.format(z=zoom, x=x, y=y)
             try:
                 response = requests.get(url, timeout=10)
                 tile = Image.open(BytesIO(response.content))
@@ -244,30 +349,21 @@ def plot_usgs_raster_map(
                 print(f"Error fetching tile {x},{y}: {e}")
 
     xmin, ymin, xmax, ymax = tile_bounds_meters(x_start, y_start, x_end, y_end, zoom)
-    lon_min, lat_min = webmercator_to_lonlat(xmin, ymin)
-    lon_max, lat_max = webmercator_to_lonlat(xmax, ymax)
-
-    scale = get_cartopy_scale(zoom)
-    if verbose:
-        print(f"Using Cartopy scale: {scale}")
 
     crs_tiles = ccrs.PlateCarree()
-    if fig is None or ax is None:
-        fig = plt.figure(figsize=figsize)
-        ax = plt.axes(projection=crs_tiles)
-
+    fig = plt.figure(figsize=(12, 8))
+    ax = plt.axes(projection=crs_tiles)
     ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=crs_tiles)
-
     ax.imshow(
-        map_img,
-        origin="upper",
-        extent=[xmin, xmax, ymin, ymax],
-        transform=ccrs.Mercator.GOOGLE,
+        map_img, origin="upper", extent=[xmin, xmax, ymin, ymax], transform=ccrs.Mercator.GOOGLE
     )
+    scale = get_cartopy_scale(zoom)
+    if verbose:
+        print(f"Using Cartopy scale: {scale}")
 
     if add_features:
         ax.add_feature(cfeature.BORDERS.with_scale(scale), linewidth=0.8)
-        ax.add_feature(cfeature.COASTLINE.with_scale(scale), zorder=1)
+        ax.add_feature(cfeature.COASTLINE.with_scale(scale))
         ax.add_feature(cfeature.STATES.with_scale(scale), linewidth=0.5)
 
     if mask_ocean: