ERA5 Data Analysis¶
This notebook explores some of the variables available in the ERA5 dataset for use in the ATM-505 project.
Required packages: xarray (http://xarray.pydata.org/en/stable/), Numpy, Datetime, Cartopy (https://scitools.org.uk/cartopy/docs/latest/), Matplotlib
Some useful links that go beyond this notebook: ERA5 documentation (https://confluence.ecmwf.int/display/CKB/ERA5%3A+data+documentation), MetPy package for computing various meteorological variables using gridded data (https://unidata.github.io/MetPy/latest/index.html)
Computing distance [meters] from gridpoint lat/lon: MetPy (https://unidata.github.io/MetPy/latest/api/generated/metpy.calc.lat_lon_grid_deltas.html), Pyproj (https://janakiev.com/blog/gps-points-distance-python/)
import xarray as xr
import numpy as np
from datetime import datetime
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import matplotlib.patches as mpatches
%matplotlib inline
Plot 500-hPa Geopotential Height¶
The first step is to load the ERA5 hourly data on pressure levels dataset for your case. In this example we'll look at the 05-08 February 2020 period.
This is where you can learn more about the ERA5 pressure level dataset, and even download data on your own:
https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-pressure-levels?tab=overview
How to cite: Copernicus Climate Change Service (C3S) (2017): ERA5: Fifth generation of ECMWF atmospheric reanalyses of the global climate. Copernicus Climate Change Service Climate Data Store (CDS), date of access. https://cds.climate.copernicus.eu/cdsapp#!/home
filename = '/home/disk/meso-home/jfinlon/ERA5_Data/20200205_20200208/era5.data_on_pressure_levels.2020020500_2020020821.nc'
ds = xr.open_dataset(filename) # This loads the NetCDF file as an xarray object
Let's take a look at how the data are structured in this file.
ds
The domain for this project ranges in latitude from 20N to 65N and in longitude from 170W to 50W. Here's how it looks on a map.
proj = ccrs.LambertConformal(central_longitude=-110.0, central_latitude=35.0)
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection=proj)
ax.stock_img()
ax.add_patch(mpatches.Rectangle(xy=[-170, 20], width=120, height=45, facecolor='red', alpha=0.2, transform=ccrs.PlateCarree()))
ax.set_extent([-180, -30, 10, 70], crs=ccrs.PlateCarree())
ax.gridlines()
ax.coastlines()
plt.show()
Read the time, pressure level, latitude, longitude, and geopotential variables.
time = ds['time'].values # numpy datetime64 object
lvl = ds['level'].values # units of hPa
lat = ds['latitude'].values # units of degrees (postitive is northern hemisphere)
lon = ds['longitude'].values # units of degrees (negative is western hemisphere)
geopot = ds['z'].values # units of m**2 / s**2
# Compute the geopotential height [m] from the geopotential variable
geopot_hght = geopot / 9.81
print(time.shape, lvl.shape, lat.shape, lon.shape, geopot_hght.shape)
The geopotential height variable needs to be trimmed if we're only interested in one time period and one pressure level...
# Obtain the time index that corresponds to 1200 UTC 07 February
time_ind = np.where(time==np.datetime64('2020-02-07T12:00:00.00'))[0][0]
# Obtain the pressure level index that corresponds to 500 hPa
lvl_ind = np.where(lvl==500.)[0][0]
geopot_hght_trimmed = geopot_hght[time_ind, lvl_ind, :, :]
To make a filled contour plot we need to transform the latitude and longitude arrays into bins.
[X, Y] = np.meshgrid(lon, lat)
print(X.shape, Y.shape)
Finally, time to plot the geopotential height! Sections of this cell block are optional and for aesthetic purposes only.
proj = ccrs.Mercator() # You can use other projections, just make sure that PlateCarree is used in any transform arguments below.
fig = plt.figure(figsize=(16, 12))
ax = fig.add_subplot(111, projection=proj)
plt1 = ax.contourf(X, Y, geopot_hght_trimmed/10., levels=np.linspace(480, 600, 41), cmap='nipy_spectral', transform=ccrs.PlateCarree())
plt2 = ax.contour(X, Y, geopot_hght_trimmed/10., levels=np.linspace(480, 600, 21), colors='black', linewidths=2, transform=ccrs.PlateCarree())
ax.clabel(plt2, np.linspace(480, 600, 11), inline=True, fmt='%d', fontsize=14)
ax.set_title('500-hPa Geopotential Height [dm]\n1200 UTC 07 Feb 2020', fontsize=16)
# Add some additional borders (optional)
states = cfeature.NaturalEarthFeature(category='cultural', name='admin_1_states_provinces_lines', scale='50m', facecolor='none')
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.LAKES, facecolor='none', edgecolor='black')
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS)
ax.add_feature(states, edgecolor='black')
# Format the gridlines (optional)
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=2, color='gray', alpha=0.5, linestyle='--')
gl.xlabels_top = False
gl.ylabels_right = False
gl.xlines = True; gl.xlocator = mticker.FixedLocator(np.linspace(-170, -50, 13)); gl.xformatter = LONGITUDE_FORMATTER; gl.xlabel_style = {'size':16, 'color':'black'}
gl.ylines = True; gl.ylocator = mticker.FixedLocator(np.linspace(20, 65, 10)); gl.yformatter = LATITUDE_FORMATTER; gl.ylabel_style = {'size':16, 'color':'black'}
# Plot the colorbar
cbar_ax = fig.add_axes([0, 0, 0.1, 0.1]) # Dummy values prior to finetuning the cbar position
pos = ax.get_position() # Get the axes position
cbar_ax.set_position([pos.x0 + pos.width + 0.01, pos.y0, 0.04, pos.height])
cbar = plt.colorbar(plt1, cax=cbar_ax)
cbar_ax.tick_params(labelsize=16)
plt.show()
Plot Pressure on Dynamic Tropopause Surface¶
The first step is to load the ERA5 hourly data on potential vorticity levels dataset for your case. In this example we'll look at the 05-08 February 2020 period.
As a rule of thumb, the dynamic tropopause is often considered to be found where potental vorticity values are around 2 PVU in the mid-latitudes.
filename = '/home/disk/meso-home/jfinlon/ERA5_Data/20200205_20200208/era5.data_on_pv_levels.2020020500_2020020821.nc'
ds2 = xr.open_dataset(filename) # This loads the NetCDF file as an xarray object
ds2
Read the time, latitude, longitude, and pressure variables.
time = ds2['time'].values # numpy datetime64 object
lat = ds2['latitude'].values # units of degrees (postitive is northern hemisphere)
lon = ds2['longitude'].values # units of degrees (negative is western hemisphere)
pres = ds2['pres'].values / 100. # units of hPa
print(time.shape, lat.shape, lon.shape, pres.shape)
The pressure variable needs to be trimmed if we're only interested in one time period...
# Obtain the time index that corresponds to 1200 UTC 07 February
time_ind = np.where(time==np.datetime64('2020-02-07T12:00:00.00'))[0][0]
pres_trimmed = pres[time_ind, :, :]
Finally, time to plot the dynamic tropopause pressure! Sections of this cell block are optional and for aesthetic purposes only.
[X, Y] = np.meshgrid(lon, lat)
proj = ccrs.Mercator() # You can use other projections, just make sure that PlateCarree is used in any transform arguments below.
fig = plt.figure(figsize=(16, 12))
ax = fig.add_subplot(111, projection=proj)
plt1 = ax.contourf(X, Y, pres_trimmed, levels=np.linspace(40, 720, 35), cmap='nipy_spectral_r', transform=ccrs.PlateCarree())
ax.set_title('Dynamic Tropopause Pressure [hPa]\n1200 UTC 07 Feb 2020', fontsize=16)
# Add some additional borders (optional)
states = cfeature.NaturalEarthFeature(category='cultural', name='admin_1_states_provinces_lines', scale='50m', facecolor='none')
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.LAKES, facecolor='none', edgecolor='black')
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS)
ax.add_feature(states, edgecolor='black')
# Format the gridlines (optional)
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=2, color='gray', alpha=0.5, linestyle='--')
gl.xlabels_top = False
gl.ylabels_right = False
gl.xlines = True; gl.xlocator = mticker.FixedLocator(np.linspace(-170, -50, 13)); gl.xformatter = LONGITUDE_FORMATTER; gl.xlabel_style = {'size':16, 'color':'black'}
gl.ylines = True; gl.ylocator = mticker.FixedLocator(np.linspace(20, 65, 10)); gl.yformatter = LATITUDE_FORMATTER; gl.ylabel_style = {'size':16, 'color':'black'}
# Plot the colorbar
cbar_ax = fig.add_axes([0, 0, 0.1, 0.1]) # Dummy values prior to finetuning the cbar position
pos = ax.get_position() # Get the axes position
cbar_ax.set_position([pos.x0 + pos.width + 0.01, pos.y0, 0.04, pos.height])
cbar = plt.colorbar(plt1, cax=cbar_ax)
cbar_ax.tick_params(labelsize=16)
plt.show()
