maps

scripts/make_spill_map.py

@@ -9,11 +9,12 @@ import os
 from pathlib import Path
 import pandas as pd
 import geopandas as gpd
-from shapely.geometry import Point
+from shapely.geometry import Point, box
 import matplotlib.pyplot as plt
 from matplotlib_scalebar.scalebar import ScaleBar
 import matplotlib.patheffects as pe
 from sqlalchemy import create_engine
+from matplotlib import patches
 
 ROOT = Path('/home/dadams/Repos/colorado_spills')
 ANALYSIS_DIR = ROOT / 'analysis' / 'new analysis Aug 2025'
@@ -127,9 +128,42 @@ def try_query_counties(engine) -> gpd.GeoDataFrame | None:
 
 
 def get_colorado_boundary(spills_gdf: gpd.GeoDataFrame) -> tuple[gpd.GeoDataFrame, gpd.GeoDataFrame | None]:
-    """Return Colorado boundary and optional counties. Prefer PostGIS; fallback to convex-hull from points."""
-    engine = get_engine()
+    """Return Colorado boundary and optional counties. Prefer local shapefile, then PostGIS; fallback to convex-hull from points."""
     counties = None
+
+    # 1) Prefer local shapefile if present
+    shp_dir = ROOT / 'colorado_shape' / 'USA_Census_Counties_7412386151545161096'
+    shp_path = shp_dir / 'dtl_cnty.shp'
+    if shp_path.exists():
+        try:
+            gdf = gpd.read_file(shp_path)
+            # Try common Colorado filters
+            col_candidates = [
+                ("STATEFP", "08"),
+                ("STATE_FIPS", "08"),
+                ("STATE", "CO"),
+                ("STATE_NAME", "Colorado"),
+            ]
+            cg = gdf
+            for col, val in col_candidates:
+                if col in gdf.columns:
+                    if val == "08":
+                        cg = gdf[gdf[col].astype(str) == "08"]
+                    else:
+                        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)
+            counties = cg.to_crs(4326)
+            # Dissolve to single state boundary
+            boundary = counties.dissolve().reset_index(drop=True)
+            return boundary[["geometry"]], counties[["geometry"]]
+        except Exception:
+            pass
+
+    # 2) PostGIS
+    engine = get_engine()
     if engine is not None:
         boundary = try_query_boundary(engine)
         if boundary is not None and not boundary.empty:
@@ -155,6 +189,36 @@ def style_axes(ax):
     for spine in ax.spines.values():
         spine.set_visible(False)
 
+def get_cities_gdf():
+    """Major Colorado cities for reference with lat/lon (EPSG:4326)."""
+    cities = [
+        {"name": "Denver", "lat": 39.7392, "lon": -104.9903},
+        {"name": "Colorado Springs", "lat": 38.8339, "lon": -104.8214},
+        {"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},
+        {"name": "Boulder", "lat": 40.01499, "lon": -105.2705},
+        {"name": "Durango", "lat": 37.2753, "lon": -107.8801},
+    ]
+    df = pd.DataFrame(cities)
+    gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df['lon'], df['lat']), crs='EPSG:4326')
+    return gdf
+
+def plot_cities(ax, cities_gdf_3857):
+    """Plot city points and labels on given axis (coordinates in EPSG:3857)."""
+    x = cities_gdf_3857.geometry.x.values
+    y = cities_gdf_3857.geometry.y.values
+    ax.scatter(x, y, s=12, facecolor='white', edgecolor='black', linewidth=0.6, zorder=6)
+    for i, row in cities_gdf_3857.iterrows():
+        ax.annotate(
+            row['name'],
+            (row.geometry.x, row.geometry.y),
+            xytext=(4, 3), textcoords='offset points',
+            fontsize=8, color='black', zorder=6,
+            path_effects=[pe.withStroke(linewidth=2.5, foreground='white')]
+        )
+
 
 def map_points(spills: gpd.GeoDataFrame, boundary: gpd.GeoDataFrame, counties: gpd.GeoDataFrame | None):
     # Project to Web Mercator for plotting
