Fill U.S. counties by value using Python & Cartopy?

Question

I'd like to know how to fill in a map of U.S. counties by value (i.e., a chloropleth map), using Python 3 and Cartopy, and I haven't yet found anything online to guide me in that. That filled value could be, for instance, highest recorded tornado rating (with counties left blank for no recorded tornadoes), or even something arbitrary such as whether I've visited (=1) or lived (=2) in the county. I found a helpful MetPy example to get the county boundaries on a map:

https://unidata.github.io/MetPy/latest/examples/plots/US_Counties.html

What I envision is somehow setting a list (or dictionary?) of county names to a certain value, and then each value would be assigned to a particular fill color. This is my current script, which generates a nice blank county map of the CONUS/lower 48 (though I'd eventually also like to add Alaska/Hawaii insets).

import cartopy
import cartopy.crs as ccrs
import matplotlib as mpl
import matplotlib.pyplot as plt
from metpy.plots import USCOUNTIES

plot_type = 'png'

borders  = cartopy.feature.BORDERS
states   = cartopy.feature.NaturalEarthFeature(category='cultural', scale='10m', facecolor='none', name='admin_1_states_provinces_lakes')
oceans   = cartopy.feature.OCEAN
lakes    = cartopy.feature.LAKES

mpl.rcParams['figure.figsize'] = (12,10)
water_color = 'lightblue'

fig = plt.figure()
ax = plt.axes(projection=ccrs.LambertConformal(central_longitude=-97.5, central_latitude=38.5, standard_parallels=(38.5,38.5)))
ax.set_extent([-120, -74, 23, 50], ccrs.Geodetic())
ax.coastlines()
ax.add_feature(borders, linestyle='-')
ax.add_feature(states, linewidth=0.50, edgecolor='black')
ax.add_feature(oceans, facecolor=water_color)
ax.add_feature(lakes, facecolor=water_color, linewidth=0.50, edgecolor='black')
ax.add_feature(USCOUNTIES.with_scale('500k'), linewidth=0.10, edgecolor='black')
plt.savefig('./county_map.'+plot_type)
plt.close()

Any ideas or tips on how to assign values to counties and fill them accordingly?

Answer 1

So Cartopy's shapereader.Reader can give you access to all of the records in the shapefile, including their attributes. Putting this together with MetPy's get_test_data to get access to the underlying shapefile you can get what you want, assuming you have a dataset that maps e.g. FIPSCODE to EF rating:

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())

cmap = plt.get_cmap('magma')
norm = plt.Normalize(0, 5)

# Fake tornado dataset with a value for each county code
tor_data = dict()

# This will only work (have access to the shapefile's database of
# attributes after it's been download by using `USCOUNTIES` or
# running get_test_data() for the .shx and .dbf files as well.
for rec in shpreader.Reader(get_test_data('us_counties_20m.shp',
                                          as_file_obj=False)).records():
    # Mimic getting data, but actually getting a random number
    # GEOID seems to be the FIPS code
    max_ef = tor_data.get(rec.attributes['GEOID'], np.random.randint(0, 5))

    # Normalize the data to [0, 1] and colormap manually
    color = tuple(cmap(norm(max_ef)))

    # Add the geometry to the plot, being sure to specify the coordinate system
    ax.add_geometries([rec.geometry], crs=ccrs.PlateCarree(), facecolor=color)

ax.set_extent((-125, -65, 25, 48))

That gives me: