diff --git a/colorado_shape/USA_Census_Counties_7412386151545161096/dtl_cnty.cpg b/colorado_shape/USA_Census_Counties_7412386151545161096/dtl_cnty.cpg new file mode 100644 index 0000000..3ad133c --- /dev/null +++ b/colorado_shape/USA_Census_Counties_7412386151545161096/dtl_cnty.cpg @@ -0,0 +1 @@ +UTF-8 \ No newline at end of file diff --git a/colorado_shape/USA_Census_Counties_7412386151545161096/dtl_cnty.dbf b/colorado_shape/USA_Census_Counties_7412386151545161096/dtl_cnty.dbf new file mode 100644 index 0000000..e00a766 Binary files /dev/null and b/colorado_shape/USA_Census_Counties_7412386151545161096/dtl_cnty.dbf differ diff --git a/colorado_shape/USA_Census_Counties_7412386151545161096/dtl_cnty.prj b/colorado_shape/USA_Census_Counties_7412386151545161096/dtl_cnty.prj new file mode 100644 index 0000000..f45cbad --- /dev/null +++ b/colorado_shape/USA_Census_Counties_7412386151545161096/dtl_cnty.prj @@ -0,0 +1 @@ +GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]] \ No newline at end of file diff --git a/colorado_shape/USA_Census_Counties_7412386151545161096/dtl_cnty.shp b/colorado_shape/USA_Census_Counties_7412386151545161096/dtl_cnty.shp new file mode 100644 index 0000000..5a3a6a1 Binary files /dev/null and b/colorado_shape/USA_Census_Counties_7412386151545161096/dtl_cnty.shp differ diff --git a/colorado_shape/USA_Census_Counties_7412386151545161096/dtl_cnty.shp.xml b/colorado_shape/USA_Census_Counties_7412386151545161096/dtl_cnty.shp.xml new file mode 100644 index 0000000..15fe718 --- /dev/null +++ b/colorado_shape/USA_Census_Counties_7412386151545161096/dtl_cnty.shp.xml @@ -0,0 +1,2 @@ + +20241126211117001.0FGDC CSDGM MetadataFALSEwork2002GeographicGCS_WGS_1984Angular Unit: Degree (0.017453)<GeographicCoordinateSystem xsi:type='typens:GeographicCoordinateSystem' xmlns:xsi='http://www.w3.org/2001/XMLSchema-instance' xmlns:xs='http://www.w3.org/2001/XMLSchema' xmlns:typens='http://www.esri.com/schemas/ArcGIS/10.7'><WKT>GEOGCS[&quot;GCS_WGS_1984&quot;,DATUM[&quot;D_WGS_1984&quot;,SPHEROID[&quot;WGS_1984&quot;,6378137.0,298.257223563]],PRIMEM[&quot;Greenwich&quot;,0.0],UNIT[&quot;Degree&quot;,0.0174532925199433],AUTHORITY[&quot;EPSG&quot;,4326]]</WKT><XOrigin>-400</XOrigin><YOrigin>-400</YOrigin><XYScale>999999999.99999988</XYScale><ZOrigin>-100000</ZOrigin><ZScale>10000</ZScale><MOrigin>-100000</MOrigin><MScale>10000</MScale><XYTolerance>8.983152841195215e-09</XYTolerance><ZTolerance>0.001</ZTolerance><MTolerance>0.001</MTolerance><HighPrecision>true</HighPrecision><LeftLongitude>-180</LeftLongitude><WKID>4326</WKID><LatestWKID>4326</LatestWKID></GeographicCoordinateSystem>20220329163624002022032916362400 Version 6.2 (Build 9200) ; Esri ArcGIS 10.7.1.11595work2This layer presents the 2020 U.S. Census County (or County Equivalent) boundaries of the United States in the 50 states and the District of Columbia. It includes 2020 U.S. Census codes and population information. U.S. Counties are polygons containing population totals from the 2020 Census.<div style="font-family:&quot;Avenir Next W01&quot;, &quot;Avenir Next W00&quot;, &quot;Avenir Next&quot;, Avenir, &quot;Helvetica Neue&quot;, sans-serif; font-size:16px;"><span style="font-family:inherit;">This layer presents the 2020 U.S. Census County (or County Equivalent) boundaries of the United States in the 50 states and the District of Columbia.</span></div><div style="font-family:&quot;Avenir Next W01&quot;, &quot;Avenir Next W00&quot;, &quot;Avenir Next&quot;, Avenir, &quot;Helvetica Neue&quot;, sans-serif; font-size:16px;"><span style="font-family:inherit;"><br /></span></div><div style="font-family:&quot;Avenir Next W01&quot;, &quot;Avenir Next W00&quot;, &quot;Avenir Next&quot;, Avenir, &quot;Helvetica Neue&quot;, sans-serif; font-size:16px;">This layer is updated annually. <span style="font-family:inherit;">The geography is sourced from <a href="https://www.census.gov/geographies/mapping-files/time-series/geo/tiger-geodatabase-file.2020.html" style="color:rgb(0, 97, 155); text-decoration-line:none; font-family:inherit;" target="_blank">U.S. Census Bureau 2020 TIGER FGDB (National Sub-State)</a> and edited using <a href="https://www.census.gov/geographies/mapping-files/time-series/geo/tiger-geodatabase-file.2020.html" style="color:rgb(0, 97, 155); text-decoration-line:none; font-family:inherit;" target="_blank">TIGER Hydrography</a> to add a detailed coastline for cartographic purposes. </span>Attribute fields include 2020 total population from the <a href="https://www.census.gov/data/datasets/2020/dec/2020-census-redistricting-summary-file-dataset.html" style="color:rgb(0, 97, 155); text-decoration-line:none; font-family:inherit;" target="_blank">U.S. Census Public Law 94 data</a>.</div><span style="font-family:&quot;Avenir Next W01&quot;, &quot;Avenir Next W00&quot;, &quot;Avenir Next&quot;, Avenir, &quot;Helvetica Neue&quot;, sans-serif; font-size:16px;"><br /></span><span style="font-family:&quot;Avenir Next W01&quot;, &quot;Avenir Next W00&quot;, &quot;Avenir Next&quot;, Avenir, &quot;Helvetica Neue&quot;, sans-serif; font-size:16px;">This ready-to-use layer can be used in ArcGIS Pro and in ArcGIS Online and its configurable apps, dashboards, StoryMaps, custom apps, and mobile apps. The data can also be exported for offline workflows. Cite the 'U.S. Census Bureau' when using this data.</span>Sources: Esri; U.S. Department of Commerce, Census Bureau; U.S. Department of Commerce (DOC), National Oceanic and Atmospheric Administration (NOAA), National Ocean Service (NOS), National Geodetic Survey (NGS)File Geodatabase Feature ClassdatasetEPSG6.14(3.0.1)0SimpleFALSE0TRUEFALSEwork2Feature Class0OBJECTIDOBJECTIDOID400Internal feature number.EsriSequential unique whole numbers that are automatically generated.ShapeShapeGeometry000Feature geometry.EsriCoordinates defining the features.tmpfipstmpfipsString500NAMENAMEString5000STATE_NAMESTATE_NAMEString2000STATE_ABBRSTATE_ABBRString200STATE_FIPSSTATE_FIPSString200COUNTY_FIPSCOUNTY_FIPSString300FIPSFIPSString500POPULATIONPOPULATION_2021Integer400POP_SQMIPOP21_SQMIDouble800POPULATION_2020POPULATION_2020Integer400POP20_SQMIPOP20_SQMIDouble800SQMIwork1.SQMIDouble800Shape_LengthShape_LengthDouble800Length of feature in internal units.EsriPositive real numbers that are automatically generated.Shape_AreaShape_AreaDouble800Area of feature in internal units squared.EsriPositive real numbers that are automatically generated.20220329 diff --git a/colorado_shape/USA_Census_Counties_7412386151545161096/dtl_cnty.shx b/colorado_shape/USA_Census_Counties_7412386151545161096/dtl_cnty.shx new file mode 100644 index 0000000..9920d62 Binary files /dev/null and b/colorado_shape/USA_Census_Counties_7412386151545161096/dtl_cnty.shx differ diff --git a/out_tables/spill_map_hex.pdf b/out_tables/spill_map_hex.pdf index 4f18e0a..8951cde 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 new file mode 100644 index 0000000..8cab549 Binary files /dev/null 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 3d033f2..fe36b68 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 3f5accd..c805445 100644 --- a/scripts/make_spill_map.py +++ b/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