Merge pull request #956 from akrherz/fema_regions

✨ Support plotting by FEMA Regions

CHANGELOG.md

@@ -3,6 +3,16 @@
 
 All notable changes to this library are documented in this file.
 
+## Unreleased Version
+
+### API Changes
+
+### New Features
+
+- Add support for plotting by FEMA Regions.
+
+### Bug Fixes
+
 ## **1.21.0** (6 Sep 2024)
 
 ### API Changes

src/pyiem/plot/geoplot.py

@@ -165,6 +165,7 @@ def __init__(self, sector="iowa", **kwargs):
         # Storage of axes within this plot
         self.state = None
         self.cwa = None
+        self.fema_region = None
         self.textmask = None  # For our plot_values magic, to prevent overlap
         self.sector = sector
         self.cax = plt.axes(CAX_BOUNDS, frameon=False, yticks=[], xticks=[])
@@ -992,6 +993,7 @@ def draw_mask(self, sector=None):
             "state",
             "iowawfo",
             "cwa",
+            "fema_region",
         ):
             return
         # in lon,lat
@@ -1003,6 +1005,10 @@ def draw_mask(self, sector=None):
             s = load_geodf("cwa")
             mask_outside_geom(self.panels[0], s.at[self.cwa, "geom"])
             return
+        if sector == "fema_region":
+            s = load_geodf("fema_regions")
+            mask_outside_geom(self.panels[0], s.at[self.fema_region, "geom"])
+            return
         if sector == "conus":
             s = load_geodf("conus")
             mask_outside_geom(self.panels[0], s.iloc[0]["geom"])
@@ -1155,6 +1161,27 @@ def fill_states(self, data, **kwargs):
             geodf.loc[st, "lat"] = geodf.loc[st, "lat"] + 0.3
         polygon_fill(self, geodf, data, **kwargs)
 
+    def draw_fema_regions(self, color: str = "k", **kwargs):
+        """Overlay FEMA Regions."""
+        kwargs["edgecolor"] = color
+        regions = load_geodf("fema_regions")
+        for gp in self.panels:
+            regions.to_crs(gp.crs).plot(
+                ax=gp.ax,
+                aspect=None,
+                zorder=Z_POLITICAL,
+                facecolor="None",
+                **kwargs,
+            )
+
+    def fill_fema_regions(self, data, **kwargs):
+        """Add overlay of filled polygons for FEMA Regions.
+
+        Data is dictionary-ish and keys should be ints!
+        """
+        geodf = load_geodf("fema_regions")
+        polygon_fill(self, geodf, data, **kwargs)
+
     def draw_cwas(self, color="k", **kwargs):
         """Overlay CWA Borders
 

src/pyiem/plot/util.py

@@ -309,6 +309,24 @@ def sector_setter(mp, axbounds, **kwargs):
             **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(

src/pyiem/reference.py

@@ -667,6 +667,18 @@
     "VEF": [-118.99, 33.83, -112.33, 38.88],
 }
 
+fema_region_bounds = {
+    1: [-73.93, 40.75, -66.69, 47.66],
+    2: [-79.96, 17.48, -64.37, 45.21],
+    3: [-83.88, 36.34, -74.50, 42.47],
+    4: [-91.84, 24.76, -75.26, 39.34],
+    5: [-97.43, 36.79, -80.32, 49.57],
+    6: [-109.25, 25.65, -88.82, 37.20],
+    7: [-104.26, 35.79, -88.91, 43.70],
+    8: [-116.26, 36.79, -96.24, 49.20],
+    9: [-178.64, -14.58, 146.28, 42.20],
+    10: [-179.43, 41.79, -110.85, 71.63],
+}
 
 # State FIPS
 state_fips = {

tests/plot/test_geoplot.py

@@ -222,6 +222,23 @@ def test_nws_sector_twitter_res():
     return mp.fig
 
 
+@pytest.mark.mpl_image_compare(tolerance=PAIN)
+def test_fema_region6():
+    """Test drawing fema regions."""
+    mp = MapPlot(
+        figsize=(8, 6),
+        nocaption=True,
+        fema_region=6,
+        sector="fema_region",
+        title="FEMA Region 6",
+        nostates=True,
+    )
+    mp.fill_fema_regions({6: 10})
+    mp.draw_fema_regions()
+    mp.draw_mask()
+    return mp.fig
+
+
 @pytest.mark.mpl_image_compare(tolerance=PAIN)
 def test_nashville():
     """Test that Benton County, TNC005 does not show for OHX."""

util/make_shapes.py

@@ -5,6 +5,7 @@
 import warnings
 
 import geopandas as gpd
+from sqlalchemy import text
 
 from pyiem.database import get_sqlalchemy_conn
 from pyiem.reference import state_bounds
@@ -102,6 +103,20 @@ def dump_iowawfo(fn):
     df.to_parquet(fn)
 
 
+def dump_fema_regions(fn):
+    """Dump fema regions."""
+    with get_sqlalchemy_conn("postgis", user="nobody") as conn:
+        df = gpd.read_postgis(
+            text(""" SELECT region, states,
+                ST_MakeValid(ST_Simplify(geom, 0.01)) as geom
+            from fema_regions"""),
+            conn,
+            geom_col="geom",
+            index_col="region",
+        )
+    df.to_parquet(fn)
+
+
 def dump_ugc(gtype, fn, is_firewx=False):
     """UGCS."""
     source_limiter = "source != 'fz'"
@@ -146,6 +161,8 @@ def getfn(prefix):
 
 def main():
     """Go Main"""
+    dump_fema_regions(getfn("fema_regions"))
+    check_file(getfn("fema_regions"))
     dump_conus(getfn("conus"))
     check_file(getfn("conus"))
     dump_iowawfo(getfn("iowawfo"))