Feature 447 stratosphere qbo (#450)

* Adding phase diagrams for Stratosphere QBO

* Updating documentation

* Adding timeseries plot

* Updating stratosphere plots

docs/Users_Guide/stratosphere_plots.rst

@@ -6,11 +6,11 @@ Description
 ===========
 
 The **stratosphere_plots.py** script contains the plotting portion for
-two Stratosphere use cases, one which creates a ME plot in latitude and pressure
-and another which creates ME and RMSE plots for lead time and pressure
-Two METplus use cases, illustrate how to use this plot.  One runs `zonal mean biases
-<https://metplus.readthedocs.io/en/develop/generated/model_applications/s2s/UserScript_obsERA_obsOnly_WeatherRegime.html#sphx-glr-generated-model-applications-s2s-userscript-fcstgfs-obsera-stratospherebias-py>`_
-and another creates bias and RMSE for `polar cap temperature and polar vortex U <https://metplus.readthedocs.io/en/develop/generated/model_applications/s2s/UserScript_obsERA_obsOnly_WeatherRegime.html#sphx-glr-generated-model-applications-s2s-userscript-fcstgfs-obsera-stratospherepolar-py>`_
+three Stratosphere use cases.  One use case creates a ME plot in latitude and pressure,
+another which creates ME and RMSE plots for lead time and pressure, and a third which 
+creates two phase diagrams and a time series of U for at 50mb and 30mb.
+The three METplus use cases, illustrate how to use these plotting scripts for `zonal mean biases
+<https://metplus.readthedocs.io/en/develop/generated/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereBias.html#sphx-glr-generated-model-applications-s2s-stratosphere-userscript-fcstgfs-obsera-stratospherebias-py>`_ , creating bias and RMSE for `polar cap temperature and polar vortex U <https://metplus.readthedocs.io/en/develop/generated/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratospherePolar.html#sphx-glr-generated-model-applications-s2s-stratosphere-userscript-fcstgfs-obsera-stratospherepolar-py>`_ and creating `phase diagrams and time series for QBO <https://metplus.readthedocs.io/en/develop/generated/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereQBO.html#sphx-glr-generated-model-applications-s2s-stratosphere-userscript-fcstgfs-obsera-stratosphereqbo-py>`_
 
 These files are used by the image comparison test:
 
@@ -32,6 +32,18 @@ These files are used by the image comparison test:
 * **RMSE_2018_02_polar_vortex_U.png**:  Run "plot_polar_rmse" in **stratosphere_plots.py**
   to create this plot.
 
+* **ERA_GFS_QBO_circuits.png**:  Run "plot_qbo_phase_circuits" in **stratosphere_plots.py**
+  to create this plot.
+
+* **ERA5_QBO_PhaseSpace.png**:  Run "plot_qbo_phase_space" in **stratosphere_plots.py**
+  to create this plot.
+
+* **ERA_GFS_timeseries_30mb_u_201710_201802.png**:  Run "plot_u_timeseries" in **stratosphere_plots.py**
+  to create this plot.
+
+* **ERA_GFS_timeseries_50mb_u_201710_201802.png**:  Run "plot_u_timeseries" in **stratosphere_plots.py**
+  to create this plot.
+
 
 Required Packages
 =================
@@ -46,6 +58,8 @@ Required Packages
 
 * numpy
 
+* xarray
+
 * pandas
 
 * cmocean
@@ -60,13 +74,13 @@ Import stratosphere_plots into the script:
 
 .. code-block:: ini
 
-   from stratosphere_plots import plot_zonal_bias,plot_polar_bias,plot_polar_rmse
+   from stratosphere_plots import plot_zonal_bias,plot_polar_bias,plot_polar_rmse, plot_qbo_phase_circuits,plot_qbo_phase_space,plot_u_timeseries
 
 For plot_zonal_bias
 -------------------
 
 In the code, generate the following as numpy
-arrays (except outfile, ptitle, and plevels).
+arrays (except outfile, ptitle, and plevs).
 
 **lats:**  A numpy array of the latitude values under consideration.
 
@@ -82,14 +96,14 @@ file, a **.png** version will be written.
 
 **ptitle:** A string containing the title of the plot.
 
-**plevels:** A list containing integers of the contour levels used in
+**plevs:** A list containing integers of the contour levels used in
 plotting the obs climatology.
 
 For plot_polar_bias
 -------------------
 
 In the code, generate the following as numpy arrays
-(except wrnum, output_plotname, and plevels).
+(except outfile, ptitle, and plevs).
 
 **leads:**  A numpy array containing the forecast lead times.
 
@@ -103,27 +117,83 @@ file, a **.png** version will be written.
 **ptitle:**  A string containing the title of the plot.
 
 **plevs:**  A list containing floats of the contour levels used in
-plotting
+plotting.
 
 For plot_polar_rmse
 -------------------
 
 In the code, generate the following as numpy arrays
-(except wrnum, output_plotname, and plevels).
+(except outfile, ptitle, and plevs).
 
 **leads:**  A numpy array containing the forecast lead times.
 
 **levels:**  A numpy array of the pressure level values under consideration.
 
-**pdata:**  A numpy array containing the bias.
+**pdata:**  A numpy array containing the RMSE.
 
 **outfile:**  The full path and filename of the output plot
 file, a **.png** version will be written.
 
 **ptitle:**  A string containing the title of the plot.
 
 **plevs:**  A list containing floats of the contour levels used in
-plotting
+plotting.
+
+For plot_qbo_phase_circuits
+---------------------------
+
+In the code, generate the following as numpy arrays
+(except inits, periods, and outfile).
+
+**inits:**  A listing of datetimes that are the start date for each plot.
+
+**periods:**  An integer containing the number of days to plot from the inits.
+
+**rean_qbo_pcs:**  An xarray dataarray containing the projected daily 
+zonal winds for the observations.
+
+**rfcst_qbo_pcs:**  An xarray dataarray containing the projected 
+daily zonal winds for the model.
+
+**outfile:**  The full path and filename of the output plot
+file, a **.png** version will be written.
+
+For plot_qbo_phase_space
+------------------------
+
+In the code, generate the following as numpy arrays
+(except ptitle and outfile).
+
+**rean_qbo_pcs:**  An xarray dataarray containing the projected 
+daily zonal winds. 
+
+**eofs:**  An xarray dataarray containing the EOFs.
+
+**ptitle:**  A string containing the title of the plot.
+
+**outfile:**  The full path and filename of the output plot
+file, a **.png** version will be written.
+
+For plot_u_timeseries
+---------------------
+
+In the code, generate the following as numpy arrays
+(except ptitle and outfile).
+
+**obs_dt:**  A numpy array of datetimes for the observations.
+
+**obs_u:**  A numpy array containing U wind values for 
+the observations.
+
+**fcst_dt:**  A numpy array of datetimes for the forecasts.
+
+**fcst_u:**  A numpy array containing U wind values for 
+the forecasts.
+
+**ptitle:**  A string containing the title of the plot.
+
+**outfile:**  The full path and filename of the output plot
+file, a **.png** version will be written.
 
 Invoke the plotting functions:
 
@@ -135,7 +205,13 @@ Invoke the plotting functions:
 
    plot_polar_rmse(leads,levels,pdata,outfile,ptitle,plevs)
 
+   plot_qbo_phase_circuits(inits,periods,rean_qbo_pcs,rfcst_qbo_pcs,outfile)
+
+   plot_qbo_phase_space(rean_qbo_pcs,eofs,ptitle,outfile)
+
+   plot_u_timeseries(obs_dt,obs_u,fcst_dt,fcst_u,ptitle,outfile)
+
 The output will be **.png** version of all requested plots and will
 be located based on what was specified (path and name) in the
-**output_plotname**.
+**outfile**.
 

metplotpy/contributed/stratosphere_diagnostics/stratosphere_plots.py

@@ -1,5 +1,6 @@
 import xarray as xr  # http://xarray.pydata.org/
 import numpy as np
+import pandas as pd
 import cmocean
 import matplotlib
 from matplotlib import cm
@@ -55,3 +56,143 @@ def plot_polar_contour(leads,levels,pdata,outfile,ptitle,plevs,ctable):
     plt.savefig(outfile)
     plt.close()
 
+
+def plot_qbo_phase_circuits(inits,periods,rean_qbo_pcs,rfcst_qbo_pcs,outfile):
+
+    fig = plt.figure(1)
+
+    for i,init in enumerate(inits):
+
+        dates = pd.date_range(init, periods=periods, freq="1D")
+        #Remove leap days
+        dates = dates[(dates.month != 2) | (dates.day != 29)] 
+        ax = fig.add_subplot(4,7,i+1)
+
+        # plot the QBO phase space time series
+        plt_handles = []
+
+        # green/red dots are the start/end of the fcst period
+        plt.plot(rfcst_qbo_pcs.sel(mode=0, time=dates),
+             rfcst_qbo_pcs.sel(mode=1, time=dates), linewidth=1.5)
+        plt.plot(rfcst_qbo_pcs.sel(mode=0, time=dates[0]),
+             rfcst_qbo_pcs.sel(mode=1, time=dates[0]), marker="o", color="green")
+        plt.plot(rfcst_qbo_pcs.sel(mode=0, time=dates[-1]),
+             rfcst_qbo_pcs.sel(mode=1, time=dates[-1]), marker="o", color="red")
+
+        # plot same thing from reanalysis in black
+        plt.plot(rean_qbo_pcs.sel(mode=0, time=dates),
+             rean_qbo_pcs.sel(mode=1, time=dates), color='black', linewidth=1.5)
+        plt.plot(rean_qbo_pcs.sel(mode=0, time=dates[0]),
+             rean_qbo_pcs.sel(mode=1, time=dates[0]), marker="o", color="green")
+        plt.plot(rean_qbo_pcs.sel(mode=0, time=dates[-1]),
+             rean_qbo_pcs.sel(mode=1, time=dates[-1]), marker="o", color="red")
+        plt.xlim((-2.5,2.5))
+        plt.ylim((-2.5,2.5))
+        plt.gca().set_aspect('equal')
+        plt.title('Start: '+init.strftime('%Y-%m-%d'), fontsize=22)
+
+        # plot the QBO phase lines
+        amp = 2.5
+        ax.plot([-amp,amp],[-amp,amp], linewidth=0.5, linestyle='--', color='k', zorder=0)
+        ax.plot([-amp,amp],[amp,-amp], linewidth=0.5, linestyle='--', color='k', zorder=0)
+        ax.plot([-amp,amp],[0,0], linewidth=0.5, linestyle='--', color='k', zorder=0)
+        ax.plot([0,0],[-amp,amp], linewidth=0.5, linestyle='--', color='k', zorder=0)
+        ax.set_xticks(np.arange(-2.,2.1,1))
+        ax.set_yticks(np.arange(-2.,2.1,1))
+
+
+    fig.set_size_inches(24,16)
+    fig.subplots_adjust(hspace = 0.25, left=0.07, right=0.92, top=0.92, bottom=0.07)
+    plt.savefig(outfile, dpi=150, facecolor='white', bbox_inches="tight")
+
+
+
+def plot_qbo_phase_space(rean_qbo_pcs,eofs,ptitle,outfile):
+
+    from mpl_toolkits.axes_grid1.inset_locator import (inset_axes, InsetPosition, mark_inset)
+
+    fig = plt.figure(4)
+
+    ax = fig.add_subplot(1,1,1)
+    ax.plot(rean_qbo_pcs.sel(mode=0), rean_qbo_pcs.sel(mode=1), alpha=0.75, linewidth=0.6)
+    ax.set_xlim((-2.5,2.5))
+    ax.set_ylim((-2.5,2.5))
+    ax.set_aspect('equal')
+
+    amp = 2.5
+    ax.plot([-amp,amp],[-amp,amp], linewidth=0.65, linestyle='--', color='k', zorder=0)
+    ax.plot([-amp,amp],[amp,-amp], linewidth=0.65, linestyle='--', color='k', zorder=0)
+    ax.plot([-amp,amp],[0,0], linewidth=0.65, linestyle='--', color='k', zorder=0)
+    ax.plot([0,0],[-amp,amp], linewidth=0.65, linestyle='--', color='k', zorder=0)
+    ax.set_xticks(np.arange(-2.,2.1,1))
+    ax.set_yticks(np.arange(-2.,2.1,1))
+
+    ax.set_xlabel('EOF1', fontsize=20)
+    ax.set_ylabel('EOF2', fontsize=20)
+    ax.tick_params(which='major', axis='both', length=7, labelsize=18)
+    ax.grid(True,alpha=0.5)
+
+    ## Plot idealized zonal wind profiles for different combinations of
+    ## EOF values
+    inset_axes = []
+    axi = ax.inset_axes([2,-0.25,0.5,0.5], transform=ax.transData)
+    axi.plot(eofs.sel(mode=0), eofs.pres, linewidth=1.5)
+    inset_axes.append(axi)
+
+    axi = ax.inset_axes([np.sqrt(2)+0.2, np.sqrt(2)+0.2, 0.5, 0.5], transform=ax.transData)
+    axi.plot(eofs.sel(mode=0)+eofs.sel(mode=1), eofs.pres, linewidth=1.5)
+    inset_axes.append(axi)
+
+    axi = ax.inset_axes([-0.25,2, 0.5, 0.5], transform=ax.transData)
+    axi.plot(eofs.sel(mode=1), eofs.pres, linewidth=1.5)
+    inset_axes.append(axi)
+
+    axi = ax.inset_axes([-np.sqrt(2)-0.2-0.5, np.sqrt(2)+0.2, 0.5, 0.5], transform=ax.transData)
+    axi.plot(-1*eofs.sel(mode=0)+eofs.sel(mode=1), eofs.pres, linewidth=1.5)
+    inset_axes.append(axi)
+
+    axi = ax.inset_axes([-2.5,-0.25, 0.5, 0.5], transform=ax.transData)
+    axi.plot(-1*eofs.sel(mode=0), eofs.pres, linewidth=1.5)
+    inset_axes.append(axi)
+
+    axi = ax.inset_axes([-np.sqrt(2)-0.2-0.5, -np.sqrt(2)-0.2-0.5, 0.5, 0.5], transform=ax.transData)
+    axi.plot(-1*eofs.sel(mode=0)-eofs.sel(mode=1), eofs.pres, linewidth=1.5)
+    inset_axes.append(axi)
+
+    axi = ax.inset_axes([-0.25,-2.5, 0.5, 0.5], transform=ax.transData)
+    axi.plot(-1*eofs.sel(mode=1), eofs.pres, linewidth=1.5)
+    inset_axes.append(axi)
+
+    axi = ax.inset_axes([np.sqrt(2)+0.2, -np.sqrt(2)-0.2-0.5, 0.5, 0.5], transform=ax.transData)
+    axi.plot(eofs.sel(mode=0)-eofs.sel(mode=1), eofs.pres, linewidth=1.5)
+    inset_axes.append(axi)
+
+    for axi in inset_axes:
+        axi.invert_yaxis()
+        axi.set_yscale('log')
+        axi.set_yticks([])
+        axi.set_xticks([])
+        axi.minorticks_off()
+        axi.set_xlim(-1.1,1.1)
+        axi.axvline(0, color='black')
+        axi.set_facecolor('lightgrey')
+
+    ax.set_title(ptitle, fontsize=24, fontweight='semibold')
+    fig.set_size_inches(10,10)
+
+    plt.savefig(outfile, dpi=150, facecolor='white')
+
+
+def plot_u_timeseries(obs_dt,obs_u,fcst_dt,fcst_u,plot_title,outfile):
+
+    fig = plt.figure(figsize=(16,8))
+    ax = fig.add_subplot(1,1,1)
+    ax.plot(fcst_dt,fcst_u,'r',linewidth = 2,label='model')
+    ax.plot(obs_dt,obs_u,'b',linewidth = 2,label='obs')
+    ax.set_title(plot_title,fontsize=20)
+    ax.set_xlabel('Date',fontsize=14)
+    ax.set_ylabel('Zonal Mean U (m/s)',fontsize=14)
+    ax.xaxis.set_major_locator(mticker.MultipleLocator(15))
+    plt.legend()
+    plt.savefig(outfile)
+    #plt.savefig(outfile,bbox_inches='tight')