Adding Provinces to Canadian Map in Python Cartopy

Question

Confused grad student here. I'm looking to add the lines for Canadian provinces the same way that the ax.add_feature(feature.STATES) line does, but just for provinces. Does anyone know how to do this? I've attached my code below and the map that it produces.

        fig, ax = plt.subplots(1, 1, figsize=(10,8))
        ax = plt.axes(projection=ccrs.PlateCarree())
        ax.set_global()
        ax.set_extent([-140,-60,40,75], crs=ccrs.PlateCarree())
        ax.coastlines()
        ax.stock_img()
        ax.set_xticks([-150, -140, -130, -120, -110, -100, -90, -80, -70, -60],crs=ccrs.PlateCarree())
        ax.set_yticks([75, 70, 65, 60, 55, 50, 45],crs=ccrs.PlateCarree())
        ax.add_feature(feature.BORDERS)
        ax.add_feature(feature.STATES)
        ax.set_title('Map of Areas with Projected Viable Summer Temperatures in 2020', size = 'xx-large')
        ax.add_feature(feature.LAND)
        ax.add_feature(feature.COASTLINE)

        # plot data
        pc = ax.pcolormesh(lon,lat,tempsummer2020viable,cmap="OrRd")

        cax,kw = mpl.colorbar.make_axes(ax,location='bottom',pad=0.05,shrink=0.7)
        out=fig.colorbar(pc,cax=cax,extend='both',**kw)
        out.set_label('Temperature Within Viable Limit (°C))',size=10)

Image of my map

user_2177413 · Answer

Here is a sample code to plot a map of Canada that has some basic features including the province boundaries, plus a sample thematic layer.

# for Canada
import matplotlib.pyplot as plt
import cartopy
import cartopy.feature as cfeature
import cartopy.crs as ccrs
import numpy as np

extent = [-130, -55, 36.5, 75]
central_lon = np.mean(extent[:2])
central_lat = np.mean(extent[2:])

plt.figure(figsize=(12, 8))

ax = plt.axes(projection=ccrs.AlbersEqualArea(central_lon, central_lat))
ax.set_extent(extent)

# data resolution
resol = '50m'

# country boundaries
country_bodr = cartopy.feature.NaturalEarthFeature(category='cultural', 
    name='admin_0_boundary_lines_land', scale=resol, facecolor='none', edgecolor='k')

# province boundaries
provinc_bodr = cartopy.feature.NaturalEarthFeature(category='cultural', 
    name='admin_1_states_provinces_lines', scale=resol, facecolor='none', edgecolor='k')

# land areas
land = cartopy.feature.NaturalEarthFeature('physical', 'land', \
    scale=resol, edgecolor='k', facecolor=cfeature.COLORS['land'])

# Ocean/seas
ocean = cartopy.feature.NaturalEarthFeature('physical', 'ocean', \
    scale=resol, edgecolor='none', facecolor=cfeature.COLORS['water'])

# Lakes
lakes = cartopy.feature.NaturalEarthFeature('physical', 'lakes', \
    scale=resol, edgecolor='b', facecolor=cfeature.COLORS['water'])

# Rivers
rivers = cartopy.feature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', \
    scale=resol, edgecolor='b', facecolor='none')

# Add all features to the map 
ax.add_feature(land, facecolor='beige', zorder=4)
ax.add_feature(ocean, linewidth=0.2 )
ax.add_feature(lakes, zorder=5)
ax.add_feature(rivers, linewidth=0.5, zorder=6)

ax.add_feature(country_bodr, linestyle='--', linewidth=0.8, edgecolor="k", zorder=10)  #USA/Canada
ax.add_feature(provinc_bodr, linestyle='--', linewidth=0.6, edgecolor="k", zorder=10)

# --- Begin: plot some raster-thematic layer
# simulate a spatial pattern across the map
xlims = (-130, -55)
ylims = (40, 70)
resolution = 0.2
y, x = np.mgrid[slice(ylims[0], ylims[1] + resolution, resolution),
               slice(xlims[0], xlims[1] + resolution, resolution)]
z = -x + np.sin(x)**2 + np.cos(y)
im = ax.pcolormesh(x, y, z, cmap='viridis_r', zorder=7, alpha=0.2, transform=ccrs.PlateCarree())
plt.colorbar(im, ax=ax, shrink=0.5)  #sample colorbar
# --- End: thematic layer plotting.

ax.gridlines(draw_labels=True, lw=1.2, edgecolor="darkblue", zorder=12)

plt.show()

And the sample map:

Adding Provinces to Canadian Map in Python Cartopy

1 Answers