Plotting Colorized Map of Specific US Counties with Cartopy - Python

Question

I am attempting to plot colored US counties based on some data, using Cartopy. I have created a dataset of FIPS codes and the count of 1950-2022 tornadoes for the state of Michigan. I have it in a pandas dataframe and a dict (below).

I assume I need to somehow compare my FIPS codes to the MetPy USCOUNTIES set, to use the geometries stored in USCOUNTIES?

Any ideas?

My data looks like this as a dict:

{
    '26049': 44,
    '26115': 33,
    '26081': 32,
    '26163': 31,
    '26021': 30,
    '26005': 30,
    '26125': 28,
    '26155': 28,
    '26091': 28,
    '26161': 26,
    '26077': 25,
    '26059': 23,
    '26045': 23,
    '26145': 23,
    '26093': 22,
    '26067': 22,
    '26065': 22,
    '26099': 22,
    '26147': 21,
    '26139': 20,
    '26151': 18,
    '26037': 18,
    '26015': 17,
    '26157': 17,
    '26063': 16,
    '26159': 16,
    '26075': 16,
    '26087': 16,
    '26023': 15,
    '26129': 15,
    '26027': 14,
    '26017': 14,
    '26025': 13,
    '26007': 13,
    '26057': 13,
    '26073': 12,
    '26123': 12,
    '26133': 12,
    '26009': 11,
    '26041': 11,
    '26149': 11,
    '26111': 10,
    '26039': 10,
    '26069': 10,
    '26107': 10,
    '26103': 9,
    '26143': 9,
    '26051': 9,
    '26071': 8,
    '26035': 8,
    '26011': 8,
    '26001': 8,
    '26121': 8,
    '26165': 8,
    '26117': 8,
    '26079': 7,
    '26043': 7,
    '26113': 7,
    '26141': 7,
    '26003': 6,
    '26031': 6,
    '26119': 6,
    '26033': 6,
    '26109': 5,
    '26127': 5,
    '26047': 5,
    '26105': 5,
    '26137': 4,
    '26095': 4,
    '26053': 4,
    '26135': 4,
    '26085': 4,
    '26055': 3,
    '26131': 3,
    '26029': 3,
    '26089': 3,
    '26097': 3,
    '26019': 3,
    '26083': 2,
    '26153': 2,
    '26013': 2,
    '26101': 2,
    '26000': 1
}

Here is a random Michigan county map with COVID-19 data I found on Google that looks like what I want, counties in Michigan colored based on my tornado count data. Ideally, this code would only plot the FIPS codes that I have in my dataset, not every US county.

enter image description here

Answer 1

Cartopy's add_geometries() method can take a styler argument, which is a function that returns the visual style arguments for a given piece of geometry. You can see it demonstrated in this example. Of course in this case, it's not obvious how to make that. In the example below, I took your dict above and called that tor_counts. I then looped over the records in the shapefile to generate a lookup from geometry to tornado count:

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.io.shapereader as shpreader
import matplotlib.pyplot as plt
from metpy.cbook import get_test_data

# Open MetPy's counties shapefile using Cartopy's shape reader
counties = shpreader.Reader(get_test_data('us_counties_20m.shp', as_file_obj=False))

# Loop over the records in the shapefile and generate a dictionary that
# maps geometry->tornado count
county_lookup = {rec.geometry: tor_counts.get(rec.attributes['GEOID'], 0)
                 for rec in counties.records() if rec.attributes['STATEFP'] == '26'}

# Create a "styler" that returns the appropriate draw colors, etc.
# for a given geometry. This uses our lookup to find the counties
# that need coloring
def color_torns(geom):
    count = county_lookup.get(geom, 0)
    if count:
        cmap = plt.get_cmap('viridis')
        norm = plt.Normalize(0, 50)
        facecolor = cmap(norm(count))
    else:
        facecolor = 'none'
    return {'edgecolor': 'black', 'facecolor': facecolor}


# Create figure with some maps
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(projection=ccrs.LambertConformal())
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.STATES)

# Add the geometries from the shapefile, pass our styler function
ax.add_geometries(counties.geometries(), ccrs.PlateCarree(), styler=color_torns)
ax.set_extent((-91, -82, 41, 48))

The result:

I should add that checking the state FIPS (of 26) above is just a way to keep the dictionary smaller than having one entry for each county in the US since we really don't care. This would pretty much scale to that though.

Answer 2

Another option where you only need to add the dict to geopandas dataframe https://geopandas.org/en/stable/docs/user_guide/mapping.html