Introdcution
Working with geographic data often requires converting latitude and longitude coordinates into administrative regions such as counties or states. This process—called reverse geocoding—is particularly useful for applications like wildfire analysis, environmental monitoring, satellite data validation, and geospatial analytics.
The workflow uses open datasets from the US Census Bureau.
Download US County Boundary Data
The most authoritative source for US county boundaries is the TIGER/Line dataset from the US Census Bureau.
Download the County shapefile:
https://www.census.gov/geographies/mapping-files/time-series/geo/tiger-line-file.html
Look for (Counties):
1 | TIGER/Line Shapefiles → Counties |
Download the nationwide dataset:
1 | tl_2023_us_county.zip |
After extracting the archive you should see files like:
1 2 3 4 | tl_2023_us_county.shp tl_2023_us_county.dbf tl_2023_us_county.shx tl_2023_us_county.prj |
Install Required Python Libraries
1 | pip install geopandas shapely folium pandas |
Load the County Shapefile
1 2 3 | import geopandas as gpd counties = gpd.read_file("tl_2023_us_county.shp") |
1 | counties.columns |
returns
1 2 3 4 5 | Index(['STATEFP', 'COUNTYFP', 'COUNTYNS', 'GEOID', 'GEOIDFQ', 'NAME', 'NAMELSAD', 'LSAD', 'CLASSFP', 'MTFCC', 'CSAFP', 'CBSAFP', 'METDIVFP', 'FUNCSTAT', 'ALAND', 'AWATER', 'INTPTLAT', 'INTPTLON', 'geometry', 'state_name'], dtype='object') |
Convert State FIPS Codes to State Names
STATEFP is the state FIPS code (a numeric identifier used by the US Census).
To get state names, you simply map the FIPS codes to a state lookup table.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 | state_fips_to_name = { "01": "Alabama", "02": "Alaska", "04": "Arizona", "05": "Arkansas", "06": "California", "08": "Colorado", "09": "Connecticut", "10": "Delaware", "11": "District of Columbia", "12": "Florida", "13": "Georgia", "15": "Hawaii", "16": "Idaho", "17": "Illinois", "18": "Indiana", "19": "Iowa", "20": "Kansas", "21": "Kentucky", "22": "Louisiana", "23": "Maine", "24": "Maryland", "25": "Massachusetts", "26": "Michigan", "27": "Minnesota", "28": "Mississippi", "29": "Missouri", "30": "Montana", "31": "Nebraska", "32": "Nevada", "33": "New Hampshire", "34": "New Jersey", "35": "New Mexico", "36": "New York", "37": "North Carolina", "38": "North Dakota", "39": "Ohio", "40": "Oklahoma", "41": "Oregon", "42": "Pennsylvania", "44": "Rhode Island", "45": "South Carolina", "46": "South Dakota", "47": "Tennessee", "48": "Texas", "49": "Utah", "50": "Vermont", "51": "Virginia", "53": "Washington", "54": "West Virginia", "55": "Wisconsin", "56": "Wyoming" } counties["state_name"] = counties["STATEFP"].map(state_fips_to_name) |
Get All Counties for a Given State
Once you have loaded the TIGER/Line counties GeoDataFrame, extracting all counties for a given state is very straightforward. Each county is associated with a state FIPS code (STATEFP), which provides a fast and reliable way to filter the data.
For example, Maryland corresponds to the FIPS code "24":
1 | maryland_counties = counties[counties["STATEFP"] == "24"] |
If you previously mapped FIPS codes to state names, you can also use:
1 | counties[counties["state_name"] == "Maryland"] |
Both approaches are equivalent, but filtering by STATEFP is faster and recommended, especially when working with large datasets such as satellite observations.
Extract County Names
Once filtered, retrieving the list of county names is simple:
1 2 | county_names = maryland_counties["NAME"].tolist() print(county_names) |
This returns:
1 2 3 4 5 | ['Harford', 'Garrett', 'Carroll', "St. Mary's", 'Caroline', 'Calvert', 'Somerset', 'Baltimore', 'Anne Arundel', 'Kent', "Prince George's", 'Dorchester', 'Montgomery', "Queen Anne's", 'Howard', 'Allegany', 'Baltimore', 'Charles', 'Frederick', 'Wicomico', 'Worcester', 'Cecil', 'Talbot', 'Washington'] |
Note: You may notice "Baltimore" appears twice. This is expected:
- Baltimore County
- Baltimore City (independent city treated as a county equivalent)
Visualize Counties with Folium
Interactive visualization is extremely useful for validating geographic datasets, especially when working with satellite-derived data such as fire detections.
Example: Plot Howard County
We first extract a single county:
1 2 3 4 | gdf_sample = counties[ (counties["STATEFP"] == "24") & (counties["NAME"] == "Howard") ] |
Convert to geographic coordinates (required by Folium):
1 | gdf_sample = gdf_sample.to_crs("EPSG:4326") |
Compute the map center:
1 2 3 4 | center = [ gdf_sample.geometry.centroid.y.mean(), gdf_sample.geometry.centroid.x.mean() ] |
Create the interactive map:
1 2 3 4 5 6 7 8 | import folium m = folium.Map( location=center, zoom_start=11, tiles='https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}', attr='Esri World Imagery' ) |
Add the county boundary:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 | tooltip_fields = [col for col in gdf_sample.columns if col != "geometry"] tooltip_fields = tooltip_fields[:5] folium.GeoJson( gdf_sample, name="Howard County", style_function=lambda x: { "color": "yellow", "weight": 3, "fillOpacity": 0.2 }, tooltip=folium.GeoJsonTooltip(fields=tooltip_fields) ).add_to(m) folium.LayerControl().add_to(m) m |
Visualize All Maryland Counties with Labels
We can extend this approach to visualize all counties in Maryland with labels.
First, compute centroids:
1 2 | maryland_counties = maryland_counties.to_crs("EPSG:4326") maryland_counties["centroid"] = maryland_counties.geometry.centroid |
Create the map:
1 2 3 4 5 6 7 8 9 10 11 | center = [ maryland_counties.centroid.y.mean(), maryland_counties.centroid.x.mean() ] m = folium.Map( location=center, zoom_start=8, tiles='https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}', attr='Esri World Imagery' ) |
Add county polygons:
1 2 3 4 5 6 7 8 9 10 | folium.GeoJson( maryland_counties.drop(columns=["centroid"]), name="Maryland Counties", style_function=lambda x: { "color": "yellow", "weight": 2, "fillOpacity": 0.15 }, tooltip=folium.GeoJsonTooltip(fields=["NAME"]) ).add_to(m) |
Add labels:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 | for _, row in maryland_counties.iterrows(): folium.Marker( location=[row.centroid.y, row.centroid.x], icon=folium.DivIcon( html=f""" <div style=" font-size:10px; font-weight:bold; color:white; text-shadow:1px 1px 2px black; "> {row['NAME']} </div> """ ) ).add_to(m) |
Finalize:
1 2 | folium.LayerControl().add_to(m) m |
Tip: For better label placement in irregular shapes, consider using:
1 | row.geometry.representative_point() |
instead of centroids.
Retrieve County Name for a Specific Latitude/Longitude
Single Coordinate Lookup
1 2 3 4 5 6 7 8 9 10 | from shapely.geometry import Point lat = 39.2037 lon = -76.8610 point = Point(lon, lat) county = counties[counties.contains(point)] print(county["NAME"].values) |
Output:
1 | ['Howard'] |
This method is ideal for quick lookups.
Spatial Join for Large Datasets
For large datasets (e.g., thousands or millions of points), use a spatial join:
1 2 3 4 5 6 7 | points = gpd.GeoDataFrame( df, geometry=gpd.points_from_xy(df.lon, df.lat), crs="EPSG:4326" ) result = gpd.sjoin(points, counties, predicate="within") |
This efficiently assigns a county to each point.
Merge County Data with Other Datasets
You can enrich county data by merging it with other geographic datasets such as state boundaries:
1 2 3 4 5 | counties_with_state = gpd.sjoin( counties, states[["name", "geometry"]], predicate="within" ) |
This allows you to combine:
- County-level data
- State-level attributes
- External datasets (e.g., emissions, fire detections)
References
| Links | Site |
|---|---|
| https://www.census.gov/geographies/mapping-files/time-series/geo/tiger-line-file.html | US Census TIGER/Line Shapefiles |
| https://geopandas.org | GeoPandas Documentation |
| https://python-visualization.github.io/folium/ | Folium Documentation |
| https://shapely.readthedocs.io | Shapely Documentation |
| https://naturalearthdata.com |
