diff --git a/out_tables/spill_map_hex.pdf b/out_tables/spill_map_hex.pdf index f8d39b3..ea62a43 100644 Binary files a/out_tables/spill_map_hex.pdf and b/out_tables/spill_map_hex.pdf differ diff --git a/out_tables/spill_map_hex.png b/out_tables/spill_map_hex.png index bdbe633..c9a40b8 100644 Binary files a/out_tables/spill_map_hex.png and b/out_tables/spill_map_hex.png differ diff --git a/out_tables/spill_map_points.png b/out_tables/spill_map_points.png index fe36b68..6df22d0 100644 Binary files a/out_tables/spill_map_points.png and b/out_tables/spill_map_points.png differ diff --git a/scripts/make_spill_map.py b/scripts/make_spill_map.py index 67eb016..8efef73 100644 --- 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))