"""pyiem.plot.util Plotting Utilities."""
# pylint: disable=import-outside-toplevel
import functools
import os
import geopandas as gpd
import matplotlib.colors as mpcolors
import matplotlib.image as mpimage
import matplotlib.patches as mpatches
import matplotlib.path as mpath
import numpy as np
import pandas as pd
from pyproj import Transformer
from shapely.geometry import MultiPolygon, Polygon
from pyiem import reference
from pyiem.plot.colormaps import stretch_cmap
from pyiem.reference import FIGSIZES, LATLON
from ._mpl import GeoPanel, SphericalMercatorPanel
DATADIR = os.sep.join([os.path.dirname(__file__), "..", "data"])
LOGO_BOUNDS = (0.005, 0.91, 0.08, 0.086)
LOGOFILES = {"dep": "deplogo.png", "iem": "logo.png", "nwa": "nwalogo.png"}
[docs]
def update_kwargs_apctx(func):
"""Decorate things provided by an autoplot context dict."""
@functools.wraps(func)
def wrapped(*args, **kwargs):
"""Compute things and update things."""
if "apctx" in kwargs:
apctx = kwargs["apctx"]
# If figsize is set, we take no other action
if "figsize" not in kwargs:
_r = apctx.get("_r", None)
if _r in FIGSIZES:
kwargs["figsize"] = FIGSIZES[_r]
# Merge in csector, this will override sector
csector = apctx.get("csector", None)
if csector is not None:
# Quasi magic
if len(csector) == 2:
kwargs["sector"] = "state"
kwargs["state"] = csector
else:
kwargs["sector"] = csector
# Merge in dpi, if not set
dpi = apctx.get("dpi", None)
if dpi is not None and "dpi" not in kwargs:
kwargs["dpi"] = dpi
return func(*args, **kwargs)
return wrapped
[docs]
def draw_features_from_shapefile(gp, name, **kwargs):
"""Add features as we need to."""
name2 = name if name != "borders" else "admin_0_boundary_lines_land"
boundpoly = gp.get_bounds_polygon()
# Life choices: which resolution to use?
threshold = 25 if gp.crs.is_geographic else 3e12
resolution = "50m" if boundpoly.area > threshold else "10m"
# Forward support cartopy_offlinedata variables (0.2 vs 0.20)
shpfn = os.path.join(
os.environ.get(
"CARTOPY_OFFLINE_SHARED", os.environ.get("CARTOPY_DATA_DIR")
),
"shapefiles",
"natural_earth",
"physical" if name != "borders" else "cultural",
f"ne_{resolution}_{name2}.shp",
)
fastfn = (
os.path.join(
"/opt/miniconda3/pyiem_data/parquet",
str(gp.crs.to_epsg()),
*shpfn.split(os.sep)[-3:],
).removesuffix(".shp")
+ ".parquet"
)
if os.path.isfile(fastfn):
df = gpd.read_parquet(fastfn).cx[
slice(*gp.get_xlim()), slice(*gp.get_ylim())
]
else:
df = (
gpd.read_file(shpfn, engine="pyogrio")
.to_crs(gp.crs)
.cx[slice(*gp.get_xlim()), slice(*gp.get_ylim())]
)
if not df.empty:
df.plot(ax=gp.ax, aspect=None, **kwargs)
[docs]
def ramp2df(name) -> pd.DataFrame:
"""Load pyIEM color ramp into a Pandas DataFrame.
Args:
name (str): the name of the bundled color ramp.
Returns
-------
pandas.DataFrame
"""
return pd.read_csv(f"{DATADIR}/ramps/{name}.txt")
[docs]
def pretty_bins(minval, maxval, bins=8):
"""Return a **smooth** binning that encloses the min and max value.
The returned array is at most the specified bins + 1 in size, but could
be smaller given this algorithm and the data range.
Args:
minval (real): minimum value to enclose.
maxval (real): maximum value to enclose.
bins (int): maximum number of bins to generate
Returns:
``np.array`` of bins
"""
center = (maxval + minval) / 2.0
return centered_bins(maxval - center, on=center, bins=bins)
[docs]
def centered_bins(absmax, on=0, bins=8):
"""Return a **smooth** binning around some number.
The returned array is +1 in size of the bins specified, since we want the
bin edges.
Args:
absmax (real): positive distance from the `on` value for bins to enclose.
on (real): where to center these bins.
bins (int): number of bins to generate
Returns:
``np.array`` of bins
"""
minval = on - absmax
maxval = on + absmax
vrange = maxval - minval
step = vrange / float(bins)
avail = [0.001, 0.005, 0.01, 0.05, 0.1, 0.25, 0.5, 1, 2, 5, 10, 25, 50]
avail.extend([100, 150, 200, 250, 500, 750, 1000, 1e4, 1e5, 1e6])
if step >= avail[-1]:
raise ValueError(f"step: {step} too large to handle")
# compute a new and round step value
step = avail[np.digitize(step, avail)]
# Create new min and max values
minval = np.floor(minval / step) * step
maxval = np.ceil(maxval / step) * step
res = np.arange(minval, maxval + 0.001, step)
# Account for some numerical instability
if len(res) == bins + 1 and abs(res[bins // 2] - on) < 0.001:
res[bins // 2] = on
return res
[docs]
def draw_logo(fig, logoname):
"""Place the logo."""
if logoname is None:
return
filename = LOGOFILES.get(logoname, "logo.png")
fn = os.path.join(DATADIR, filename)
# Create a fake axes to place this Logo
logo = mpimage.imread(fn)
# imshow messes with the aspect, so about the best we can do here is
# pin it to the upper edge
ax = fig.add_axes(
LOGO_BOUNDS, frameon=False, yticks=[], xticks=[], anchor="NW"
)
ax.imshow(logo, aspect="equal", zorder=-1)
[docs]
def fontscale(ratio, fig):
"""Return a font size suitable for this NDC ratio.
Args:
ratio (float): value between 0 and 1
fig (matplotlib.Figure,optional): The Figure of interest
Returns
-------
float: font size
"""
bbox = fig.get_window_extent().transformed(fig.dpi_scale_trans.inverted())
# point size is 72 pts per inch
return bbox.height * 72.0 * ratio
[docs]
def fitbox(fig, text, x0, x1, y0, y1, **kwargs):
"""Fit text into a NDC box.
Args:
textsize (int, optional): First attempt this textsize to see if it fits.
"""
if text is None:
return None
figbox = fig.get_window_extent().transformed(
fig.dpi_scale_trans.inverted()
)
# need some slop for decimal comparison below
px0 = x0 * fig.dpi * figbox.width - 0.15
px1 = x1 * fig.dpi * figbox.width + 0.15
py0 = y0 * fig.dpi * figbox.height - 0.15
py1 = y1 * fig.dpi * figbox.height + 0.15
xanchor = x0
if kwargs.get("ha", "") == "center":
xanchor = x0 + (x1 - x0) / 2.0
yanchor = y0
if kwargs.get("va", "") == "center":
yanchor = y0 + (y1 - y0) / 2.0
fontsize = int(kwargs.get("textsize", kwargs.get("fontsize", 50)))
txt = fig.text(
xanchor,
yanchor,
text,
fontsize=fontsize,
ha=kwargs.get("ha", "left"),
va=kwargs.get("va", "bottom"),
color=kwargs.get("color", "k"),
bbox=kwargs.get("bbox"),
)
def _fits(txt):
"""Test for fitting."""
tb = txt.get_window_extent(fig.canvas.get_renderer())
return tb.x0 >= px0 and tb.x1 < px1 and tb.y0 >= py0 and tb.y1 <= py1
if not _fits(txt):
for size in range(fontsize, 1, -2):
txt.set_fontsize(size)
if _fits(txt):
break
return txt
[docs]
def make_panel(
ndc_axbounds, fig, extent, crs, aspect, is_geoextent=False, **kwargs
) -> GeoPanel:
"""Make a GeoPanel.
Args:
ndc_axbounds (list): the NDC coordinates of axes to create
extent (list): x0,x1,y0,y1 *in projected space* plot extent, unless
`is_geoextent` is based as True, then it is Geodetic.
crs (pyproj.CRS): the crs of the axes
aspect (str): matplotlib's aspect of axes
is_geoextent(bool): is the passed extent Geodetic?
sector_label (bool): A Label that tracks what this is called
background (str): background to use.
Returns
-------
GeoPanel: the panel
"""
# https://jdhao.github.io/2017/06/03/change-aspect-ratio-in-mpl/
# aspect is data ratio divided by axes ratio
sector_label = kwargs.get("sector_label", "")
hacky = "spherical_mercator"
cls = SphericalMercatorPanel if sector_label == hacky else GeoPanel
gp = cls(
fig,
ndc_axbounds,
crs,
aspect=aspect,
adjustable="datalim",
facecolor=(0.4471, 0.6235, 0.8117),
xticks=[],
yticks=[],
sector_label=sector_label,
)
# Get the frame at the proper zorder
for spine in gp.ax.spines.values():
spine.set_zorder(reference.Z_FRAME)
# Turn off autoscale so that we can control the axis limits
gp.ax.autoscale(False)
# Set the extent
gp.set_extent(extent, crs=None if is_geoextent else crs)
gp.draw_background(kwargs.get("background"))
return gp
[docs]
def sector_setter(mp, axbounds, **kwargs):
"""Use kwargs to set the MapPlot sector."""
aspect = kwargs.get("aspect", "auto") # !important
if mp.sector == "cwa":
mp.cwa = kwargs.get("cwa", "DMX")
gp = make_panel(
axbounds,
mp.fig,
[
reference.wfo_bounds[mp.cwa][0],
reference.wfo_bounds[mp.cwa][2],
reference.wfo_bounds[mp.cwa][1],
reference.wfo_bounds[mp.cwa][3],
],
reference.EPSG[3857],
aspect,
is_geoextent=True,
sector_label=f"cwa_{mp.cwa}",
**kwargs,
)
mp.panels.append(gp)
elif mp.sector == "fema_region":
mp.fema_region = int(kwargs.get("fema_region", 7))
gp = make_panel(
axbounds,
mp.fig,
[
reference.fema_region_bounds[mp.fema_region][0],
reference.fema_region_bounds[mp.fema_region][2],
reference.fema_region_bounds[mp.fema_region][1],
reference.fema_region_bounds[mp.fema_region][3],
],
reference.EPSG[3857],
aspect,
is_geoextent=True,
sector_label=f"fema_region_{mp.fema_region}",
**kwargs,
)
mp.panels.append(gp)
elif mp.sector == "state":
mp.state = kwargs.get("state", "IA")
gp = make_panel(
axbounds,
mp.fig,
[
reference.state_bounds[mp.state][0],
reference.state_bounds[mp.state][2],
reference.state_bounds[mp.state][1],
reference.state_bounds[mp.state][3],
],
reference.EPSG[3857 if mp.state != "AK" else 3467],
aspect,
is_geoextent=True,
sector_label=f"state_{mp.state}",
**kwargs,
)
mp.panels.append(gp)
elif mp.sector == "iowawfo":
gp = make_panel(
axbounds,
mp.fig,
[-99.6, -89.0, 39.8, 45.5],
reference.EPSG[3857],
aspect,
is_geoextent=True,
sector_label=mp.sector,
**kwargs,
)
mp.panels.append(gp)
elif mp.sector == "custom":
# 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"]],
kwargs.get("projection", reference.LATLON),
aspect,
is_geoextent="projection" not in kwargs,
sector_label=mp.sector,
**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,
mp.fig,
[-4.5e6, 4.3e6, -3.9e6, 3.8e6],
reference.EPSG[9311],
"auto",
sector_label=mp.sector,
**kwargs,
)
mp.panels.append(gp)
elif mp.sector in ["conus", "nws"]:
gp = make_panel(
axbounds,
mp.fig,
[-2400000, 2300000, 27600, 3173000],
reference.EPSG[5070],
aspect,
sector_label=mp.sector,
**kwargs,
)
mp.panels.append(gp)
if mp.sector == "nws":
width, height = mp.fig.canvas.get_width_height()
# Create PR
ln = mp.panels[0].plot(
[-80.5, -74.55],
[21.5, 22.8],
color="None",
)[0]
bbox = ln.get_window_extent(mp.fig.canvas.get_renderer())
mp.panels.append(
make_panel(
[
bbox.x0 / width,
bbox.y0 / height,
(bbox.x1 - bbox.x0) / width,
(bbox.y1 - bbox.y0) / height,
],
mp.fig,
[-68.0, -65.0, 17.5, 18.6],
reference.LATLON,
aspect,
is_geoextent=True,
sector_label="state_PR",
**kwargs,
)
)
fitbox(
mp.fig,
"Puerto Rico",
bbox.x0 / width,
bbox.x1 / width,
bbox.y1 / height,
bbox.y1 / height + 0.1,
color="white",
bbox=dict(boxstyle="round,pad=0", color="k"),
)
# Create AK
ln = mp.panels[0].plot(
[-118.0, -105.5],
[20.0, 30.25],
color="None",
)[0]
bbox = ln.get_window_extent(mp.fig.canvas.get_renderer())
mp.panels.append(
make_panel(
[
bbox.x0 / width,
bbox.y0 / height,
(bbox.x1 - bbox.x0) / width,
(bbox.y1 - bbox.y0) / height,
],
mp.fig,
[
reference.state_bounds["AK"][0],
reference.state_bounds["AK"][2],
reference.state_bounds["AK"][1],
reference.state_bounds["AK"][3],
],
reference.EPSG[3467],
aspect,
is_geoextent=True,
sector_label="state_AK",
**kwargs,
)
)
# Guam
ln = mp.panels[0].plot(
[-120.5, -117.5],
[28.3, 31.7],
color="None",
)[0]
bbox = ln.get_window_extent(mp.fig.canvas.get_renderer())
mp.panels.append(
make_panel(
[
bbox.x0 / width,
bbox.y0 / height,
(bbox.x1 - bbox.x0) / width,
(bbox.y1 - bbox.y0) / height,
],
mp.fig,
[
reference.state_bounds["GU"][0],
reference.state_bounds["GU"][2],
reference.state_bounds["GU"][1],
reference.state_bounds["GU"][3],
],
reference.LATLON,
aspect,
is_geoextent=True,
sector_label="state_GU",
**kwargs,
)
)
fitbox(
mp.fig,
"Guam",
bbox.x0 / width,
bbox.x1 / width,
bbox.y1 / height,
bbox.y1 / height + 0.1,
color="white",
bbox=dict(boxstyle="round,pad=0", color="k"),
)
# Hawaii
ln = mp.panels[0].plot(
[-95.8, -85.24],
[22.9, 27.7],
color="None",
)[0]
bbox = ln.get_window_extent(mp.fig.canvas.get_renderer())
mp.panels.append(
make_panel(
[
bbox.x0 / width,
bbox.y0 / height,
(bbox.x1 - bbox.x0) / width,
(bbox.y1 - bbox.y0) / height,
],
mp.fig,
[-161.0, -154.0, 18.5, 22.5],
reference.LATLON,
aspect,
is_geoextent=True,
sector_label="state_HI",
**kwargs,
)
)
# Do last in case of name overlaps above
elif mp.sector in reference.SECTORS:
gp = make_panel(
axbounds,
mp.fig,
reference.SECTORS[mp.sector],
reference.EPSG[3857],
aspect,
is_geoextent=True,
sector_label=mp.sector,
**kwargs,
)
mp.panels.append(gp)
# Kindof a hack for now
if mp.panels:
mp.ax = mp.panels[0].ax
def _mask_with_path(gp, path):
"""Make a path."""
# Removes any external data
patch = mpatches.PathPatch(
path,
facecolor="white",
edgecolor="none",
zorder=reference.Z_CLIP,
)
gp.ax.add_patch(patch)
# Then gives a nice semitransparent look
patch = mpatches.PathPatch(
path,
facecolor="black",
edgecolor="none",
zorder=reference.Z_CLIP2,
alpha=0.65,
)
gp.ax.add_patch(patch)
[docs]
def mask_outside_polygon(poly_verts, gp):
"""Make outside of a polygon.
POLY_VERTS is in CCW order, as this is the interior of the polygon
"""
# Get current plot limits
xlim = gp.get_xlim()
ylim = gp.get_ylim()
# Verticies of the plot boundaries in clockwise order
bound_verts = np.array(
[
(xlim[0], ylim[0]),
(xlim[0], ylim[1]),
(xlim[1], ylim[1]),
(xlim[1], ylim[0]),
(xlim[0], ylim[0]),
]
)
# A series of codes (1 and 2) to tell matplotlib whether to draw a lineor
# move the "pen" (So that there's no connecting line)
bound_codes = [mpath.Path.MOVETO] + (len(bound_verts) - 1) * [
mpath.Path.LINETO
]
poly_codes = [mpath.Path.MOVETO] + (len(poly_verts) - 1) * [
mpath.Path.LINETO
]
# Plot the masking patch
path = mpath.Path(
np.concatenate([bound_verts, poly_verts]), bound_codes + poly_codes
)
_mask_with_path(gp, path)
# Reset the plot limits to their original extents
gp.ax.set_xlim(xlim)
gp.ax.set_ylim(ylim)
[docs]
def polygon_fill(mymap, geodf, data, **kwargs):
"""Generalize function for overlaying filled polygons on the map.
Args:
mymap (MapPlot): The MapPlot instance
geodf (GeoDataFrame): A GeoDataFrame with a `geom` column.
data (dict): The dictionary of keys and values used for picking colors
Keyword Args:
ilabel (Optional[bool]): should values be labelled? Defaults to `False`
lblfmt (str,optional): format string for labels. Defaults to %s.
plotmissing (bool): should geometries not included in the `data`
be mapped? Defaults to `True`
color (str or dict): Providing an explicit color (used for both edge
and fill). Either provide one color or a dictionary to lookup a
color by the mapping key.
fc (str or dict): Same as `color`, but controls the fill color.
Providing this value will over-ride any `color` setting.
ec (str or dict): Same as `color`, but controls the edge color.
Providing this value will over-ride any `color` setting.
zorder (int): The zorder to use for this layer, default `Z_FILL`
lw (float): polygon outline width
"""
bins = kwargs.get("bins", np.arange(0, 101, 10))
cmap = stretch_cmap(kwargs.get("cmap"), bins, extend=kwargs.get("extend"))
ilabel = kwargs.get("ilabel", False)
norm = mpcolors.BoundaryNorm(bins, cmap.N)
lblformat = kwargs.get("lblformat", "%s")
plotmissing = kwargs.get("plotmissing", True)
color = kwargs.get("color")
ec = kwargs.get("ec")
fc = kwargs.get("fc")
labels = kwargs.get("labels", {})
zorder = kwargs.get("zorder", reference.Z_FILL)
to_label = {"x": [], "y": [], "vals": []}
# Merge data into the data frame
geodf["val"] = pd.Series(data)
# Sort the frame, placing the missing values at the beginning
geodf = geodf.sort_values("val", na_position="first")
if not plotmissing:
geodf = geodf[pd.notna(geodf["val"])].copy()
for gp in mymap.panels:
# Reproject data into plot native, found out that clip() sometimes
# returns a GeometryCollection, which geopandas hates to plot
native = geodf.to_crs(gp.crs).cx[
slice(*gp.get_xlim()), slice(*gp.get_ylim())
]
if native.empty:
continue
native = native.copy()
native["fc"] = pd.Series(fc) if fc else "white"
native["ec"] = pd.Series(ec) if ec else "k"
# Compute colors and such
for idx, row in native.iterrows():
lbl = (
labels.get(idx, "-")
if pd.isna(row["val"])
else labels.get(idx, lblformat % (row["val"],))
)
# How we compute a fill and edge color
_fc, _ec = (None, None)
if color is not None:
if isinstance(color, str):
_fc, _ec = (color, color)
else:
_fc = color.get(idx, "white")
_ec = color.get(idx, "k")
if fc is not None:
_fc = fc if isinstance(fc, str) else fc.get(idx, "white")
if ec is not None:
_ec = ec if isinstance(ec, str) else ec.get(idx, "k")
_ec = "k" if _ec is None else _ec
if _fc is None:
_fc = (
"white"
if pd.isna(row["val"])
else mpcolors.to_hex(cmap(norm([row["val"]]))[0])
)
native.at[idx, "fc"] = _fc
native.at[idx, "ec"] = _ec
if ilabel:
# prefer our stored centroid vs calculated one
mx = row.get("lon", native.at[idx, "geom"].centroid.x)
my = row.get("lat", native.at[idx, "geom"].centroid.y)
to_label["x"].append(mx)
to_label["y"].append(my)
to_label["vals"].append(lbl)
native.plot(
ax=gp.ax,
fc=native["fc"].values,
ec=native["ec"].values,
aspect=None,
zorder=zorder,
lw=kwargs.get("lw", 0.1),
)
if to_label:
mymap.plot_values(
to_label["x"],
to_label["y"],
to_label["vals"],
labelbuffer=kwargs.get("labelbuffer", 1),
textsize=12,
textoutlinewidth=2,
clip_on=True,
)
kwargs.pop("cmap", None)
kwargs.pop("bins", None)
if kwargs.pop("draw_colorbar", True):
mymap.draw_colorbar(bins, cmap, norm, **kwargs)
[docs]
def mask_outside_geom(gp, geom):
"""Create a white patch over the plot for what we want to ask out.
Args:
gp (GeoPanel): The GeoPanel instance
geom (geometry):
"""
xlim = gp.get_xlim()
ylim = gp.get_ylim()
# Verticies of the plot boundaries in clockwise order
verts = np.array(
[
(xlim[0], ylim[0]),
(xlim[0], ylim[1]),
(xlim[1], ylim[1]),
(xlim[1], ylim[0]),
(xlim[0], ylim[0]),
]
)
codes = [mpath.Path.MOVETO] + (len(verts) - 1) * [mpath.Path.LINETO]
if isinstance(geom, Polygon):
geom = MultiPolygon([geom])
trans = Transformer.from_crs(LATLON, gp.crs, always_xy=True)
for geo in geom.geoms:
ccw = np.asarray(geo.exterior.coords)
if not geo.exterior.is_ccw:
ccw = ccw[::-1]
points = trans.transform(ccw[:, 0], ccw[:, 1])
# pylint: disable=unsubscriptable-object
ar = np.column_stack([*points])[:-1]
verts = np.concatenate([verts, ar])
codes = np.concatenate(
[
codes,
[mpath.Path.MOVETO] + (ar.shape[0] - 1) * [mpath.Path.LINETO],
]
)
path = mpath.Path(verts, codes)
_mask_with_path(gp, path)
# Reset the plot limits to their original extents
gp.ax.set_xlim(xlim)
gp.ax.set_ylim(ylim)