{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Example 2: Analysis of ERA5 profiles\n",
    "====================\n",
    "\n",
    "Using hourly data from ERA5 (<strong>ea</strong>) model-level (<strong>ml</strong>) and surface (<strong>sfc</strong>) analyses (<strong>an</strong>).\n",
    "\n",
    "This is for single-point profile data in <strong>.nc</strong> format - here the LAFE field campaign site at the DOE-ARM Southern Great Plains facility in Oklahoma.\n",
    "\n",
    "Calculates thermodynamic variables, indices for profile, and L-A metrics.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datetime import datetime\n",
    "\n",
    "import sys\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "from metpy.plots import SkewT\n",
    "from metpy.units import units\n",
    "import metpy.calc as calc\n",
    "import metpy.constants as c\n",
    "#import metpy.interpolate as interp\n",
    "import numpy as np\n",
    "import pandas as pd # Gonna use dataframes to manage profile data\n",
    "from netCDF4 import Dataset, date2index"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Path to ERA5 data. \n",
    "\n",
    "Set requested time using a datetime object (in range of MJJAS 2011-2019)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#epath = '/Users/dirmeyer/MyFiles/Data Sets/CLASP/SGP_hourly/'\n",
    "epath = '/scratch/dirmeyer/clasp/'\n",
    "dt = datetime(2017, 8, 19, 12) # YYYY, M, D, H(UTC) for \"sounding\"\n",
    "\n",
    "# This is a table of the pseudo-levels for ERA5 model-level data \n",
    "plevels = pd.read_csv(epath+'ERA5_levels.csv') \n",
    "plevs = plevels.loc[::-1,'plevel'].tolist() * units.pascals\n",
    "ptops = plevels.loc[::-1,'ptop'].tolist() * units.pascals\n",
    "#print(ptops)\n",
    "\n",
    "# Surface and model-level data files for indicated year\n",
    "sfcfile = epath+'ea_an_sfc_lafe_'+str(dt.year)+'0501_'+str(dt.year)+'0930.nc'\n",
    "fcfile = epath+'ea_fc_sfc_lafe_'+str(dt.year)+'0501_'+str(dt.year)+'0930.nc'\n",
    "mlfile = epath+'ea_an_ml_lafe_'+str(dt.year)+'0501_'+str(dt.year)+'0930.nc'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract Data\n",
    "----------------\n",
    "\n",
    "NetCDF data files; <strong>sfc</strong> has only a time dimension that varies, although there are 1-D latitude and longitude dimensions.\n",
    "\n",
    "The <strong>ml</strong> file has varying time and pseudo-pressure levels that need to be scaled by surface pressure (sp) to get real levels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MetPy does not define a constant for mean sea level pressure - why not?\n",
    "mean_slp = 101325.0 * units.pascal\n",
    "\n",
    "# Read sounding data based on time (dt)\n",
    "sfcdata = Dataset(sfcfile,\"r\")\n",
    "fcdata = Dataset(fcfile,\"r\")\n",
    "mldata = Dataset(mlfile,\"r\")\n",
    "#print (mldata)\n",
    "\n",
    "# Find the time index for the requested date\n",
    "all_times = mldata.variables['time']\n",
    "dindex = date2index(dt,all_times)\n",
    "#print(dindex)\n",
    "\n",
    "\n",
    "mps=units.m / units.sec\n",
    "kgkg=units.kg / units.kg\n",
    "# Standard sounding data \n",
    "ml_t = np.squeeze(mldata.variables[\"t\"][dindex]) * units.kelvin\n",
    "ml_u = np.squeeze(mldata.variables[\"u\"][dindex]) * mps #units.m / units.sec\n",
    "ml_v = np.squeeze(mldata.variables[\"v\"][dindex]) * mps # units.m / units.sec\n",
    "ml_q = np.squeeze(mldata.variables[\"q\"][dindex]) * kgkg #units.kg / units.kg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate derived fields\n",
    "---------------------------\n",
    "\n",
    "Model levels convolved with plevels and surface pressure to get actual pressure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[97563.734]]\n"
     ]
    },
    {
     "ename": "AttributeError",
     "evalue": "Neither Quantity object nor its magnitude (97563.734375) has attribute 'ndim'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m~/.conda/envs/aoes/lib/python3.6/site-packages/pint/quantity.py\u001b[0m in \u001b[0;36m__getattr__\u001b[0;34m(self, item)\u001b[0m\n\u001b[1;32m   1706\u001b[0m         \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1707\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mgetattr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_magnitude\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mitem\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1708\u001b[0m         \u001b[0;32mexcept\u001b[0m \u001b[0mAttributeError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mAttributeError\u001b[0m: 'float' object has no attribute 'ndim'",
      "\nDuring handling of the above exception, another exception occurred:\n",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-5-7d619a944d4b>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msfcdata\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvariables\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"sp\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mdindex\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;31m#sfc_p = np.squeeze(sfcdata.variables[\"sp\"][dindex]) * pa #* units.pascal\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m \u001b[0msfc_p\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msqueeze\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msfcdata\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvariables\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"sp\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mdindex\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0munits\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpascal\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      6\u001b[0m \u001b[0mml_p\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mplevs\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0msfc_p\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0mmean_slp\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0mml_top\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mptops\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0msfc_p\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0mmean_slp\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.conda/envs/aoes/lib/python3.6/site-packages/numpy/ma/core.py\u001b[0m in \u001b[0;36m__mul__\u001b[0;34m(self, other)\u001b[0m\n\u001b[1;32m   4100\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_delegate_binop\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mother\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   4101\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mNotImplemented\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 4102\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mmultiply\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mother\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   4103\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   4104\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m__rmul__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mother\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.conda/envs/aoes/lib/python3.6/site-packages/numpy/ma/core.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, a, b, *args, **kwargs)\u001b[0m\n\u001b[1;32m   1027\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1028\u001b[0m         \u001b[0;31m# Case 1. : scalar\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1029\u001b[0;31m         \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mresult\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndim\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1030\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mm\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1031\u001b[0m                 \u001b[0;32mreturn\u001b[0m \u001b[0mmasked\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.conda/envs/aoes/lib/python3.6/site-packages/pint/quantity.py\u001b[0m in \u001b[0;36m__getattr__\u001b[0;34m(self, item)\u001b[0m\n\u001b[1;32m   1709\u001b[0m             raise AttributeError(\n\u001b[1;32m   1710\u001b[0m                 \u001b[0;34m\"Neither Quantity object nor its magnitude ({}) \"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1711\u001b[0;31m                 \u001b[0;34m\"has attribute '{}'\"\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_magnitude\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mitem\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1712\u001b[0m             )\n\u001b[1;32m   1713\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mAttributeError\u001b[0m: Neither Quantity object nor its magnitude (97563.734375) has attribute 'ndim'"
     ]
    }
   ],
   "source": [
    "pa=units.pascal\n",
    "\n",
    "print(sfcdata.variables[\"sp\"][dindex])\n",
    "#sfc_p = np.squeeze(sfcdata.variables[\"sp\"][dindex]) * pa #* units.pascal\n",
    "sfc_p = np.squeeze(sfcdata.variables[\"sp\"][dindex]) * units.pascal\n",
    "ml_p = plevs * sfc_p / mean_slp\n",
    "ml_top = ptops * sfc_p / mean_slp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make profile of dew points, mixing ratios"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml_td = calc.dewpoint_from_specific_humidity(ml_q, ml_t, ml_p).to('kelvin')\n",
    "ml_w = calc.mixing_ratio_from_specific_humidity(ml_q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make profiles of various potential temperatures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml_pot = calc.potential_temperature(ml_p, ml_t)\n",
    "ml_pote = calc.equivalent_potential_temperature(ml_p, ml_t, ml_td)\n",
    "ml_potv = calc.virtual_potential_temperature(ml_p, ml_t, ml_q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make the stepwise calculations work, we have to add the surface pressure as an additional \"level\", then do the same for the other variables used.  This is because profile - based calculations must be bounded by profile levels at top and bottom to do the calculations (e.g., can't extrapolate from center of lowest model level to ground).\n",
    "\n",
    "<u>Seems that the numpy <strong>append</strong> confounds MetPy unit handling - have to redeclare the units</u>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mls_p = np.append(ml_p,sfc_p) * units.pascal\n",
    "#print(ml_p, ml_p[54])\n",
    "\n",
    "# Tack on the other \"surface\" variables\n",
    "sfc_t = np.squeeze(sfcdata.variables[\"t2m\"][dindex]) * units.kelvin # Will use this again below\n",
    "mls_t = np.append(ml_t,sfc_t) * units.kelvin\n",
    "sfc_td = np.squeeze(sfcdata.variables[\"d2m\"][dindex]) * units.kelvin # Will use this again below\n",
    "mls_td = np.append(ml_td,sfc_td) * units.kelvin\n",
    "mls_u = np.append(ml_u,np.squeeze(sfcdata.variables[\"u10\"][dindex])) * units('m/s')\n",
    "ml_v = np.append(ml_v,np.squeeze(sfcdata.variables[\"v10\"][dindex])) * units('m/s')\n",
    "\n",
    "# Gotta calcualte a surface mixing ratio\n",
    "sfc_q = calc.specific_humidity_from_dewpoint(sfc_td,sfc_p)\n",
    "sfc_w = calc.mixing_ratio_from_specific_humidity(sfc_q)\n",
    "mls_w = np.append(ml_w,sfc_w) * units.kg / units.kg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "More profile variables - these have to be built up from the bottom. Because extrapolation is forbidden, the top layer doesn't get a value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create empty profiles of new variables to fill\n",
    "#ml_thick = np.full_like(ml_p, np.nan) * units.meter\n",
    "#ml_z = np.full_like(ml_p, np.nan) * units.meter\n",
    "ml_thick = np.zeros_like(ml_p) * units.meter\n",
    "ml_z = np.zeros_like(ml_p) * units.meter\n",
    "\n",
    "bot_p = sfc_p \n",
    "for i in range(len(ml_p)-1):\n",
    "    ii = len(ml_p) - i - 1\n",
    "    top_p = ml_top[ii]\n",
    "    mid_p = ml_p[ii]\n",
    "    del_p = bot_p - top_p\n",
    "    ml_thick[ii] = calc.thickness_hydrostatic(mls_p, mls_t, mixing=mls_w, bottom=bot_p, depth=del_p)\n",
    "    ml_z[ii] = np.sum(ml_thick) - ml_thick[ii]/2.0 # Clever bit of summing to find z at middle of layer\n",
    "#    print(bot_p,mid_p,top_p,ml_z[ii])\n",
    "    bot_p = top_p\n",
    "#print(ml_z)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once we have the heights, we can calculate more profile variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml_mse = calc.moist_static_energy(ml_z, ml_t, ml_q)\n",
    "ml_mse[0] = ml_mse[0] * 0.0\n",
    "\n",
    "ml_parcel = calc.parcel_profile(np.flip(ml_p), sfc_t, sfc_td) # Note the flipping; goes from surface up.\n",
    "#print(ml_parcel)\n",
    "mls_parcel = calc.parcel_profile(np.flip(mls_p), sfc_t, sfc_td) # Note the flipping; goes from surface up.\n",
    "#print(mls_parcel)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lastly, single value (s_*) indices calculated from the profiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Some come in the surface analysis\n",
    "\"\"\"\n",
    "s_pbl_e = np.squeeze(sfcdata.variables[\"blh\"][dindex]) * units.meter\n",
    "s_cape_e = np.squeeze(sfcdata.variables[\"cape\"][dindex]) * units.joule / units.kg\n",
    "s_pw_e = np.squeeze(sfcdata.variables[\"tcw\"][dindex]) * units.kg / units.meter / units.meter\n",
    "s_tclw_e = np.squeeze(sfcdata.variables[\"tclw\"][dindex]) * units.kg / units.meter / units.meter\n",
    "s_tciw_e = np.squeeze(sfcdata.variables[\"tciw\"][dindex]) * units.kg / units.meter / units.meter\n",
    "s_pbl_f = np.squeeze(fcdata.variables[\"blh\"][dindex]) * units.meter \n",
    "s_cape_f = np.squeeze(fcdata.variables[\"cape\"][dindex]) * units.joule / units.kg\n",
    "s_cin_f = np.squeeze(fcdata.variables[\"cin\"][dindex]) * units.joule / units.kg\n",
    "s_ustar_f = np.squeeze(fcdata.variables[\"zust\"][dindex]) * units.meter / units.sec\n",
    "s_cbh_f = np.squeeze(fcdata.variables[\"cbh\"][dindex]) * units.meter \n",
    "s_sh_f = np.squeeze(fcdata.variables[\"msshf\"][dindex]) * units.meter \n",
    "s_lh_f = np.squeeze(fcdata.variables[\"mslhf\"][dindex]) * units.meter \n",
    "s_swd_f = np.squeeze(fcdata.variables[\"msdwswrf\"][dindex]) * units.watt / units.meter / units.meter\n",
    "s_lwd_f = np.squeeze(fcdata.variables[\"msdwlwrf\"][dindex]) * units.watt / units.meter / units.meter\n",
    "s_swn_f = np.squeeze(fcdata.variables[\"msnswrf\"][dindex]) * units.watt / units.meter / units.meter\n",
    "s_lwn_f = np.squeeze(fcdata.variables[\"msnlwrf\"][dindex]) * units.watt / units.meter / units.meter\n",
    "s_precip_f = np.squeeze(fcdata.variables[\"tp\"][dindex]) * units.meter / units.hour\n",
    "s_kx_f = np.squeeze(fcdata.variables[\"kx\"][dindex]) * units.kelvin\n",
    "s_ttx_f = np.squeeze(fcdata.variables[\"totalx\"][dindex]) * units.kelvin\n",
    "\"\"\"\n",
    "\n",
    "# The rest we calculate\n",
    "s_el_p, s_el_t = calc.el(mls_p, mls_t, mls_td, mls_parcel)\n",
    "s_lcl_p, s_lcl_t = calc.lcl(sfc_p, sfc_t, sfc_td, max_iters=50, eps=1e-05)\n",
    "\n",
    "s_lcl_z = (sfc_t - sfc_td) / c.dry_adiabatic_lapse_rate.to('K/m') + 2.0 * units.m\n",
    "s_lfc = calc.lfc(mls_p, mls_t, mls_td)\n",
    "s_cape, s_cin = calc.cape_cin(np.flip(mls_p), np.flip(mls_t), np.flip(mls_td), mls_parcel)\n",
    "s_pw = calc.precipitable_water(mls_td, mls_p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make Skew-T Plot\n",
    "----------------\n",
    "\n",
    "The code below makes a basic skew-T plot using the MetPy plot module\n",
    "that contains a SkewT class.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Change default to be better for skew-T\n",
    "fig = plt.figure(figsize=(14, 14))\n",
    "\n",
    "# Initiate the skew-T plot type from MetPy class loaded earlier\n",
    "skew = SkewT(fig, rotation=40)\n",
    "\n",
    "# Plot the data using normal plotting functions, in this case using\n",
    "# log scaling in Y, as dictated by the typical meteorological plot\n",
    "skew.plot(ml_p.to('hectopascal'), ml_t.to('degC'), 'crimson')\n",
    "skew.plot(ml_p.to('hectopascal'), ml_td.to('degC'), 'teal')\n",
    "skew.plot(sfc_p.to('hectopascal'), sfc_t.to('degC'), color='maroon', marker='x')\n",
    "skew.plot(sfc_p.to('hectopascal'), sfc_td.to('degC'), color='darkslategrey', marker='x')\n",
    "\n",
    "skew.plot(mls_p.to('hectopascal'), np.flip(mls_parcel).to('degC'), 'hotpink')\n",
    "skew.plot_barbs(ml_p[::3].to('hectopascal'), ml_u[::3], ml_v[::3], y_clip_radius=0.03)\n",
    "\n",
    "# Set some appropriate axes limits for x and y\n",
    "skew.ax.set_xlim(-30, 40)\n",
    "skew.ax.set_ylim(1020, 250)\n",
    "\n",
    "# Shade areas of CAPE and CIN\n",
    "skew.shade_cin(mls_p.to('hectopascal'), mls_t.to('degC'), np.flip(mls_parcel).to('degC'))\n",
    "skew.shade_cape(mls_p.to('hectopascal'), mls_t.to('degC'), np.flip(mls_parcel).to('degC'))\n",
    "\n",
    "# Add the relevant special lines to plot throughout the figure\n",
    "skew.plot_dry_adiabats(t0=np.arange(233, 533, 10) * units.K,\n",
    "                       alpha=0.25, color='orangered')\n",
    "skew.plot_moist_adiabats(t0=np.arange(233, 400, 5) * units.K,\n",
    "                         alpha=0.25, color='tab:green')\n",
    "skew.plot_mixing_lines(p=np.arange(1000, 99, -20) * units.hPa,\n",
    "                       linestyle='dotted', color='tab:blue')\n",
    "\n",
    "# Add some descriptive titles\n",
    "plt.title('ERA5 Pseudo-Sounding at LAFE site', loc='left')\n",
    "plt.title('Valid Time: {}'.format(dt), loc='right')\n",
    "\n",
    "plt.show()\n",
    "#print(\"EL:\",s_el_p, s_el_t)\n",
    "print(\"LCL:\",s_lcl_z, s_lcl_p, s_lcl_t)\n",
    "print(\"CAPE:\",s_cape,s_cape_e)\n",
    "print(\"CIN:\",s_cin)\n",
    "print(\"PW:\",s_pw,s_pw_e)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Red is CAPE area, blue is CIN; surface values are 'x'."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (aoes)",
   "language": "python",
   "name": "aoes"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}