county labels for map

2025-11-01 15:32:39 -07:00
parent 33a3a96126
commit f476eb26a3
4 changed files with 101 additions and 7 deletions
--- a/scripts/make_spill_map.py
+++ b/scripts/make_spill_map.py
@@ -153,12 +153,24 @@ def get_colorado_boundary(spills_gdf: gpd.GeoDataFrame) -> tuple[gpd.GeoDataFram
                        cg = gdf[gdf[col].astype(str).str.upper() == str(val).upper()]
                    if not cg.empty:
                        break
-            if cg.empty:
-                cg = gdf  # If filtering failed, use all (last resort)
+            # If attribute filtering failed or is clearly not limited to CO, use spatial filter based on spills envelope
+            if cg.empty or len(cg) > 200:  # heuristic: many features implies nationwide
+                # Build a coarse CO mask from spill points
+                mask_geom = spills_gdf.to_crs(3857).unary_union.convex_hull.buffer(40_000)
+                mask_gdf = gpd.GeoDataFrame(geometry=[mask_geom], crs='EPSG:3857').to_crs(4326)
+                try:
+                    cg = gdf[gdf.to_crs(4326).intersects(mask_gdf.geometry.iloc[0])]
+                    if cg.empty:
+                        cg = gdf  # last resort
+                except Exception:
+                    cg = gdf
+
+            # Use the (possibly filtered) counties for overlay and boundary
            counties = cg.to_crs(4326)
            # Dissolve to single state boundary
            boundary = counties.dissolve().reset_index(drop=True)
-            return boundary[["geometry"]], counties[["geometry"]]
+            # Return boundary geometry only, but keep ALL county attributes for labeling
+            return boundary[["geometry"]], counties
        except Exception:
            pass

@@ -197,7 +209,7 @@ def get_cities_gdf():
        {"name": "Fort Collins", "lat": 40.5853, "lon": -105.0844},
        {"name": "Pueblo", "lat": 38.2544, "lon": -104.6091},
        {"name": "Grand Junction", "lat": 39.0639, "lon": -108.5506},
-        {"name": "Greeley", "lat": 40.4233, "lon": -104.7091},
+    # Removed Greeley label per request
        {"name": "Boulder", "lat": 40.01499, "lon": -105.2705},
        {"name": "Durango", "lat": 37.2753, "lon": -107.8801},
    ]
@@ -219,6 +231,73 @@ def plot_cities(ax, cities_gdf_3857):
            path_effects=[pe.withStroke(linewidth=2.5, foreground='white')]
        )

+def label_selected_counties(ax, counties_gdf: gpd.GeoDataFrame | None):
+    """Label Weld, Garfield, and Rio Blanco counties if available.
+    Expects counties in EPSG:4326; will project to 3857 for plotting.
+    """
+    if counties_gdf is None or counties_gdf.empty:
+        return
+    gdf = counties_gdf.copy()
+    # Try common county name fields
+    name_col = None
+    for c in ['NAME', 'Name', 'COUNTY', 'County', 'county', 'NAMELSAD', 'NAMELSAD10', 'NAME10', 'name']:
+        if c in gdf.columns:
+            name_col = c
+            break
+    if name_col is None:
+        return
+    targets = {n.upper() for n in ['Weld', 'Garfield', 'Rio Blanco']}
+    gdf['__nm__'] = gdf[name_col].astype(str).str.strip()
+    gdf['__NMU__'] = gdf['__nm__'].str.upper()
+    sel = gdf[gdf['__NMU__'].isin(targets)]
+    if sel.empty:
+        # Fallback: contains matching (handles NAMELSAD like "Weld County")
+        mask = False
+        for t in targets:
+            mask = mask | gdf['__NMU__'].str.contains(t)
+        sel = gdf[mask]
+    if sel.empty:
+        return
+    sel = sel.to_crs(3857)
+    # Centroids for labels
+    for _, row in sel.iterrows():
+        centroid = row.geometry.representative_point()
+        ax.annotate(
+            row['__nm__'],
+            (centroid.x, centroid.y),
+            xytext=(0, 0), textcoords='offset points',
+            fontsize=8, color='black', ha='center', va='center', zorder=7,
+            path_effects=[pe.withStroke(linewidth=2.5, foreground='white')]
+        )
+
+def bold_target_county_borders(ax, counties_gdf: gpd.GeoDataFrame | None):
+    """Emphasize borders for Weld, Garfield, and Rio Blanco counties by drawing thicker dark outlines."""
+    if counties_gdf is None or counties_gdf.empty:
+        return
+    gdf = counties_gdf.copy()
+    # Try common county name fields
+    name_col = None
+    for c in ['NAME', 'Name', 'COUNTY', 'County', 'county', 'NAMELSAD', 'NAMELSAD10', 'NAME10', 'name']:
+        if c in gdf.columns:
+            name_col = c
+            break
+    if name_col is None:
+        return
+    targets = {n.upper() for n in ['Weld', 'Garfield', 'Rio Blanco']}
+    gdf['__nm__'] = gdf[name_col].astype(str).str.strip()
+    gdf['__NMU__'] = gdf['__nm__'].str.upper()
+    sel = gdf[gdf['__NMU__'].isin(targets)]
+    if sel.empty:
+        # Fallback contains match for NAMELSAD variants
+        mask = False
+        for t in targets:
+            mask = mask | gdf['__NMU__'].str.contains(t)
+        sel = gdf[mask]
+    if sel.empty:
+        return
+    sel = sel.to_crs(3857)
+    sel.boundary.plot(ax=ax, color='#222222', linewidth=1.4, zorder=5)
+

 def map_points(spills: gpd.GeoDataFrame, boundary: gpd.GeoDataFrame, counties: gpd.GeoDataFrame | None):
    # Project to Web Mercator for plotting
