Merge pull request #697 from akrherz/gh304_interstates_lakes

introduce map background for MapPlot

CHANGELOG.md

@@ -28,6 +28,7 @@ database has.
 - Jabber channels for METAR wind gust alerts were enhanced (#683).
 - 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).
 - 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/data/backgrounds/README.md

@@ -0,0 +1,15 @@
+# pyIEM MapPlot backgrounds
+
+The goal is to have pretty things in the background of plots and to remember
+how I set this mess up!  So the directory nomenclature is as such:
+
+    <name>/<sector>_<epsg>.{png,wld}
+
+There should be a `default_4326.png` that covers the NW quad of the globe
+in a modest fashion.  We then go down the rabbit hole of providing pre-projected
+scenes for various common sector + epsg combinations.
+
+The ne2 source is [here](https://www.naturalearthdata.com/downloads/10m-raster-data/10m-natural-earth-2/),
+"Natural Earth II with Shaded Relief, Water, and Drainages".
+
+There's a `util/cookie_cutter.py` that generates scenes we want.

src/pyiem/data/backgrounds/ne2/conus_5070.wld

@@ -0,0 +1,6 @@
+5076.8400653595
+0.0000000000
+0.0000000000
+-3378.8476562500
+-3296639.2217973866
+3328580.5761718750

src/pyiem/data/backgrounds/ne2/default_4326.wld

@@ -0,0 +1,6 @@
+0.1546875000
+0.0000000000
+0.0000000000
+-0.1029509479
+-188.9226562500
+97.6594098267

src/pyiem/plot/_mpl.py

@@ -1,15 +1,79 @@
 """A custom axes implementation."""
 # pylint: disable=unsubscriptable-object,unpacking-non-sequence
+import os
 
 # Third party libraries
+import numpy as np
 from shapely.geometry import Polygon
 import geopandas as gpd
 from pyproj import Transformer
+import rasterio
+from rasterio.warp import reproject, Resampling
 
 # Local imports
 from pyiem.reference import LATLON
 
 
+def draw_background(panel, background):
+    """Draw the background for this plot!"""
+    if background is None:
+        return
+    datadirs = [
+        os.sep.join(
+            [
+                os.path.dirname(__file__),
+                "..",
+                "data",
+                "backgrounds",
+                background,
+            ]
+        ),
+        f"/opt/miniconda3/pyiem_data/backgrounds/{background}",
+    ]
+    src_epsg = str(panel.crs)
+    rasterfn = f"{panel.get_sector_label()}_{src_epsg.split(':')[1]}.png"
+    full = os.sep.join([datadirs[0], rasterfn])
+    if not os.path.isfile(full):
+        full = os.sep.join([datadirs[1], rasterfn])
+        if not os.path.isfile(full):
+            rasterfn = "default_4326.png"
+            full = os.sep.join([datadirs[0], rasterfn])
+            src_epsg = "EPSG:4326"
+    worldfn = f"{full[:-4]}.wld"
+    with open(worldfn, encoding="ascii") as fh:
+        (dx, _, _, dy, west, north) = [float(x) for x in fh]
+    src_aff = rasterio.Affine(dx, 0, west, 0, dy, north)
+    src_crs = {"init": src_epsg}
+    (px0, px1, py0, py1) = panel.get_extent()
+    pbbox = panel.ax.get_window_extent()
+    pdx = (px1 - px0) / pbbox.width
+    pdy = (py1 - py0) / pbbox.height
+    dest_aff = rasterio.Affine(pdx, 0, px0, 0, pdy, py0)
+    res = np.zeros((int(pbbox.height), int(pbbox.width), 3), dtype=np.uint8)
+    band = np.zeros((int(pbbox.height), int(pbbox.width)), dtype=np.uint8)
+    with rasterio.open(full) as src:
+        data = src.read()
+        for i in range(3):
+            reproject(
+                data[i, :, :],
+                band,
+                src_transform=src_aff,
+                src_crs=src_crs,
+                src_nodata=0,
+                dst_transform=dest_aff,
+                dst_crs={"init": str(panel.crs)},
+                resampling=Resampling.nearest,
+            )
+            res[:, :, i] = band
+    panel.ax.imshow(
+        res / 255.0,
+        interpolation="nearest",  # prevents artifacts
+        extent=(px0, px1, py0, py1),
+        origin="lower",
+        zorder=-1,
+    ).set_rasterized(True)
+
+
 class GeoPanel:
     """
     A container class holding references to a matplotlib axes.
@@ -19,10 +83,17 @@ def __init__(self, fig, rect, crs, **kwargs):
         """
         Initialize the axes.
         """
+        self.sector_label = kwargs.pop("sector_label", "")
         self.figure = fig
         self.ax = self.figure.add_axes(rect, **kwargs)
         self.crs = crs
 
+    draw_background = draw_background
+
+    def get_sector_label(self) -> str:
+        """Return the property"""
+        return self.sector_label
+
     def get_bounds_polygon(self):
         """Return the axes extent as a shapely polygon bounds."""
         xlim = self.ax.get_xlim()

src/pyiem/plot/geoplot.py

@@ -146,6 +146,8 @@ def __init__(self, sector="iowa", **kwargs):
               [left, bottom, width, height] for the main axes.
             stateborderwidth (float,optional): how wide to make the
               state borders (default: 1.).
+            background (str,optional): Background imagery to use `ne2` is the
+              only option currently.
         """
         self.debug = kwargs.get("debug", False)
         self.fig = kwargs.get("fig")
@@ -170,12 +172,16 @@ def __init__(self, sector="iowa", **kwargs):
             self.state = "IA"
         axes_position = kwargs.pop("axes_position", MAIN_AX_BOUNDS)
         sector_setter(self, axes_position, **kwargs)
+        has_background = kwargs.get("background") is not None
         for gp in self.panels:
             # legacy usage of axisbg here
             _c = kwargs.get(
                 "axisbg", kwargs.get("continentalcolor", "#EEEEEE")
             )
-            draw_features_from_shapefile(gp, "land", facecolor=_c, zorder=Z_CF)
+            if not has_background:
+                draw_features_from_shapefile(
+                    gp, "land", facecolor=_c, zorder=Z_CF
+                )
             # NB we neeed both borders (lines between countries) and
             # coastlines (lines between land and water)
             draw_features_from_shapefile(
@@ -192,13 +198,14 @@ def __init__(self, sector="iowa", **kwargs):
                 ec="k",
                 zorder=Z_POLITICAL,
             )
-            draw_features_from_shapefile(
-                gp,
-                "lakes",
-                edgecolor=(0.4471, 0.6235, 0.8117),
-                facecolor=(0.4471, 0.6235, 0.8117),
-                zorder=Z_CF,
-            )
+            if not has_background:
+                draw_features_from_shapefile(
+                    gp,
+                    "lakes",
+                    edgecolor=(0.4471, 0.6235, 0.8117),
+                    facecolor=(0.4471, 0.6235, 0.8117),
+                    zorder=Z_CF,
+                )
             if "nostates" not in kwargs:
                 xlim = gp.ax.get_xlim()
                 ylim = gp.ax.get_ylim()

src/pyiem/plot/util.py

@@ -218,7 +218,7 @@ def _fits(txt):
 
 
 def make_panel(
-    ndc_axbounds, fig, extent, crs, aspect, is_geoextent=False
+    ndc_axbounds, fig, extent, crs, aspect, is_geoextent=False, **kwargs
 ) -> GeoPanel:
     """Factory for making a GeoPanel.
 
@@ -229,6 +229,8 @@ def make_panel(
       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
@@ -244,6 +246,7 @@ def make_panel(
         facecolor=(0.4471, 0.6235, 0.8117),
         xticks=[],
         yticks=[],
+        sector_label=kwargs.get("sector_label", ""),
     )
     # Get the frame at the proper zorder
     for _k, spine in gp.ax.spines.items():
@@ -252,6 +255,7 @@ def make_panel(
     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
 
 
@@ -273,6 +277,8 @@ def sector_setter(mp, axbounds, **kwargs):
             reference.EPSG[3857],
             aspect,
             is_geoextent=True,
+            sector_label=f"cwa_{mp.cwa}",
+            **kwargs,
         )
         mp.panels.append(gp)
     elif mp.sector == "state":
@@ -289,6 +295,8 @@ def sector_setter(mp, axbounds, **kwargs):
             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":
@@ -299,6 +307,8 @@ def sector_setter(mp, axbounds, **kwargs):
             reference.EPSG[3857],
             aspect,
             is_geoextent=True,
+            sector_label=mp.sector,
+            **kwargs,
         )
         mp.panels.append(gp)
     elif mp.sector == "custom":
@@ -313,6 +323,8 @@ def sector_setter(mp, axbounds, **kwargs):
             kwargs.get("projection", reference.LATLON),
             aspect,
             is_geoextent="projection" not in kwargs,
+            sector_label=mp.sector,
+            **kwargs,
         )
         mp.panels.append(gp)
 
@@ -323,6 +335,8 @@ def sector_setter(mp, axbounds, **kwargs):
             [-4.5e6, 4.3e6, -3.9e6, 3.8e6],
             reference.EPSG[2163],
             "auto",
+            sector_label=mp.sector,
+            **kwargs,
         )
         mp.panels.append(gp)
 
@@ -333,6 +347,8 @@ def sector_setter(mp, axbounds, **kwargs):
             [-2400000, 2300000, 27600, 3173000],
             reference.EPSG[5070],
             aspect,
+            sector_label=mp.sector,
+            **kwargs,
         )
         mp.panels.append(gp)
 
@@ -345,6 +361,8 @@ def sector_setter(mp, axbounds, **kwargs):
                 reference.LATLON,
                 aspect,
                 is_geoextent=True,
+                sector_label="state_PR",
+                **kwargs,
             )
             mp.fig.text(
                 0.78,
@@ -368,6 +386,8 @@ def sector_setter(mp, axbounds, **kwargs):
                 reference.EPSG[3467],
                 aspect,
                 is_geoextent=True,
+                sector_label="state_AK",
+                **kwargs,
             )
             mp.panels.append(gp)
             # Guam
@@ -383,6 +403,8 @@ def sector_setter(mp, axbounds, **kwargs):
                 reference.LATLON,
                 aspect,
                 is_geoextent=True,
+                sector_label="state_GU",
+                **kwargs,
             )
             mp.fig.text(
                 0.015,
@@ -409,6 +431,8 @@ def sector_setter(mp, axbounds, **kwargs):
                 reference.LATLON,
                 aspect,
                 is_geoextent=True,
+                sector_label="state_HI",
+                **kwargs,
             )
             mp.panels.append(gp)
     # Do last in case of name overlaps above
@@ -420,6 +444,8 @@ def sector_setter(mp, axbounds, **kwargs):
             reference.EPSG[3857],
             aspect,
             is_geoextent=True,
+            sector_label=mp.sector,
+            **kwargs,
         )
         mp.panels.append(gp)
     # Kindof a hack for now
@@ -529,7 +555,7 @@ def polygon_fill(mymap, geodf, data, **kwargs):
     # Merge data into the data frame
     geodf["val"] = pd.Series(data)
     if not plotmissing:
-        geodf = geodf[~pd.isna(geodf["val"])].copy()
+        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

tests/plot/test_geoplot.py

@@ -1,5 +1,5 @@
 """Test plots made by pyiem.plot.geoplot"""
-# pylint disable=too-many-lines
+# pylint: disable=too-many-lines
 import datetime
 import tempfile
 import os
@@ -47,6 +47,31 @@ def test_invalid_file():
     assert load_bounds("this shall not work") is None
 
 
+@pytest.mark.mpl_image_compare(tolerance=PAIN)
+def test_conus_background():
+    """Test that a conus sector plot with a background!"""
+    mp = MapPlot(
+        nocaption=True, sector="conus", background="ne2", twitter=True
+    )
+    return mp.fig
+
+
+@pytest.mark.mpl_image_compare(tolerance=PAIN)
+def test_custom_background():
+    """Test we get a background when a custom sector is picked"""
+    mp = MapPlot(
+        sector="custom",
+        background="ne2",
+        twitter=True,
+        north=34,
+        east=-84,
+        south=28,
+        west=-90,
+        nocaption=True,
+    )
+    return mp.fig
+
+
 @pytest.mark.mpl_image_compare(tolerance=PAIN)
 def test_force_cities():
     """Test that we can force cities underneath some data."""