I am looking to add a solid fill underneath my contoured data (contourf in Matplotlib). Specifically, I would like the solid fill to be constrained to specific shapefiles, in this case individual states in the continental United States.
My current plot
My Current Code
Please note that I've included the key inner workings to plot the contours and shapefiles, I've excluded extraneous code (like the city plots) to help those working on a solution to this question.
import os, sys, re, ast, tifffile
import matplotlib.pyplot as plt
import matplotlib.colors as clr
import matplotlib.figure as fig
import matplotlib as mpl
import cartopy.crs as crs
import cartopy.feature as cfeature
from matplotlib.offsetbox import AnnotationBbox, OffsetImage
from cartopy.io.shapereader import Reader
from cartopy.feature import ShapelyFeature
states = '/filepath/gis/st_us.shp'
counties = '/filepath/gis/cnt_us.shp'
lakes_hires = '/filepath/gis/lk_us.shp'
states_feature = ShapelyFeature(Reader(states).geometries(), crs.LambertConformal(), facecolor='none', edgecolor='black')
states_fill = ShapelyFeature(Reader(states).geometries(), crs.LambertConformal(), facecolor='#C7C6C1', edgecolor='none')
counties_feature = ShapelyFeature(Reader(counties).geometries(), crs.LambertConformal(), facecolor='none', edgecolor='lightgray')
lakeshi_feature = ShapelyFeature(Reader(lakes_hires).geometries(), crs.LambertConformal(), facecolor='none', edgecolor='navy')
file = '/filepath/gis/globalpop_conus_1km.tif'
tiffimage = tifffile.imread(file)
lat0 = 52.817021799999999
lon0 = -130.122944799999999
dlat = 0.008332904595765
dlon = 0.008332442558080
tiffimage[tiffimage<-1e100]=np.nan
tiffimage[tiffimage>1e100]=np.nan
tiflats = np.linspace(lat0,(lat0 - (tiffimage.shape[0]*dlat)),num=tiffimage.shape[0])
tiflons = np.linspace(lon0,(lon0 + (tiffimage.shape[1]*dlon)),num=tiffimage.shape[1])
tiffgrid = np.meshgrid(tiflats,tiflons)
tiffimage = tiffimage.transpose()
xtif,ytif = tiffgrid[1],tiffgrid[0]
fig = plt.figure(figsize=(14,8))
ax = plt.axes([0.25, 0.05, 0.95, 0.9],projection=crs.LambertConformal())
ax.set_adjustable('datalim')
ax.set_extent(lit_domain, crs=crs.LambertConformal())
ax.add_feature(counties_feature, linewidth=0.75)
ax.add_feature(states_feature, linewidth=1.75)
vals=[1, 10, 20, 30, 40, 50, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 1500, 2000, 2500, 3000, 3500, 4000, 4500, 5000, 6000, 7000, 8000, 9000, 10000, 15000, 20000, 25000, 30000, 35000, 40000, 45000, 50000]
cmap = clr.LinearSegmentedColormap.from_list('Population',[(0, '#1e6eeb'), (0.15, '#e1ffff'),(0.25, '#0fa00f'),(0.255, '#fffaaa'),(0.40, '#e11400'),(0.55, '#78000c'),(0.79, '#e96f58'),(0.80, '#643c32'),(1, '#f0dcd2')], N=256)
bounds=[1, 10, 20, 30, 40, 50, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 1500, 2000, 2500, 3000, 3500, 4000, 4500, 5000, 6000, 7000, 8000, 9000, 10000, 15000, 20000, 25000, 30000, 35000, 40000, 45000, 50000]
norm = mpl.colors.BoundaryNorm(bounds, cmap.N)
im = plt.contourf(xtif, ytif, tiffimage, levels=vals, cmap=cmap, norm=norm, alpha=0.85, transform=crs.LambertConformal())
plt.colorbar(im, pad=0.01, boundaries=[1, 10, 20, 30, 40, 50, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 1500, 2000, 2500, 3000, 3500, 4000, 4500, 5000, 6000, 7000, 8000, 9000, 10000, 15000, 20000, 25000, 30000, 35000, 40000, 45000, 50000], ticks=[1, 10, 20, 30, 40, 50, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 1500, 2000, 2500, 3000, 3500, 4000, 4500, 5000, 6000, 7000, 8000, 9000, 10000, 15000, 20000, 25000, 30000, 35000, 40000, 45000, 50000], spacing='uniform', aspect=40)
My Initial Approach
My initial approach was to execute the following lines of code (not shown above since it did not produce the desired result)
states_fill = ShapelyFeature(Reader(states).geometries(), crs.LambertConformal(), facecolor='#C7C6C1', edgecolor='none')
ax.add_feature(states_fill, linewidth=1.75)
im = plt.contourf(xtif, ytif, tiffimage, levels=vals, cmap=cmap, norm=norm, alpha=0.85, transform=crs.LambertConformal())
The addition and placement of ax.add_feature(states_fill, linewidth=1.75) plotted a gray filled shapefile of all U.S states, which was expected and desired. However, it placed this shapefile over the contoured plot, even though I had ordered this correctly before the execution of plt.contourf.
What is the best approach to place a filled contour underneath the contoured plot? Any help is greatly appreciated!
