<jupyter_start><jupyter_text>Xarray 5: Time Series  Overview1. Work with an ERA-5 `Dataset` hosted on NCAR's Research Data Archive2. Subset the Dataset along its dimensions3. Perform diagnostic calculations and unit conversions4. Perform Split-Apply-Combine5. Create and refine a time series plot of minimum SLP and maximum windspeed over the subsetted region  Prerequisites| Concepts | Importance | Notes || --- | --- | --- || Xarray Lessons 1-4| Necessary | |* **Time to learn**: 30 minutes***  Imports<jupyter_code>import xarray as xr
import pandas as pd
import numpy as np
from datetime import datetime as dt
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
from matplotlib.dates import DateFormatter, AutoDateLocator,HourLocator,DayLocator,MonthLocator
from matplotlib.ticker import (MultipleLocator, AutoMinorLocator)
import requests<jupyter_output><empty_output><jupyter_text>Work with an ERA-5 Dataset hosted on NCAR's [Remote Data Archive](https://rda.ucar.edu/datasets/ds633.0/description) Create OpenDAP URLs pointing to the variables as they are stored in RDA's THREDDS server Note: Single-level fields are stored in one monthly file. Fields on pressure levels are stored in daily files. In this case, we will retrieve SLP and 10-meter winds, so we will be using the sfc directory for all three variables.<jupyter_code>ds_urlu = 'https://rda.ucar.edu/thredds/dodsC/files/g/ds633.0/e5.oper.an.sfc/201210/e5.oper.an.sfc.128_165_10u.ll025sc.2012100100_2012103123.nc'
ds_urlv = 'https://rda.ucar.edu/thredds/dodsC/files/g/ds633.0/e5.oper.an.sfc/201210/e5.oper.an.sfc.128_166_10v.ll025sc.2012100100_2012103123.nc'
ds_urlp = 'https://rda.ucar.edu/thredds/dodsC/files/g/ds633.0/e5.oper.an.sfc/201210/e5.oper.an.sfc.128_151_msl.ll025sc.2012100100_2012103123.nc'<jupyter_output><empty_output><jupyter_text>Use your RDA credentials to access these individual NetCDF files. You should have previously registered with RDA, and then saved your user id and password in your home directory as .rdarc , with permissions such that only you have read/write access.<jupyter_code>#Retrieve login credintial for RDA.
from pathlib import Path
HOME = str(Path.home())
credFile = open(HOME+'/.rdarc','r')
userId, pw = credFile.read().split()<jupyter_output><empty_output><jupyter_text>Connect to the RDA THREDDS server and point to the files of interest.<jupyter_code>session = requests.Session()
session.auth = (userId, pw)
storep = xr.backends.PydapDataStore.open(ds_urlp, session=session)
storeu = xr.backends.PydapDataStore.open(ds_urlu, session=session)
storev = xr.backends.PydapDataStore.open(ds_urlv, session=session)<jupyter_output><empty_output><jupyter_text>If you get an error after the cell, check to be sure you have registered @ RDA and you have entered and saved your password file as directed. Now open the three files as separate Datasets<jupyter_code>dsSLP = xr.open_dataset(storep)
dsU10 = xr.open_dataset(storeu)
dsV10 = xr.open_dataset(storev)<jupyter_output><empty_output><jupyter_text>Examine one of the `Dataset`s.<jupyter_code>dsU10<jupyter_output><empty_output><jupyter_text>Create `DataArray` objects for each of the three variables in the `Dataset`s.<jupyter_code>slp = dsSLP['MSL']
u10 = dsU10['VAR_10U']
v10 = dsV10['VAR_10V']<jupyter_output><empty_output><jupyter_text>Subset the Dataset along its dimensions. As we did before, we will perform temporal and areal subsetting.  Take advantage of Pandas' date/time methods. In this case, we will use Pandas [date_range](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.date_range.html). Let's define a function whose output we will use to label our plotted time series.<jupyter_code>timeRange = pd.date_range(start="2012-10-27 00:00",end="2012-10-31 00:00",freq='1H')
timeRange
def areaBox(lonW,lonE,latS,latN):
    '''
    Returns a string containing the bounds of the selected region
    '''
    outputStr = ''
    for val in [lonW, lonE]:
        if (val < 0):
            val *= -1
            valStr = str(val) + 'W '
        else:
            valStr = str(val) + 'E '
        outputStr = outputStr + valStr
    for val in [latS, latN]:
        if (val< 0): 
            val *= -1
            valStr = str(val) + 'S '
        else:
            valStr = str(val) + 'N '
        outputStr = outputStr + valStr
    
    return outputStr