@@ -171,6 +235,10 @@ def map_points(spills: gpd.GeoDataFrame, boundary: gpd.GeoDataFrame, counties: g
     # Plot points with alpha and small size to avoid overplotting
     spills_3857.plot(ax=ax, markersize=1.2, alpha=0.35, color='#2c7fb8')
 
+    # Cities
+    cities = get_cities_gdf().to_crs(3857)
+    plot_cities(ax, cities)
+
     # Tighten to boundary
     minx, miny, maxx, maxy = boundary_3857.total_bounds
     ax.set_xlim(minx, maxx)
@@ -193,8 +261,39 @@ def map_points(spills: gpd.GeoDataFrame, boundary: gpd.GeoDataFrame, counties: g
             pass
     ax.set_title(title, fontsize=12, pad=10)
 
+    # Front Range inset (approx bounding box Denver–Boulder–Fort Collins corridor)
+    # Define bbox in EPSG:4326 and project to 3857
+    fr_minx, fr_miny, fr_maxx, fr_maxy = -106.0, 39.0, -104.5, 40.9
+    bbox_poly_3857 = gpd.GeoSeries([box(fr_minx, fr_miny, fr_maxx, fr_maxy)], crs='EPSG:4326').to_crs(3857)
+    bxmin, bymin, bxmax, bymax = bbox_poly_3857.total_bounds
+
+    # Locator rectangle on main map
+    rect = patches.Rectangle((bxmin, bymin), bxmax - bxmin, bymax - bymin,
+                             linewidth=1.2, edgecolor='#333333', facecolor='none', zorder=7)
+    ax.add_patch(rect)
+
+    # Inset axes: [left, bottom, width, height] in figure fraction
+    inset_ax = fig.add_axes([0.62, 0.62, 0.32, 0.32])
+    boundary_3857.plot(ax=inset_ax, color='#f5f5f5', edgecolor='#888888', linewidth=0.6)
+    if counties is not None and not counties.empty:
+        counties_inset_3857 = counties.to_crs(3857)
+        counties_inset_3857.boundary.plot(ax=inset_ax, color='#aaaaaa', linewidth=0.25)
+    spills_3857.plot(ax=inset_ax, markersize=1.2, alpha=0.35, color='#2c7fb8')
+    # Cities on inset
+    cities_inset = get_cities_gdf().to_crs(3857)
+    plot_cities(inset_ax, cities_inset)
+    inset_ax.set_xlim(bxmin, bxmax)
+    inset_ax.set_ylim(bymin, bymax)
+    style_axes(inset_ax)
+    # Thin frame around inset
+    for spine in inset_ax.spines.values():
+        spine.set_visible(True)
+        spine.set_linewidth(0.6)
+        spine.set_edgecolor('#444444')
+
     out_png = OUT / 'spill_map_points.png'
-    fig.savefig(out_png, dpi=300)
+    # Tight crop to reduce whitespace
+    fig.savefig(out_png, dpi=300, bbox_inches='tight', pad_inches=0.05)
     plt.close(fig)
     return out_png
 
@@ -219,6 +318,10 @@ def map_hex(spills: gpd.GeoDataFrame, boundary: gpd.GeoDataFrame, counties: gpd.
     cb = fig.colorbar(hb, ax=ax, fraction=0.028, pad=0.01)
     cb.set_label('Spill count')
 
+    # Cities
+    cities = get_cities_gdf().to_crs(3857)
+    plot_cities(ax, cities)
+
     # Scale bar and north arrow
     ax.add_artist(ScaleBar(1, units='m', location='lower left', box_alpha=0.7))
     add_north_arrow(ax)
@@ -228,20 +331,12 @@ def map_hex(spills: gpd.GeoDataFrame, boundary: gpd.GeoDataFrame, counties: gpd.
     ax.set_xlim(minx, maxx)
     ax.set_ylim(miny, maxy)
     style_axes(ax)
-    # Dynamic title from available year fields
-    title = 'Spill Density (Hexbin)'
-    year_cols = [c for c in spills.columns if c.lower() in ('year', 'yr')]
-    if year_cols:
-        try:
-            yrs = spills[year_cols[0]].dropna().astype(int)
-            if not yrs.empty:
-                title = f"Spill Density (Hexbin), {yrs.min()}–{yrs.max()}"
-        except Exception:
-            pass
-    ax.set_title(title, fontsize=12, pad=10)
+    # No title per manuscript layout; keep figure clean
 
     out_pdf = OUT / 'spill_map_hex.pdf'
+    out_png = OUT / 'spill_map_hex.png'
     fig.savefig(out_pdf)
+    fig.savefig(out_png, dpi=300, bbox_inches='tight', pad_inches=0.05)
     plt.close(fig)
     return out_pdf