From 4836ea297d235b2b9253edde084dd6cc1579446c Mon Sep 17 00:00:00 2001 From: FaugereE Date: Tue, 8 Jul 2025 11:00:48 +0200 Subject: [PATCH 1/4] [EJP] update to ccrs.PlateCarree() --- bluemath_tk/core/plotting/WMS_USGSTopo.py | 31 ++++++++++++++++++++--- 1 file changed, 28 insertions(+), 3 deletions(-) diff --git a/bluemath_tk/core/plotting/WMS_USGSTopo.py b/bluemath_tk/core/plotting/WMS_USGSTopo.py index a97818e..3f61e67 100644 --- a/bluemath_tk/core/plotting/WMS_USGSTopo.py +++ b/bluemath_tk/core/plotting/WMS_USGSTopo.py @@ -162,6 +162,27 @@ def get_cartopy_scale(zoom: int) -> str: else: 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, @@ -218,14 +239,18 @@ 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) + + print(f"Converted bounds to degrees: lon_min={lon_min}, lon_max={lon_max}, lat_min={lat_min}, lat_max={lat_max}") - 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: From 724b9947e390e3ed30c3f52ea89d86cde9d7395a Mon Sep 17 00:00:00 2001 From: FaugereE Date: Tue, 8 Jul 2025 12:11:07 +0200 Subject: [PATCH 2/4] [EJP] update to ccrs.PlateCarree( --- bluemath_tk/core/plotting/WMS_USGSTopo.py | 27 +++++++++++++++-------- 1 file changed, 18 insertions(+), 9 deletions(-) diff --git a/bluemath_tk/core/plotting/WMS_USGSTopo.py b/bluemath_tk/core/plotting/WMS_USGSTopo.py index 3f61e67..34b42b2 100644 --- a/bluemath_tk/core/plotting/WMS_USGSTopo.py +++ b/bluemath_tk/core/plotting/WMS_USGSTopo.py @@ -162,6 +162,7 @@ def get_cartopy_scale(zoom: int) -> str: else: return "110m" + def webmercator_to_lonlat(x: float, y: float) -> Tuple[float, float]: """ Converts Web Mercator projection coordinates (meters) to lon/lat. @@ -184,6 +185,7 @@ def webmercator_to_lonlat(x: float, y: float) -> Tuple[float, float]: 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, @@ -194,6 +196,9 @@ def plot_usgs_raster_map( 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. @@ -242,27 +247,31 @@ def plot_usgs_raster_map( lon_min, lat_min = webmercator_to_lonlat(xmin, ymin) lon_max, lat_max = webmercator_to_lonlat(xmax, ymax) - print(f"Converted bounds to degrees: lon_min={lon_min}, lon_max={lon_max}, lat_min={lat_min}, lat_max={lat_max}") + scale = get_cartopy_scale(zoom) + if verbose: + print(f"Using Cartopy scale: {scale}") crs_tiles = ccrs.PlateCarree() - fig = plt.figure(figsize=(12, 8)) - ax = plt.axes(projection=crs_tiles) + if fig is None or ax is None: + fig = plt.figure(figsize=figsize) + 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)) + ax.add_feature(cfeature.COASTLINE.with_scale(scale), zorder=1) 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 From 7e6e430f491eeaad57615c3dc8c5b2cfd8fd6786 Mon Sep 17 00:00:00 2001 From: FaugereE Date: Tue, 8 Jul 2025 15:25:33 +0200 Subject: [PATCH 3/4] [EJP] add sources --- bluemath_tk/core/plotting/WMS_USGSTopo.py | 194 ++++++++++++++++------ 1 file changed, 145 insertions(+), 49 deletions(-) diff --git a/bluemath_tk/core/plotting/WMS_USGSTopo.py b/bluemath_tk/core/plotting/WMS_USGSTopo.py index 34b42b2..e26de43 100644 --- a/bluemath_tk/core/plotting/WMS_USGSTopo.py +++ b/bluemath_tk/core/plotting/WMS_USGSTopo.py @@ -163,42 +163,17 @@ 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. @@ -206,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 @@ -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: From 67bb5e59af034668adb7c3168bd52951c3992830 Mon Sep 17 00:00:00 2001 From: FaugereE Date: Tue, 8 Jul 2025 15:28:32 +0200 Subject: [PATCH 4/4] [EJP] add sources --- bluemath_tk/core/plotting/WMS_USGSTopo.py | 1 + 1 file changed, 1 insertion(+) diff --git a/bluemath_tk/core/plotting/WMS_USGSTopo.py b/bluemath_tk/core/plotting/WMS_USGSTopo.py index e26de43..b075e2e 100644 --- a/bluemath_tk/core/plotting/WMS_USGSTopo.py +++ b/bluemath_tk/core/plotting/WMS_USGSTopo.py @@ -386,3 +386,4 @@ def plot_usgs_raster_map( mask_ocean=True, ) fig.show() +