diff --git a/bluemath_tk/core/plotting/WMS_USGSTopo.py b/bluemath_tk/core/plotting/WMS_USGSTopo.py index a97818e..b075e2e 100644 --- a/bluemath_tk/core/plotting/WMS_USGSTopo.py +++ b/bluemath_tk/core/plotting/WMS_USGSTopo.py @@ -169,6 +169,7 @@ def plot_usgs_raster_map( lon_min: float, lon_max: float, zoom: int = None, + source: str = "arcgis", verbose: bool = True, mask_ocean: bool = False, add_features: bool = True, @@ -180,15 +181,39 @@ def plot_usgs_raster_map( 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 @@ -203,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)) @@ -219,13 +350,12 @@ def plot_usgs_raster_map( xmin, ymin, xmax, ymax = tile_bounds_meters(x_start, y_start, x_end, y_end, zoom) - crs_tiles = ccrs.Mercator.GOOGLE + crs_tiles = ccrs.PlateCarree() fig = plt.figure(figsize=(12, 8)) ax = plt.axes(projection=crs_tiles) - ax.set_extent([xmin, xmax, ymin, ymax], crs=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=crs_tiles + map_img, origin="upper", extent=[xmin, xmax, ymin, ymax], transform=ccrs.Mercator.GOOGLE ) scale = get_cartopy_scale(zoom) if verbose: @@ -237,7 +367,7 @@ def plot_usgs_raster_map( ax.add_feature(cfeature.STATES.with_scale(scale), linewidth=0.5) if mask_ocean: - ax.add_feature(cfeature.OCEAN.with_scale(scale), facecolor="w", zorder=3) + ax.add_feature(cfeature.OCEAN.with_scale(scale), facecolor="w", zorder=1) return fig, ax @@ -256,3 +386,4 @@ def plot_usgs_raster_map( mask_ocean=True, ) fig.show() +