{
 "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": {
    "tags": []
   },
   "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": {
    "tags": []
   },
   "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(f\"Current time is: {now}\")\n",
    "print(f\"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(f\"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": {
    "tags": []
   },
   "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": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#current = True\n",
    "current = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "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": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "station = 'ALB'\n",
    "\n",
    "df = WyomingUpperAir.request_data (query_date,station)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"admonition alert alert-warning\">\n",
    "    <p class=\"admonition-title\" style=\"font-weight:bold\">Possible data retrieval errors:</p>\n",
    "    <ul>\n",
    "        <li>If you get a <code>No data available for YYYY-MM-DD HHZ for station XXX </code>, you may have specified an incorrect date or location, or there just might not have been a successful RAOB at that date and time.</code></li>\n",
    "        <li>If you get a <code>503 Server Error</code>, try requesting the sounding again. The Wyoming server does not handle simultaneous requests well. Usually the request will eventually work. </li>\n",
    "    </ul>\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "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": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Examine one of the columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df['pressure']"
   ]
  },
  {
   "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": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df.units"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Attach units to the various columns using MetPy's `pandas_dataframe_to_unit_arrays` function, as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df = pandas_dataframe_to_unit_arrays(df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that units are attached, how does the `pressure` column look?"
   ]
  },
  {
   "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. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P = df['pressure']\n",
    "Z = df['height']\n",
    "T = df['temperature']\n",
    "Td = df['dewpoint']\n",
    "u = df['u_wind']\n",
    "v = df['v_wind']"
   ]
  },
  {
   "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": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "wind_speed = mpcalc.wind_speed(u,v)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "wind_speed"
   ]
  },
  {
   "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": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "lclp, lclt = mpcalc.lcl(P[0], T[0], Td[0])\n",
    "print (f'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": {
    "tags": []
   },
   "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": {
    "tags": []
   },
   "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": {
    "tags": []
   },
   "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": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "logTop = np.log10(pTop)\n",
    "logBot = np.log10(pBot)\n",
    "interval = np.logspace(logTop,logBot) * units('hPa')"
   ]
  },
  {
   "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": {
    "tags": []
   },
   "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": {
    "tags": []
   },
   "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",
    "# Shade areas between a certain temperature (in this case, the DGZ)\n",
    "skew.shade_area(P.m,-18,-12,color=['purple'])\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(f\"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'])\n",
    "plt.savefig (f\"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": {
    "tags": []
   },
   "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(f\"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": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"admonition alert alert-success\">\n",
    "    <p class=\"admonition-title\" style=\"font-weight:bold\">What's next:</p>\n",
    "    Try different RAOB locations and dates (for example, draw a skew-t relevant to your class project). <br>\n",
    "    Here is a link to <a href=\"https://www.weather.gov/upperair/nws_upper\">a map of NWS RAOB sites</a>.<br>\n",
    "    Depending on the situation, you may want to omit plotting the DGZ, LCL, parcel path, CAPE/CIN, freezing line, and/or move the location of the hodograph.<br>\n",
    "</div>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 Jan. 2024 Environment",
   "language": "python",
   "name": "jan24"
  },
  "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.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}