{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Skew-T's\n",
    "## Use `functions` to loop over dates and stations\n",
    "### Demonstrate Python's `try ... except` utility"
   ]
  },
  {
   "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",
    "from dateutil.relativedelta import relativedelta\n",
    "from time import sleep\n",
    "import pandas as pd\n",
    "from siphon.simplewebservice.wyoming import WyomingUpperAir"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define a function that will retrieve the skew-t from the Wyoming data server."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def retrieveSkewt (t, site):\n",
    "    '''\n",
    "    This function obtains sounding data for a particular time and location, and draws a Skew-T\n",
    "    '''\n",
    "    \n",
    "    curtim = t.strftime(\"%y%m%d%H\")\n",
    "    print (f\"Processing {t} for {site}\")\n",
    "\n",
    "#   Query the UWyoming Data Server. Handle the situation when a data retrieval error occurs.\n",
    "    try:\n",
    "        df = WyomingUpperAir.request_data (t ,site)\n",
    "        drawSkewt(df, site, current, lev)\n",
    "    except Exception as errorMsg:\n",
    "        print (f\"Error returning sounding for {site} at {t}: {errorMsg}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define a function that will process and draw the skew-t from the Wyoming data server."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def drawSkewt(df, site, dattim, lev):\n",
    "    '''\n",
    "    drawSkewt processes and draws a Skew-t that has been downloaded from the UWyoming server.\n",
    "    '''\n",
    "    \n",
    "    raob_timeStr = dattim.strftime(\"%Y%m%d%H\")\n",
    "    raob_timeTitle = dattim.strftime(\"%H00 UTC %-d %b %Y\")\n",
    "\n",
    "### Attach units to the various columns using MetPy's `pandas_dataframe_to_unit_arrays` function, as follows:\n",
    "\n",
    "    df = pandas_dataframe_to_unit_arrays(df)\n",
    "\n",
    "# As with any Pandas Dataframe, we can select columns, so we do so now for all the relevant variables for our Skew-T. \n",
    "\n",
    "    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']\n",
    "\n",
    "# Now that U and V are units aware, we can use one of MetPy's diagnostic routines to \n",
    "# calculate the wind speed from the vector's components.\n",
    "\n",
    "    wind_speed = mpcalc.wind_speed(U, V)\n",
    "\n",
    "# Calculate the LCL\n",
    "\n",
    "    lclp, lclt = mpcalc.lcl(P[lev], T[lev], Td[lev])\n",
    "\n",
    "# Calculate the parcel profile. Here, we pass in the array of all pressure values, and specify the index of the temperature and dewpoint arrays as the starting point of the parcel path.\n",
    "\n",
    "    parcel_prof = mpcalc.parcel_profile(P[lev:], T[lev], Td[lev])\n",
    "\n",
    "# Set the limits for the x- and y-axes so they can be changed and easily re-used when desired.\n",
    "\n",
    "    pTop = 100\n",
    "    pBot = 1050\n",
    "    tMin = -60\n",
    "    tMax = 40\n",
    "\n",
    "# Define the log-p axis.\n",
    "\n",
    "    logTop = np.log10(pTop)\n",
    "    logBot = np.log10(pBot)\n",
    "    interval = np.logspace(logTop,logBot) * units('hPa')\n",
    "\n",
    "# This interval will tend to oversample points near the top of the atmosphere relative to the bottom. \n",
    "# One effect of this will be that wind barbs will be hard to read, especially at higher levels. \n",
    "# Use a resampling method from NumPy to deal with this.\n",
    "\n",
    "    idx = mpcalc.resample_nn_1d(P, interval)\n",
    "\n",
    "# Plot the skew-t and save it as a PNG file.\n",
    "\n",
    "# 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",
    "\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[lev:], parcel_prof, 'k', linestyle='--', linewidth=2)\n",
    "\n",
    "# Shade areas of CAPE and CIN\n",
    "    skew.shade_cin(P[lev:], T[lev:], parcel_prof)\n",
    "    skew.shade_cape(P[lev:], T[lev:], 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",
    "\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\"{site}_{raob_timeStr}_skewt.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Set the stations, and the date/time range over which you wish to produce skew-t's."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stationList = ['ALB', 'BUF']\n",
    "stationList = ['ALB']\n",
    "start = datetime(2024,4,2,12)\n",
    "\n",
    "end = datetime(2024,4,2,12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Set the level from which you wish to raise the parcel from (0 = lowest level)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lev = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Loop over the times; retrieve and draw the Skew-t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "current = start\n",
    "\n",
    "while (current <= end):\n",
    "    for station in stationList:\n",
    "        # Sleep 5 seconds as an attempt to avoid UWyoming 503 Server Busy errors\n",
    "        sleep(5)\n",
    "        df = retrieveSkewt(current, station)\n",
    "        current = current + relativedelta(hours=12)"
   ]
  },
  {
   "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>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>The Wyoming server does not handle simultaneous requests well. If indeed there was a sounding taken for the site and time, usually the request will eventually work. </li>\n",
    "    </ul>\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"admonition alert alert-success\">\n",
    "    <p class=\"admonition-title\" style=\"font-weight:bold\">Things to try:</p>\n",
    "    <ol>\n",
    "        <li><b>Basic: </b>Modify the stations and times to fit your final project</li>\n",
    "        <li><b>Basic: </b>Use comment lines to modify the <code>drawSkewt</code>function to change or omit what's drawn on the Skew-t (e.g., omit the DGZ if irrelevant to your case)</li>\n",
    "        <li><b>Advanced: </b>Modify the <code>drawSkewt</code>function so what's drawn or omitted on the Skew-t is passed as an argument (<i>hint:</i> use a <b>dictionary</b>)</li>\n",
    "        <li><b>Advanced: </b>Refine the <code>retrieveSkewt</code>function so error exceptions are handled better (e.g., automatically attempt to retrieve a sounding if a prior attempt fails with a <code>503</code> error</li>\n",
    "        <li><b>Advanced: </b>Specify a particular pressure level, rather than a row index, and for each sounding, find the row index that is closest to the specified pressure value</li>\n",
    "    </ol>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References:\n",
    "1. [Python Functions](https://www.geeksforgeeks.org/python-functions/)\n",
    "1. [Try ... except in Python](https://www.freecodecamp.org/news/python-print-exception-how-to-try-except-print-an-error/)\n",
    "1. [Python `sleep` function](https://www.programiz.com/python-programming/time/sleep)"
   ]
  }
 ],
 "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
}