My example of Python code for mapping variables overlay with geographic information. Elements in the map include ocean, coastline,borders,lakes,rivers,province borders, contour, contourf, pcolormesh, scatter, labels, province (or state) names. Also includes information about how to set cmap, what range of the data (minimum, maximum) to show, how many levels to display the data. How to set title or subtitles of the map.
Import packages
# Ignore warnings
import sys
import warnings
if not sys.warnoptions:
warnings.simplefilter('ignore')
#Process data
import numpy as np
import xarray as xr
#Display data
import cartopy
import cartopy.crs as ccrs
from cartopy.mpl.geoaxes import GeoAxes
from cartopy.vector_transform import vector_scalar_to_grid
from matplotlib.axes import Axes
import cartopy.feature as cfeature
from cartopy import config
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from mpl_toolkits.axes_grid1 import AxesGrid
import matplotlib.image as image
import matplotlib.pyplot as plt
#Writing data files
import pandas as pd
Read data
im = image.imread('figures/LampsYorkuLogo.gif')
fn1 = 'data/ncfiles/reanalysis-era5-land-monthly-means.nc'
ds=xr.open_dataset(fn1)
t1 = ds["tp"][0:360]*1000*365.25 # monthly averages data for 1986-2005
t1m=t1.mean(dim='time')
Plot map
provinces= ['Ontario', 'Quebec','Manitoba','Wisconsin','Vermont','Nebraska','New York','Kansas','Illinois',
'Delaware','Connecticut','Indiana','Missouri','Michigan','New Jersey','Kentucky','Minnesota',
'Ohio','Iowa','Pennsylvania','Maryland','Virginia','West Virginia','North Dakota','South Dakota','Nunavut']
provinces1=['ON','QC','MB','WI','VT','NE','NY','KS','IL','DE','CT','IN','MO','MI','NJ','KY','MN','OH','IA','PA','MD','VA','WV','ND','SD','NU']
Latitudes= [50,53,56.4,44.5,44,41.5,43,38.5,40,39,41.6,40.3,38.6,44.2,39.8,37.8,46.4,
40.4,42,41.2,39,38,39,47.7,44.5,62]
Longitudes=[-85,-76,-98.7,-89.5,-72,-100,-75,-98,-89,-73.5,-72.7,-86,-92.6,-84.5,-74.9,-84.3,-94.6,
-83,-93.6,-77.2,-76.6,-78,-80.5,-99,-99,-98]
lats = t1.coords['latitude'][:]
lons = t1.coords['longitude'][:]
ct_x=[-79.617,-75.717,-86.917,-93.783]
ct_y=[43.667,45.383,49.767,51.067]
ct_n=['Toronto','Ottawa','Geraldton','Red Lake']
projection = ccrs.PlateCarree()
provinc_bodr = cartopy.feature.NaturalEarthFeature(category='cultural',
name='admin_1_states_provinces_lines', scale='50m', facecolor='none', edgecolor='k')
axes_class = (GeoAxes,
dict(map_projection=projection))
# lons, lats = np.meshgrid(lons, lats)
title_text=["season=DJF(m/s)", "season=MAM(m/s)", "season=JJA(m/s)", "season=SON(m/s)"]
fig = plt.figure(figsize=(15,15))
axgr = AxesGrid(fig, 111, axes_class=axes_class,
nrows_ncols=(1, 1),
axes_pad=0.6,
cbar_location='right',
cbar_mode='single',
cbar_pad=0.2,
cbar_size='3%',
label_mode='') # note the empty labe
for i, ax in enumerate(axgr):
#*************************************************************
#add ocean, coastline,borders,lakes,rivers,provinc_bodr
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.COASTLINE,linewidth=0.3)
ax.add_feature(cfeature.BORDERS, linestyle='-', alpha=0.5)
ax.add_feature(cfeature.LAKES, alpha=0.5)
ax.add_feature(cfeature.RIVERS)
ax.add_feature(provinc_bodr, linestyle='-', linewidth=1, edgecolor="k", zorder=10, alpha=0.5)
ax.set_title('Annual Total Precipitation:1981-2010 (mm)')
ax.set_xticks(np.linspace(-100, -70, 5), crs=projection)
ax.set_yticks(np.linspace(35, 65, 5), crs=projection)
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
ax.imshow(im, aspect='auto',extent=(-75,-70,35,37),zorder=-1)
ax.text(-75,37.5,"yorku.ca/ocdp",fontsize=18)
ax.text(-85,60,"Hudson Buy", fontsize=20,color='b',alpha=0.5)
ax.annotate(ct_n[0], (ct_x[0], ct_y[0]),va="top", ha="center",fontsize=14)
ax.annotate(ct_n[1], (ct_x[1], ct_y[1]),va="top", ha="center",fontsize=14)
ax.annotate(ct_n[2], (ct_x[2], ct_y[2]),va="top", ha="center",fontsize=14)
ax.annotate(ct_n[3], (ct_x[3], ct_y[3]),va="top", ha="center",fontsize=14)
for iprov in range(26):
if Longitudes[iprov]<-80 and Latitudes[iprov]<47:
ax.text(Longitudes[iprov],Latitudes[iprov],provinces1[iprov],fontsize=18)
p = ax.contourf(lons, lats, t1m,50,
transform=projection,
cmap='BrBG')
# X,Y = np.meshgrid(lons,lats)
# p = ax.pcolormesh(X, Y, t1m,
# vmin=-10,
# vmax=10,
# transform=projection,
# cmap='Reds')
ax.scatter(ct_x, ct_y,transform=ccrs.PlateCarree(),color='r')
l = ax.contour(lons, lats, t1m,25,colors=['black'], linewidth=0.3,
transform=ccrs.PlateCarree())
ax.clabel(
l, # Typically best results when labelling line contours.
colors=['black'],
manual=False, # Automatic placement vs manual placement.
inline=True, # Cut the line where the label will be placed.
fmt=' {:.0f} '.format, # Labes as integers, with some extra space.
)
axgr.cbar_axes[0].colorbar(p)
plt.tight_layout()
plt.show()
fig.savefig('Annual_Total_Precipitation_1981to2010_30Years.png',dpi=150)