How do you plot US county borders in Cartopy?
It's very straight forward to plot state and country boundaries
ax.add_feature(cfeature.BORDERS.with_scale('50m'))
ax.add_feature(cfeature.STATES.with_scale('50m'))
But I can't seem to find a similar method to add county boundaries. This was one of the nice things about Basemap.
Cartopy is a Python package designed for geospatial data processing in order to produce maps and other geospatial data analyses. Cartopy makes use of the powerful PROJ, NumPy and Shapely libraries and includes a programmatic interface built on top of Matplotlib for the creation of publication quality maps.
Given cartopy's ability to draw shapefiles, this question essentially boils down to "where can I find US county outlines?".
A similar question was asked on the Natural Earth forum at http://www.naturalearthdata.com/forums/topic/u-s-county-shape-file/. It pointed to a location at http://nationalatlas.gov/mld/countyp.html, which has unfortunately had a dose of bit-rot. A quick Google suggests this can now be found at:
https://nationalmap.gov/small_scale/atlasftp.html?openChapters=chpbound#chpbound
I decided to download one of those county shapefiles:
https://prd-tnm.s3.amazonaws.com/StagedProducts/Small-scale/data/Boundaries/countyl010g_shp_nt00964.tar.gz
With that in place, I used cartopy's shapereader to get the geometries out, and created a custom feature that can then be added to the axes:
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.io.shapereader as shpreader
import matplotlib.pyplot as plt
reader = shpreader.Reader('countyl010g.shp')
counties = list(reader.geometries())
COUNTIES = cfeature.ShapelyFeature(counties, ccrs.PlateCarree())
plt.figure(figsize=(10, 6))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.add_feature(cfeature.LAND.with_scale('50m'))
ax.add_feature(cfeature.OCEAN.with_scale('50m'))
ax.add_feature(cfeature.LAKES.with_scale('50m'))
ax.add_feature(COUNTIES, facecolor='none', edgecolor='gray')
ax.coastlines('50m')
ax.set_extent([-83, -65, 33, 44])
plt.show()
This is derived from the example at https://scitools.org.uk/cartopy/docs/v0.14/examples/feature_creation.html which constructs a NaturalEarthFeature
, rather than a ShapelyFeature
, but otherwise, the principle is pretty much the same.
Hope that is useful to you.
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With