{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Skew-T's\n",
    "## Requesting Upper Air Data from Siphon and Plotting a Skew-T with MetPy\n",
    "### Based on <a href = \"http://www.unidata.ucar.edu\"> Unidata's </a> <i> MetPy Monday </i> \\# 15 -17 YouTube links: <br>https://www.youtube.com/watch?v=OUTBiXEuDIU; <br> https://www.youtube.com/watch?v=oog6_b-844Q; <br> https://www.youtube.com/watch?v=b0RsN9mCY5k"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import the libraries we need"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.axes_grid1.inset_locator import inset_axes\n",
    "import metpy.calc as mpcalc\n",
    "from metpy.plots import Hodograph, SkewT\n",
    "from metpy.units import units, pandas_dataframe_to_unit_arrays\n",
    "from datetime import datetime, timedelta\n",
    "import pandas as pd\n",
    "from siphon.simplewebservice.wyoming import WyomingUpperAir"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Next we use methods from the datetime library to determine the latest time, and then use it to select the most recent date and time that we'll use for our selected Skew-T."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "now = datetime.utcnow()\n",
    "\n",
    "curr_year = now.year\n",
    "curr_month = now.month\n",
    "curr_day = now.day\n",
    "curr_hour = now.hour\n",
    "\n",
    "print(\"Current time is: \", now)\n",
    "print(\"Current year, month, date, hour: \" ,curr_year, curr_month,curr_day,curr_hour)\n",
    "\n",
    "if (curr_hour > 13) :\n",
    "    raob_hour = 12\n",
    "    hr_delta = curr_hour - 12\n",
    "elif (curr_hour > 1):\n",
    "    raob_hour = 0\n",
    "    hr_delta = curr_hour - 0\n",
    "else:\n",
    "    raob_hour = 12\n",
    "    hr_delta = curr_hour + 12\n",
    "\n",
    "raob_time = now - timedelta(hours=hr_delta)\n",
    "\n",
    "raob_year = raob_time.year\n",
    "raob_month = raob_time.month\n",
    "raob_day = raob_time.day\n",
    "raob_hour = raob_time.hour\n",
    "\n",
    "print(\"Time of RAOB is: \", raob_year, raob_month, raob_day, raob_hour)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Construct a datetime object that we will use in our query to the data server. Note what it looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query_date = datetime(raob_year,raob_month,raob_day,raob_hour)\n",
    "query_date"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### If desired, we can choose a past date and time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#current = True\n",
    "current = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if (current == False):\n",
    "    query_date = datetime(2012, 10, 29, 12)\n",
    "    \n",
    "raob_timeStr = query_date.strftime(\"%Y-%m-%d %H UTC\")\n",
    "raob_timeTitle = query_date.strftime(\"%H00 UTC %-d %b %Y\")\n",
    "print(raob_timeStr)\n",
    "print(raob_timeTitle)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Select our station and query the data server."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "station = 'ALB'\n",
    "\n",
    "df = WyomingUpperAir.request_data (query_date,station)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What does the returned data file look like? Well, it looks like a Pandas Dataframe!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### We need to deal with an issue where two levels may be at the same pressure, but different heights.  While rare, it won't allow a parcel path to be calculated later in this notebook:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check to see if this is occurring by looking more closely at the data.  Let's look at the last 30 rows of this dataframe (this issue usually occurs near the very top, or lowest pressures, of our sounding data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.tail(30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If any rows have the same pressure, this will cause an error when we call MetPy's `parcel_path` function. Likely this will get resolved in a future release of MetPy, but for now, let's use a Pandas function to eliminate rows whose pressure is the same as a previous row."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.drop_duplicates(inplace=True,subset='pressure',ignore_index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['pressure']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### As with any Pandas Dataframe, we can select columns, so we do so now for all the relevant variables for our Skew-T. Note now that we are using MetPy's <i>units</i> library to attach units to each variable's array of values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = df['pressure'].values * units(df.units['pressure'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Let's split up the code above, to see more detail:\n",
    "df['pressure'].values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The University of Wyoming sounding reader includes a <i>units</i> attribute, which produces a Python dictionary of units specific to Wyoming sounding data.  What does it look like? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.units"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#From the dictionary, we can specify pressure units.  In this case, hectopascal (hPa):\n",
    "units(df.units['pressure'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Putting it all together...\n",
    "p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#set this for use later:\n",
    "p_units=units(df.units['pressure'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Set remaining variables and apply units:\n",
    "Z = df['height'].values * units(df.units['height'])\n",
    "T = df['temperature'].values * units(df.units['temperature'])\n",
    "Td = df['dewpoint'].values * units(df.units['dewpoint'])\n",
    "u = df['u_wind'].values * units(df.units['u_wind'])\n",
    "v = df['v_wind'].values * units(df.units['v_wind'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: we can also attach units to the various columns using MetPy's `pandas_dataframe_to_unit_arrays` function, as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pandas_dataframe_to_unit_arrays(df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that `u` and `v` are *units aware*, we can use one of MetPy's diagnostic routines to calculate the wind speed from the vector's components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wind_speed = mpcalc.wind_speed(u,v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thermodynamic Calculations\n",
    "--------------------------\n",
    "\n",
    "Often times we will want to calculate some thermodynamic parameters of a\n",
    "sounding. The MetPy calc module has many such calculations already implemented!\n",
    "\n",
    "* **Lifting Condensation Level (LCL)** - The level at which an air parcel's\n",
    "  relative humidity becomes 100% when lifted along a dry adiabatic path.\n",
    "* **Parcel Path** - Path followed by a hypothetical parcel of air, beginning\n",
    "  at the surface temperature/pressure and rising dry adiabatically until\n",
    "  reaching the LCL, then rising moist adiabatially."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate the LCL. Note that the units will appear at the end when we print them out!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lclp, lclt = mpcalc.lcl(p[0], T[0], Td[0])\n",
    "print ('LCL Pressure = ',lclp, 'LCL Temperature = ',lclt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate the parcel profile. Here, we pass in the array of all pressure values, and specify the first (index 0, remember Python's zero-based indexing) index of the temperature and dewpoint arrays as the starting point of the parcel path."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "parcel_prof = mpcalc.parcel_profile(p, T[0], Td[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What does the parcel profile look like? Well, it's just an array of all the temperatures of the parcel at each of the pressure levels as it ascends."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "parcel_prof"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set the limits for the x- and y-axes so they can be changed and easily re-used when desired."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pTop = 100\n",
    "pBot = 1050\n",
    "tMin = -60\n",
    "tMax = 40"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Since we are plotting a skew-T, define the log-p axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logTop = np.log10(pTop)\n",
    "logBot = np.log10(pBot)\n",
    "interval = np.logspace(logTop,logBot) * p_units"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### This interval will tend to oversample points near the top of the atmosphere relative to the bottom. One effect of this will be that wind barbs will be hard to read, especially at higher levels. Use a resampling method from NumPy to deal with this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx = mpcalc.resample_nn_1d(p, interval)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now we are ready to create our skew-T. Plot it and save it as a PNG file."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"admonition alert alert-warning\">\n",
    "    <p class=\"admonition-title\" style=\"font-weight:bold\">Note:</p>\n",
    "    Annoyingly, when specifying the color of the dry/moist adiabats, and mixing ratio lines, the parameter must be <b>color<i>s</i></b>, not <b>color</b>. In other Matplotlib-related code lines in the cell below where we set the color, such as <code>skew.ax.axvline</code> ,<b> color</b> must be singular.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a new figure. The dimensions here give a good aspect ratio\n",
    "fig = plt.figure(figsize=(9, 9))\n",
    "skew = SkewT(fig, rotation=30)\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(p, T, 'r', linewidth=3)\n",
    "skew.plot(p, Td, 'g', linewidth=3)\n",
    "skew.plot_barbs(p[idx], u[idx], v[idx])\n",
    "skew.ax.set_ylim(pBot, pTop)\n",
    "skew.ax.set_xlim(tMin, tMax)\n",
    "\n",
    "# Plot LCL as black dot\n",
    "skew.plot(lclp, lclt, 'ko', markerfacecolor='black')\n",
    "\n",
    "# Plot the parcel profile as a black line\n",
    "skew.plot(p, parcel_prof, 'k', linestyle='--', linewidth=2)\n",
    "\n",
    "# Shade areas of CAPE and CIN\n",
    "skew.shade_cin(p, T, parcel_prof)\n",
    "skew.shade_cape(p, T, parcel_prof)\n",
    "\n",
    "# Plot a zero degree isotherm\n",
    "skew.ax.axvline(0, color='c', linestyle='-', linewidth=2)\n",
    "\n",
    "# Add the relevant special lines\n",
    "skew.plot_dry_adiabats()\n",
    "skew.plot_moist_adiabats(colors='#58a358', linestyle='-')\n",
    "skew.plot_mixing_lines()\n",
    "\n",
    "plt.title(\"RAOB for \" + station + ' valid at: '+ raob_timeTitle)\n",
    "\n",
    "# Create a hodograph\n",
    "# Create an inset axes object that is 30% width and height of the\n",
    "# figure and put it in the upper right hand corner.\n",
    "ax_hod = inset_axes(skew.ax, '30%', '30%', loc=1)\n",
    "h = Hodograph(ax_hod, component_range=100.)\n",
    "h.add_grid(increment=30)\n",
    "#h.plot_colormapped(u, v, df['height'].values)  # Plot a line colored by a parameter\n",
    "h.plot_colormapped(u, v, df['height'])\n",
    "plt.savefig (station + raob_timeStr + '_skewt.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's how the plot would appear if we hadn't done the resampling. Note that the wind barbs are plotted at every level, and are harder to read:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a new figure. The dimensions here give a good aspect ratio\n",
    "fig = plt.figure(figsize=(9, 9))\n",
    "skew = SkewT(fig, rotation=30)\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(p, T, 'r', linewidth=3)\n",
    "skew.plot(p, Td, 'g', linewidth=3)\n",
    "skew.plot_barbs(p, u, v)\n",
    "skew.ax.set_ylim(pBot, pTop)\n",
    "skew.ax.set_xlim(tMin, tMax)\n",
    "\n",
    "# Plot LCL as black dot\n",
    "skew.plot(lclp, lclt, 'ko', markerfacecolor='black')\n",
    "\n",
    "# Plot the parcel profile as a black line\n",
    "skew.plot(p, parcel_prof, 'k', linestyle='--', linewidth=2)\n",
    "\n",
    "# Shade areas of CAPE and CIN\n",
    "skew.shade_cin(p, T, parcel_prof)\n",
    "skew.shade_cape(p, T, parcel_prof)\n",
    "\n",
    "# Plot a zero degree isotherm\n",
    "skew.ax.axvline(0, color='c', linestyle='-', linewidth=2)\n",
    "\n",
    "# Add the relevant special lines\n",
    "skew.plot_dry_adiabats()\n",
    "skew.plot_moist_adiabats(colors='#58a358', linestyle='-')\n",
    "skew.plot_mixing_lines()\n",
    "\n",
    "plt.title(\"RAOB for \" + station + ' valid at: '+ raob_timeTitle)\n",
    "\n",
    "# Create a hodograph\n",
    "# Create an inset axes object that is 30% width and height of the\n",
    "# figure and put it in the upper right hand corner.\n",
    "ax_hod = inset_axes(skew.ax, '30%', '30%', loc=1)\n",
    "h = Hodograph(ax_hod, component_range=100.)\n",
    "h.add_grid(increment=30)\n",
    "h.plot_colormapped(u, v, df['height'])  # Plot a line colored by a parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 August 2022 Environment",
   "language": "python",
   "name": "aug22"
  },
  "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.10.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}