02_Xarray: Plotting

02_Xarray: Plotting#

Overview#

  1. Work with an Xarray Dataset

  2. Select a variable from the dataset

  3. Create a contour plot of gridded ERA5 reanalysis data

Imports#

import xarray as xr
import pandas as pd
import numpy as np
from metpy.units import units
import metpy.calc as mpcalc
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt

Work with an Xarray Dataset#

An Xarray Dataset consists of one or more DataArrays. Let’s continue to work with the ERA5 NetCDF file from the previous notebook.

ds = xr.open_dataset('/spare11/atm350/data/199303_era5.nc')

Examine the SLP dataset

ds

<xarray.Dataset> Size: 7GB
Dimensions:                  (latitude: 721, level: 13, longitude: 1440,
                              time: 124)
Coordinates:
  * latitude                 (latitude) float32 3kB 90.0 89.75 ... -89.75 -90.0
  * level                    (level) int64 104B 50 100 150 200 ... 850 925 1000
  * longitude                (longitude) float32 6kB 0.0 0.25 ... 359.5 359.8
  * time                     (time) datetime64[ns] 992B 1993-03-01 ... 1993-0...
Data variables:
    geopotential             (time, level, latitude, longitude) float32 7GB ...
    mean_sea_level_pressure  (time, latitude, longitude) float32 515MB ...

One attribute of a dataset is its size. We can access that via its nbytes attribute.

print (f'Size of dataset: {ds.nbytes / 1e9} GB')

Size of dataset: 7.2095483 GB

Select a variable from the dataset and assign it to a DataArray object

slp = ds['mean_sea_level_pressure']
slp

<xarray.DataArray 'mean_sea_level_pressure' (time: 124, latitude: 721,
                                             longitude: 1440)> Size: 515MB
[128741760 values with dtype=float32]
Coordinates:
  * latitude   (latitude) float32 3kB 90.0 89.75 89.5 ... -89.5 -89.75 -90.0
  * longitude  (longitude) float32 6kB 0.0 0.25 0.5 0.75 ... 359.2 359.5 359.8
  * time       (time) datetime64[ns] 992B 1993-03-01 ... 1993-03-31T18:00:00
Attributes:
    long_name:      Mean sea level pressure
    short_name:     msl
    standard_name:  air_pressure_at_mean_sea_level
    units:          Pa

Create a contour plot of gridded data#

We got a quick plot in the previous notebook. Let’s now customize the map, and focus on a particular region. We will use Cartopy and Matplotlib to plot our SLP contours. Matplotlib has two contour methods:

  1. Line contours: ax.contour

  2. Filled contours: ax.contourf

We will draw contour lines in hPa, with a contour interval of 4.

As we’ve done before, let’s first define some variables relevant to Cartopy.#

The dataset’s x- and y- coordinates are longitude and latitude … thus, its projection is Plate Carree. However, we wish to display the data on a regional map that uses Lambert Conformal. So we will define two Cartopy projection objects here. We will need to transform the data from its native projection to the map’s projection when it comes time to draw the contours.

latN = 55
latS = 20
lonW = -90
lonE = -60

cLon = (lonW + lonE ) / 2
cLat = (latS + latN ) / 2

proj_data = ccrs.PlateCarree() # Our data is lat-lon; thus its native projection is Plate Carree.
proj_map = ccrs.LambertConformal(central_longitude=cLon, central_latitude=cLat) # Map projection
res = '50m'

Now define the range of our contour values and a contour interval. 4 hPa is standard for a sea-level pressure map on a synoptic-scale region.#

minVal = 900
maxVal = 1076
cint = 4
cintervals = np.arange(minVal, maxVal, cint)
cintervals

array([ 900,  904,  908,  912,  916,  920,  924,  928,  932,  936,  940,
        944,  948,  952,  956,  960,  964,  968,  972,  976,  980,  984,
        988,  992,  996, 1000, 1004, 1008, 1012, 1016, 1020, 1024, 1028,
       1032, 1036, 1040, 1044, 1048, 1052, 1056, 1060, 1064, 1068, 1072])

Matplotlib’s contour methods require three arrays to be passed to them … x- and y- arrays (longitude and latitude in our case), and a 2-d array (corresponding to x- and y-) of our data variable. So we need to extract the latitude and longitude coordinate variables from our DataArray. We’ll also extract the third coordinate value, time.#

lats = ds.latitude
lons = ds.longitude
times = ds.time

Let’s use Xarray’s sel method to specify one time from the DataArray.

slp_singleTime = slp.sel(time='1993-03-14-18:00:00')

slp_singleTime

<xarray.DataArray 'mean_sea_level_pressure' (latitude: 721, longitude: 1440)> Size: 4MB
[1038240 values with dtype=float32]
Coordinates:
  * latitude   (latitude) float32 3kB 90.0 89.75 89.5 ... -89.5 -89.75 -90.0
  * longitude  (longitude) float32 6kB 0.0 0.25 0.5 0.75 ... 359.2 359.5 359.8
    time       datetime64[ns] 8B 1993-03-14T18:00:00
Attributes:
    long_name:      Mean sea level pressure
    short_name:     msl
    standard_name:  air_pressure_at_mean_sea_level
    units:          Pa

Note the units are in Pascals. We will exploit MetPy’s unit conversion library soon, but for now let’s just divide the array by 100.

slp_singleTimeHPA = slp_singleTime / 100

We’re set to make our map. After creating our figure, setting the bounds of our map domain, and adding cartographic features, we will plot the contours. We’ll assign the output of the contour method to an object, which will then be used to label the contour lines.#

fig = plt.figure(figsize=(18,12))
ax = plt.subplot(1,1,1,projection=proj_map)
ax.set_extent ([lonW,lonE,latS,latN])
ax.add_feature(cfeature.COASTLINE.with_scale(res))
ax.add_feature(cfeature.STATES.with_scale(res))
CL = ax.contour(lons,lats,slp_singleTimeHPA,cintervals,transform=proj_data,linewidths=1.25,colors=['green','red','red'])
ax.clabel(CL, inline_spacing=0.2, fontsize=11, fmt='%.0f');
../../_images/f31ada6788d626348d844d3f6c2ff99b4e5531d5ea3f795596408ba78219df11.png 
fig
../../_images/f31ada6788d626348d844d3f6c2ff99b4e5531d5ea3f795596408ba78219df11.png

There we have it! There are definitely some things we can improve about this plot (such as the lack of an informative title, and the fact that our contour labels seem not to be appearing on every contour line), but we’ll get to that in the next lesson!#