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?

score 2 · Answer 1 · answered Nov 24 '21 at 00:20

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:

sin tribu · Answer 2 · 2021-11-21T07:03:58.123

I'm not sure about passing in a dict, but you can pass in a list to facecolor.

ax.add_feature(USCOUNTIES.with_scale('500k'), linewidth=0.10, edgecolor='black', facecolor=["red", "blue", "green"])

If you know how many counties there are you can make a list that long by:

import matplotlib.cm as cm
import numpy as np 

number_of_counties = 3000
color_scale = list(cm.rainbow(np.linspace(0, 1, number_of_counties)))

ax.add_feature(USCOUNTIES.with_scale('500k'), linewidth=.10, edgecolor="black", facecolor=color_scale)

but they didn't make it easy to extract the names from USCOUNTIES. You can see where it is defined in your source code:

from metpy import plots
print(plots.__file__)

If you go inside the directory printed there is a file named cartopy_utils.py and inside the class definition for class MetPyMapFeature(Feature): you will see USCOUNTIES. You might have better luck than I did mapping county names to the geometric shapes.

EDIT: Also, I just used cm.rainbow as an example, you can choose from any color map https://matplotlib.org/stable/tutorials/colors/colormaps.html. Not sure if it even goes up to 3000, but you get the idea.

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

2 Answers2