add ESRI basemap backgrounds

refs #304

CHANGELOG.md

@@ -30,6 +30,8 @@ database has.
 - Generate a TextProduct.warning message for a VTEC product that should contain
 a polygon, but does not (#660).
 - Introduce a natural earth background option for MapPlot (#304).
+- Introduce hacky `sector="spherical_mercator"` that brings in ESRI basemaps
+for the background.  My implementation stinks and will change (#304).
 - Support `cartopy_offlinedata` version 0.20+.
 - Support new CLI format diction from NWS Sacramento.
 - Workaround autoplot context fun with mixed 3-4 character WFOs.

src/pyiem/plot/_mpl.py

@@ -1,17 +1,135 @@
 """A custom axes implementation."""
 # pylint: disable=unsubscriptable-object,unpacking-non-sequence
+import math
 import os
+from io import BytesIO
 
 # Third party libraries
 import numpy as np
+from PIL import Image
 from shapely.geometry import Polygon
 import geopandas as gpd
+from pymemcache.client import Client
 from pyproj import Transformer
 import rasterio
+import requests
 from rasterio.warp import reproject, Resampling
 
 # Local imports
-from pyiem.reference import LATLON
+from pyiem.reference import LATLON, EPSG, Z_FILL
+from pyiem.util import LOG, load_geodf
+
+# Zoom 0 through 24
+METERS_PER_PIXEL = [
+    float(x)
+    for x in (
+        "156543 78271.5 39135.8 19567.88 9783.94 4891.97 2445.98 1222.99 "
+        "611.5 305.75 152.87 76.44 38.219 19.109 9.555 4.777 2.3887 1.1943 "
+        "0.5972 0.2986 0.14929 0.074646 0.037323 0.0186615 0.00933075"
+    ).split()
+]
+
+
+def get_tile_lat_lng(zoom, x, y):
+    """convert Google-style Mercator tile coordinate to
+    (lat, lng) of top-left corner of tile"""
+
+    # "map-centric" latitude, in radians:
+    lat_rad = math.pi - 2 * math.pi * y / (2**zoom)
+    # true latitude:
+    lat_rad = 2 * math.atan(math.exp(lat_rad)) - math.pi / 2
+    lat = lat_rad * 180.0 / math.pi
+
+    # longitude maps linearly to map, so we simply scale:
+    lng = -180.0 + 360.0 * x / (2**zoom)
+
+    return (lat, lng)
+
+
+def get_lat_lng_tile(lat, lng, zoom):
+    """convert lat/lng to Google-style Mercator tile coordinate (x, y)
+    at the given zoom level"""
+
+    lat_rad = lat * math.pi / 180.0
+    # "map-centric" latitude, in radians:
+    lat_rad = math.log(math.tan((lat_rad + math.pi / 2) / 2))
+
+    x = 2**zoom * (lng + 180.0) / 360.0
+    y = 2**zoom * (math.pi - lat_rad) / (2 * math.pi)
+
+    return (x, y)
+
+
+def get_tile_data(url):
+    """Fetch the tile and hope memcached has it."""
+    key = url.split("//")[1]
+    mc = Client("iem-memcached:11211", timeout=5)
+    res = mc.get(key)
+    if res is None:
+        LOG.info("Fetching %s", url)
+        req = requests.get(url, timeout=10)
+        res = req.content
+        mc.set(key, res)
+    mc.close()
+    bio = BytesIO(res)
+    bio.seek(0)
+    with Image.open(bio) as pilimg:
+        im = np.asarray(pilimg)
+    return im
+
+
+def draw_wmts(panel, background):
+    """todo"""
+    xmin, xmax, ymin, ymax = panel.get_extent(crs=EPSG[4326])
+    tx0, ty0 = get_lat_lng_tile(ymax, xmin, panel.zoom)
+    tx1, ty1 = get_lat_lng_tile(ymin, xmax, panel.zoom)
+    # resolution = 2445.98
+    transform = Transformer.from_crs(EPSG[4326], EPSG[3857], always_xy=True)
+    for y in range(int(ty0), int(ty1) + 1):
+        for x in range(int(tx0), int(tx1) + 1):
+            minlat, minlon = get_tile_lat_lng(panel.zoom, x, y + 1)
+            maxlat, maxlon = get_tile_lat_lng(panel.zoom, x + 1, y)
+
+            x0, y0 = transform.transform(minlon, maxlat)
+            x1, y1 = transform.transform(maxlon, minlat)
+            try:
+                im = get_tile_data(
+                    "https://services.arcgisonline.com/arcgis/rest/services/"
+                    f"{background}/MapServer/tile/{panel.zoom}/{y}/{x}"
+                )
+            except Exception as exp:
+                LOG.info(exp)
+                continue
+
+            panel.ax.imshow(
+                im / 255.0,
+                interpolation="nearest",  # prevents artifacts
+                extent=(x0, x1, y0, y1),
+                origin="lower",
+                zorder=Z_FILL,
+            ).set_rasterized(True)
+
+    panel.ax.annotate(
+        "Basemap Courtesy ESRI",
+        xy=(1, 0),
+        bbox=dict(color="white"),
+        ha="right",
+        va="bottom",
+        zorder=1000,
+        xycoords="axes fraction",
+    )
+
+    # Don't ask
+    (
+        load_geodf("SLC")
+        .to_crs(panel.crs)
+        .plot(
+            ax=panel.ax,
+            aspect=None,
+            color="#f4f2db",
+            zorder=Z_FILL,
+        )
+    )
 
 
 def draw_background(panel, background):
@@ -87,6 +205,7 @@ def __init__(self, fig, rect, crs, **kwargs):
         self.figure = fig
         self.ax = self.figure.add_axes(rect, **kwargs)
         self.crs = crs
+        self.zoom = None
 
     draw_background = draw_background
 
@@ -231,3 +350,69 @@ def pcolormesh(self, x, y, vals, **kwargs):
         transform = Transformer.from_crs(crs, self.crs, always_xy=True)
         x, y = transform.transform(x, y)
         return self.ax.pcolormesh(x, y, vals, **kwargs)
+
+
+class SphericalMercatorPanel(GeoPanel):
+    """Specialized panel that maintains aspect."""
+
+    def draw_background(self, background):
+        """call"""
+        if background is None:
+            return
+        draw_wmts(self, background)
+
+    def set_extent(self, extent, crs=None):
+        """
+        Set the extent of the axes.
+
+        Args:
+            extent: (xmin, xmax, ymin, ymax)
+            crs: (optional) the CRS of the extent. default is 4326.
+        """
+        transform = Transformer.from_crs(
+            LATLON if crs is None else crs,
+            self.crs,
+            always_xy=True,
+        )
+        # Get the bounds in our local coordinates
+        ll = transform.transform(extent[0], extent[2])
+        lr = transform.transform(extent[1], extent[2])
+        ur = transform.transform(extent[1], extent[3])
+        ul = transform.transform(extent[0], extent[3])
+        xmin = min(ll[0], ul[0])
+        xmax = max(lr[0], ur[0])
+        ymin = min(ll[1], lr[1])
+        ymax = max(ul[1], ur[1])
+        # rectify to square pixels
+        bbox = self.ax.get_window_extent().transformed(
+            self.figure.dpi_scale_trans.inverted(),
+        )
+        display_ratio = bbox.height / bbox.width
+        dx = xmax - xmin
+        dy = ymax - ymin
+        data_ratio = dy / dx
+        if display_ratio > data_ratio:
+            # We need to expand the data in y
+            yadd = dx * display_ratio - dy
+            ymin -= yadd / 2
+            ymax += yadd / 2
+        else:
+            # We need to expand the data in x
+            xadd = dy / display_ratio - dx
+            xmin -= xadd / 2
+            xmax += xadd / 2
+
+        # Now we muck to get an int zoom level
+        dx = (xmax - xmin) / (bbox.width * 100.0)
+        self.zoom = np.digitize(dx, METERS_PER_PIXEL) - 1
+
+        # Now we recompute the bounds
+        centerx = (xmax + xmin) / 2.0
+        centery = (ymax + ymin) / 2.0
+        xmin = centerx - METERS_PER_PIXEL[self.zoom] * (bbox.width * 50.0)
+        xmax = centerx + METERS_PER_PIXEL[self.zoom] * (bbox.width * 50.0)
+        ymin = centery - METERS_PER_PIXEL[self.zoom] * (bbox.height * 50.0)
+        ymax = centery + METERS_PER_PIXEL[self.zoom] * (bbox.height * 50.0)
+
+        self.ax.set_xlim(xmin, xmax)
+        self.ax.set_ylim(ymin, ymax)

src/pyiem/plot/util.py

@@ -15,7 +15,7 @@
 from pyiem import reference
 from pyiem.plot.colormaps import stretch_cmap
 from pyiem.reference import LATLON, FIGSIZES
-from ._mpl import GeoPanel
+from ._mpl import GeoPanel, SphericalMercatorPanel
 
 DATADIR = os.sep.join([os.path.dirname(__file__), "..", "data"])
 LOGO_BOUNDS = (0.005, 0.91, 0.08, 0.086)
@@ -237,7 +237,10 @@ def make_panel(
     """
     # https://jdhao.github.io/2017/06/03/change-aspect-ratio-in-mpl/
     # aspect is data ratio divided by axes ratio
-    gp = GeoPanel(
+    sector_label = kwargs.get("sector_label", "")
+    hacky = "spherical_mercator"
+    cls = SphericalMercatorPanel if sector_label == hacky else GeoPanel
+    gp = cls(
         fig,
         ndc_axbounds,
         crs,
@@ -246,7 +249,7 @@ def make_panel(
         facecolor=(0.4471, 0.6235, 0.8117),
         xticks=[],
         yticks=[],
-        sector_label=kwargs.get("sector_label", ""),
+        sector_label=sector_label,
     )
     # Get the frame at the proper zorder
     for _k, spine in gp.ax.spines.items():
@@ -327,7 +330,23 @@ def sector_setter(mp, axbounds, **kwargs):
             **kwargs,
         )
         mp.panels.append(gp)
-
+    elif mp.sector == "spherical_mercator":
+        # Crude check if we are crossing the anti-meridian
+        if kwargs["west"] > 90 and kwargs["east"] < -90:
+            # Convert to 0 to 360 longitude?
+            kwargs["east"] += 360
+        gp = make_panel(
+            axbounds,
+            mp.fig,
+            [kwargs["west"], kwargs["east"], kwargs["south"], kwargs["north"]],
+            reference.EPSG[3857],
+            "auto",
+            is_geoextent="projection" not in kwargs,
+            sector_label=mp.sector,
+            background=kwargs.pop("background", "World_Street_Map"),
+            **kwargs,
+        )
+        mp.panels.append(gp)
     elif mp.sector == "north_america":
         gp = make_panel(
             axbounds,