Plotting Wind from ERA5 Reanalysis

This article demonstrates how one can plot wind maps from ERA5 renalaysis using python. ERA5 is a 5th generation reanalysis dataset produced by the European Centre for Medium Range Weather Forecasts (ECMWF). The data can be downloaded from the Copernicus Climate Data Store (https://cds.climate.copernicus.eu/#!/home). Here, July 2020 winds over North America are plotted.

ERA5 wind speed and direction over North America for selected days within the hurricane month, July 2020.

The code for making the plot is written below;

#Import libraries
import xarray as xr
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeature
#Read-in file
file = xr.open_dataset('C:/python_work/phd/atm_science/wind_north_america.nc')
#Select the variables
wspd = np.mean(np.sqrt(np.add(np.square(file.u),np.square(file.v))),axis=0)
u = np.mean(file.u,axis=0)
v = np.mean(file.v,axis=0)
#setting up the quiver arguments (these display the wind vectors on the plot)
xx = wspd.longitude.values
yy = wspd.latitude.values
X,Y =np.meshgrid(xx, yy)
U = u.data
V = v.data
# Set figure parameters
fig = plt.figure(figsize=(12, 8), dpi=500)
mpl.rcParams['font.family'] = 'Times New Roman'
plt.rcParams['font.size'] = 18
plt.rcParams['axes.linewidth'] = 1
# Plot the map
ax = plt.axes(projection=ccrs.PlateCarree())
plot = plt.pcolormesh(wspd['longitude'], wspd['latitude'], wspd, cmap='Blues',
                          vmin=0, vmax=16.3)
plt.xlim([-147.4, -43.75])
plt.ylim([13.64, 71.7])
plt.quiver(X[::15, ::15], Y[::15, ::15], U[::15, ::15], V[::15, ::15], 
           transform=ccrs.PlateCarree(), color='k', scale=200, width=0.004)
# Add geographical features to the map
states = cfeature.NaturalEarthFeature('cultural', scale="50m",
                             facecolor="none",
                             name="admin_1_states_provinces_lines")
ax.add_feature(states, linewidth=0.5, edgecolor="black")
ax.coastlines(resolution='10m', color='black', linewidth=0.9)
ax.add_feature(cartopy.feature.BORDERS, linewidth=0.9)
colorbar_axes = plt.gcf().add_axes([0.2, 0.125, 0.6, 0.03])
cb = plt.colorbar(plot, colorbar_axes, 
                  orientation='horizontal', 
                  label='m/s')
plt.show()