# Areal extent
lonW = -100
lonE = -60
latS = 20
latN = 50
domainStr = areaBox(lonW,lonE,latS,latN)
latRange = np.arange(latS-5,latN+5,.25) # expand the data range a bit beyond the plot range
lonRange = np.arange((lonW-5+360),(lonE+5+360),.25) # Need to match longitude values to those of the coordinate variable

# Vertical level specificaton
vlevel = 10
levelStr = str(vlevel)

# Data variable selection
SLP = dsSLP['MSL'].sel(time=timeRange,latitude=latRange,longitude=lonRange)
U10 = dsU10['VAR_10U'].sel(time=timeRange,latitude=latRange,longitude=lonRange)
V10 = dsV10['VAR_10V'].sel(time=timeRange,latitude=latRange,longitude=lonRange)<jupyter_output><empty_output><jupyter_text>Examine a couple of the DataArrays following our subsetting<jupyter_code>SLP
U10<jupyter_output><empty_output><jupyter_text>Xarray has [Split-Apply-Combine](http://xarray.pydata.org/en/stable/groupby.html) functionality similar to Pandas. In this next cell, we are grouping by date/time and then calculating the minimum value of SLP over the other two dimensions.<jupyter_code>SLP.min(dim=['latitude','longitude'])<jupyter_output><empty_output><jupyter_text>Call Xarray's built-in `plot` method to get a quick look at the time series. Note that we use the Python convention where `_` is an  object corresponding to the output of the previous line of code.<jupyter_code>_.plot()<jupyter_output><empty_output><jupyter_text>Perform unit conversions Convert meters per second to knots, and Pascals to hectoPascals. Both are straightforward with MetPy's `convert_units` method.<jupyter_code>SLP = SLP.metpy.convert_units('hPa')
U10 = U10.metpy.convert_units('kts')<jupyter_output><empty_output><jupyter_text>V10 = V10.metpy.convert_units('kts') Calculate wind speed from the u- and v- components using MetPy's diagnostic function, `wind_speed`.<jupyter_code>WSPD = mpcalc.wind_speed(U10,V10)<jupyter_output><empty_output><jupyter_text>Perform Split-Apply-Combine Compute the minimum SLP, and then the maximum wind speed, for all latitude and longitude points ... for each time.<jupyter_code>SLPmin = SLP.min(dim=['latitude','longitude'])
print(SLPmin)
WSPDmax = WSPD.max(dim=['latitude','longitude'])<jupyter_output><empty_output><jupyter_text>We'll take a deeper dive into Xarray's *Split-Apply-Combine* implementation in the next notebook. Here, we've effectively *reduced* the data array down to one dimension (time)!  Create and refine a time series plot of minimum SLP and maximum windspeed over the subsetted region. Create a meaningful title string.<jupyter_code>tl1 = "ERA-5 maximum " + levelStr + " m windspeed (kts) and minimum SLP (hPa)"
tl2 = str('Valid within: '+ domainStr)
title_line = (tl1 + '\n' + tl2 + '\n')
plt.style.use('seaborn')<jupyter_output><empty_output><jupyter_text>Plot two y-axes; one for SLP, the other for windspeed<jupyter_code>fig = plt.figure(figsize=(18,12))
ax1 = plt.subplot(1,1,1)
color = 'tab:red'
ax1.plot(SLP.time,SLPmin,color=color,label = 'Min SLP')

ax1.set_title (title_line,fontsize=14)

ax2 = ax1.twinx()  # instantiate a second axes that shares the same x-axis
color = 'tab:blue'
ax2.plot(WSPD.time,WSPDmax,color=color,label = 'Max Wind');<jupyter_output><empty_output><jupyter_text>What's the meaning of `tab:` in our `color` definition? It's one of many [color table specifications](https://matplotlib.org/stable/gallery/color/colormap_reference.html) available in Matplotlib; in this case, *tab* refers to a [color table](https://matplotlib.org/stable/tutorials/colors/colors.html) derived from a commercial dataviz package, [Tableau](https://www.tableau.com).  That's not a bad first plot, but much can be done to improve it. Here we use several Matplotlib `Axes` methods to customize the plot.  Often the best way to learn about these methods (and others) is to browse the relevant documentation at https://matplotlib.org and play around with them. You may be the type who always wants to keep making it look better, but avoid the temptation to keep trying for perfection.<jupyter_code>fig = plt.figure(figsize=(18,12),dpi=150)
ax1 = plt.subplot(1,1,1)

# First Axes: SLP
color = 'tab:red'
ax1.plot(SLP.time,SLPmin,color=color,label = 'Min SLP')
ax1.grid(color='black',linewidth=0.5)

ax1.xaxis.set_major_locator(HourLocator(interval=6))
ax1.xaxis.set_minor_locator(HourLocator(interval=1))
dateFmt = DateFormatter('%H')
ax1.xaxis.set_major_formatter(dateFmt)
ax1.set_xlabel('Date/Time (UTC)')
ax1.set_xlim(SLP.time[0],SLP.time[-1])

ax1.set_ylabel('SLP (hPa)',color=color, fontsize=14)
ax1.yaxis.set_minor_locator(MultipleLocator(5))
ax1.tick_params(which='minor', length=6)
ax1.tick_params(which='major',axis='y', labelcolor=color,direction='out', length=6, width=2, colors=color,
               grid_color=color, grid_alpha=0.5)
ax1.tick_params(which='minor', axis = 'y', length=4, color=color)

ax1.set_title (title_line,fontsize=14)

# Second Axes: Windspeed
ax2 = ax1.twinx()  # instantiate a second axes that shares the same x-axis
color = 'tab:blue'
ax2.plot(WSPD.time,WSPDmax,color='tab:blue',label = 'Max Wind')
# Do not draw grid lines (makes plot more readable)
ax2.grid(b=None)

ax2.set_ylabel('Windspeed (kts)', color=color, fontsize=14)  # we already handled the x-label with ax1
ax2.tick_params(axis='y', labelcolor=color,direction='out', length=6, width=2, colors=color,
               grid_color=color, grid_alpha=0.5) # grid-related specs will be ignored since we are not plotting grid lines for this axes

fig.legend(loc='upper right',frameon=True,bbox_to_anchor=(0.9,1.1), bbox_transform=ax1.transAxes,fontsize=16,shadow=True,edgecolor='purple')
fig.autofmt_xdate(rotation=45);<jupyter_output><empty_output><jupyter_text>Save our figure to our current directory.<jupyter_code>fig.savefig('ERA5_Time_Series.png')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11416552533750733, 'fraction_numerical': 0.03198982586577969, 'mean_word_length': 5.604005167958657, 'pattern_counts': {'":': 0, '<': 75, '<?xml version=': 0, '>': 73, 'https://': 9, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '580542', 'license': 'CC0-1.0', 'n_tokens_mistral': 3738, 'n_tokens_neox': 3512, 'n_words': 927, 'path': '05_Xarray_TimeSeries.ipynb', 'repo_name': 'ktyle/binder_xarray', 'size': '10222'}