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.
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.