@@ -227,10 +306,10 @@ def map_points(spills: gpd.GeoDataFrame, boundary: gpd.GeoDataFrame, counties: g

    # Figure aesthetics
    fig, ax = plt.subplots(figsize=(7, 9), dpi=300, constrained_layout=True)
-    boundary_3857.plot(ax=ax, color='#f5f5f5', edgecolor='#888888', linewidth=0.8)
+    boundary_3857.plot(ax=ax, color='#f5f5f5', edgecolor='#888888', linewidth=0.8, zorder=1)
    if counties is not None and not counties.empty:
        counties_3857 = counties.to_crs(3857)
-        counties_3857.boundary.plot(ax=ax, color='#aaaaaa', linewidth=0.3)
+        counties_3857.boundary.plot(ax=ax, color='#bbbbbb', linewidth=0.4, zorder=2)

    # Plot points with alpha and small size to avoid overplotting
    spills_3857.plot(ax=ax, markersize=1.2, alpha=0.35, color='#2c7fb8')
@@ -238,6 +317,10 @@ def map_points(spills: gpd.GeoDataFrame, boundary: gpd.GeoDataFrame, counties: g
    # Cities
    cities = get_cities_gdf().to_crs(3857)
    plot_cities(ax, cities)
+    # County labels
+    label_selected_counties(ax, counties)
+    # Bold borders for target counties
+    bold_target_county_borders(ax, counties)

    # Tighten to boundary
    minx, miny, maxx, maxy = boundary_3857.total_bounds
@@ -282,6 +365,9 @@ def map_points(spills: gpd.GeoDataFrame, boundary: gpd.GeoDataFrame, counties: g
    # Cities on inset
    cities_inset = get_cities_gdf().to_crs(3857)
    plot_cities(inset_ax, cities_inset)
+    # County labels on inset
+    label_selected_counties(inset_ax, counties)
+    bold_target_county_borders(inset_ax, counties)
    inset_ax.set_xlim(bxmin, bxmax)
    inset_ax.set_ylim(bymin, bymax)
    style_axes(inset_ax)
@@ -314,13 +400,21 @@ def map_hex(spills: gpd.GeoDataFrame, boundary: gpd.GeoDataFrame, counties: gpd.
        counties_3857.boundary.plot(ax=ax, color='#aaaaaa', linewidth=0.3)

    # Use hexbin; gridsize tuned to CO scale; mincnt to suppress isolated cells
-    hb = ax.hexbin(x, y, gridsize=70, mincnt=3, cmap='viridis', alpha=0.9)
+    hb = ax.hexbin(x, y, gridsize=70, mincnt=3, cmap='viridis', alpha=0.9, zorder=3)
    cb = fig.colorbar(hb, ax=ax, fraction=0.028, pad=0.01)
    cb.set_label('Spill count')

+    # Draw county boundaries on top for visibility
+    if counties is not None and not counties.empty:
+        counties_3857.boundary.plot(ax=ax, color='#666666', linewidth=0.6, zorder=4)
+
    # Cities
    cities = get_cities_gdf().to_crs(3857)
    plot_cities(ax, cities)
+    # County labels
+    label_selected_counties(ax, counties)
+    # Bold borders for target counties
+    bold_target_county_borders(ax, counties)

    # Scale bar and north arrow
    ax.add_artist(ScaleBar(1, units='m', location='lower left', box_alpha=0.7))