From 93d9152d415cfa8dddda5a8b3387308cf1eae981 Mon Sep 17 00:00:00 2001 From: akrherz Date: Fri, 3 Feb 2023 08:18:02 -0600 Subject: [PATCH] add ESRI basemap backgrounds refs #304 --- CHANGELOG.md | 2 + src/pyiem/data/geodf/SLC.parquet | Bin 0 -> 8295 bytes src/pyiem/plot/_mpl.py | 187 ++++++++++++++++++++++++++++++- src/pyiem/plot/util.py | 27 ++++- 4 files changed, 211 insertions(+), 5 deletions(-) create mode 100644 src/pyiem/data/geodf/SLC.parquet diff --git a/CHANGELOG.md b/CHANGELOG.md index ca8231967..beccebcaf 100644 --- a/CHANGELOG.md +++ b/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. diff --git a/src/pyiem/data/geodf/SLC.parquet b/src/pyiem/data/geodf/SLC.parquet new file mode 100644 index 0000000000000000000000000000000000000000..751a80294319c184d37b5f2b51aa2b263f54d33d GIT binary patch literal 8295 zcmeHN32+-{6_yj@q>!rxr*`U6RkR_Yj=a0l%2r6|N?KixwT@+3S5(}+EUmQKRi)LL z21+PQFPIe4a!iw6p+IQqka8!G1qNsdZ5f~~SIRvtG(Cm}7(xTQe^<8U5ROcTfim%o zrTwq>fA78TdvD{(R!o9*JZJEEw;c1{H}rDds?5L znCQ#sAH4O6-ajX3s#dso{h!az=+x9ByB~ZjO+^ z^3b`i3s)UKzw63}_djyZH8geT>6>pQFUjb!RqOXY_fDF+?8w1uuDUd%A1aoI@f&Dr z&F%N>xO`nkAJBU*Ty#B6{eV9)@wqnvvIidAliNvCjj4k|crc@nJgg7yZPC<{_p9C; zhBEqt{*~1y-%eAa`n=d zJn}P|S`~ga_1i#3-*wxHi#I<6_;#}AcENb}GfT8{K0;F??>m3VT^lmGe)Pnhhw3!- z^zg5lcfnY=``Sf+m1t_)3G6Bw{25&HwQsCvY3ePzy1MstZT{3H{g9T-cYgGdmzKHj zyn1itjyJz?Nk+f-@bi}+6lm(hj?X;$(?4y~Gei6SunxRF7dxKlI`&}K^55^Fsq_n< zFE-E3==9zB@7=kVrtsyq@r{c!dXf&kL%yG;p8QyM;nh7Ez3@H& zdF=49PaZfzQ%c~-T~iat_l~T4Eu`ehKN7W@VZ86_ zcYSLmkRe+=`M&qi)RB!ZeI)_o%s{u}n(Jul3oD-b@}=j2v(Iz8e|IfVBK4-%cCP~b z*KP@JhLnc~zx4gz!1!2|*muv>H1);5{4xC$oYybgwPGcN_{f3oT_Y#6WLvB;v8WtT3R*Q8~2cmF{D zz*)imf4qC173?csu&6_JuU!gV@^uzV|E9j45~~PorSCI~x7#E|5F69HDp!k2#kSG9 zZM&@ibMRyL=JVA`8DjEilv z3$4sZqAW}=oHJofYt0h6&Q?_twT{iA5{!9&uXd_Jf5YoP^-PQQ&QxO#kfcx)_kUB5CVUg#W!0!C?n> zCfi4b9cQu~$H_A?kGqDQXJmiaMV!fg0(ZD!KN7Jh%TlSLN~jog#A+^mkTz+$$mZ4Z zG~19+5?qeq5!^8hXC(#kR-B%ZSxplaQaI1mPfE4n%R5=c{FwH}pW*F~o`X zZERJ`sbxtsmi${MgJ{jD2Iv;TwF4dJRkaL)XEkxUVi3dtTB&lg9GvrjmC^Y&$_}c` z5e#5TTGcEzEXS3_nglbZCdXK#)(BD=s!n@OQOjBm*@Pm6ilI#qvt?0SNUM44s*?T3 z=1tu!FfXzdcs3WP#lS2I62}bN!K8N#c{hE za}p!45ht7=7}$r~he;RVb~{NoIqV`wt}qzXQCu;*(_|oj`g#oz&5Dd-y5Jm6*a>Jg zmcC`4XmmWW(JY94-|ihkKMhaV!RTq9v}#OXHG9ies8D zEY-8&BGcrXU=#L?hL}7XjBBY$;|^5_MeBS^t&KUnO^#5joMTiObF}v)BQ+uKsWOvp z!05h;RREtvwo*h+<`Vv@IX19iV$v%kT;T%vMhMMR(yOuwhMtlm(wNU9am9d^Ht$SE z@kn-@o-pRq!2l*C@EO3~N+)oM3r+!c?{K(C=C}lU20sXQlM7B{qeQ%x4#qJy5y?(D z}9nj9$-j+bLnF35~a-tMgunLM@rwQL zqf)BrjfqJHODC{wBvx=TK93r0@vVrjAk9qHQNF6F1a8Ji7HWL4Q3(4{tfI0Il=I;L zL-LNO%uRY4nPsBAg?$%O?QjM40IWZ_3nql%0$2je8F->Go zt_p}BM0_S)Y~(-(G6!>Mq!T(!j{BTVCdM}-`3TPB8OVW0N+p=QkQ~jW90h}y5y@FM z^(7bDn#>8ws0!RD2_bV&H13T>o3r|j8u6*lQ@afDn*CF+KLQpsG&(Do@D@6-wh8e>6qX|QAZF)H!JIKd`hKH;vK zx`i0{44efvIVz_dfQfJ#esV>3QwU~_vqFL_!9G-<3h9K(EkyD^6V@-)yL5~I|r~?Y@ zlfw!wp!>rLljj_>H7Ex$m0cV`Vi2cnig0T@*oPW{g@f*5m}dOp(5&8^?KpzZ?fPVD zyFFhVt%+^E0rzkrn49U0C)6ixJmx`uInsbLNi6Ue;)FHEX=5GLR*Vg}E6DG5Jb-J3 zI?72p5oigp2AGxET!C#k>os`7CXnB#4g?Se;YB70#lc7@%##hfWw7<%iK4q_N3s+TZejINF||7k4?cg6l8D;5swDPVp(bC6zr#~M<`cZB95WlAwSbaS>}S#nT{IX`96bq z+s+3wBZTC7``tkjlTjJ$qxCbxA0%-u*l72Xf6gmDBtyeASUw`n7;lxvtKd|dZ}6Q(@0ak-XuL;GAX1t?e!sEPXN-(Ta~|>4J0BnO zO%>KVj(R+sEq&G!TfV9)Y8&`cX@l{%3v2r|*lPaOqJLlj{yO~n)fUTE_%8;31AsH( A{r~^~ literal 0 HcmV?d00001 diff --git a/src/pyiem/plot/_mpl.py b/src/pyiem/plot/_mpl.py index 40046a192..7668dbf01 100644 --- a/src/pyiem/plot/_mpl.py +++ b/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) diff --git a/src/pyiem/plot/util.py b/src/pyiem/plot/util.py index a8090571e..5c5cf2a0d 100644 --- a/src/pyiem/plot/util.py +++ b/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,