pages/gallery/Observational_Data_Cross_Section.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "Real Data Cross-Section Example\n",
    "===============================\n",
    "\n",
    "Cross-section using real data from soundings.\n",
    "\n",
    "This example uses actual soundings to create a cross-section. There are\n",
    "two functions defined to help interpolate radiosonde observations, which\n",
    "won’t all be at the same level, to a standard grid. The vertical\n",
    "interpolation assumes a log-linear relationship. Each radisosonde\n",
    "vertical profile is interpolated first, then the\n",
    "``scipy.interpolate.griddata`` function is used to generate a full 2D\n",
    "(x, p) grid between each station. Pyproj is used to calculate the\n",
    "distance between each station and the standard atmosphere is used to\n",
    "convert the elevation of each station to a pressure value for plotting\n",
    "purposes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datetime import datetime\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import metpy.calc as mpcalc\n",
    "from metpy.units import units\n",
    "import numpy as np\n",
    "from pyproj import Geod\n",
    "from scipy.interpolate import griddata\n",
    "from scipy.ndimage import gaussian_filter\n",
    "from siphon.simplewebservice.wyoming import WyomingUpperAir"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_marker": "######################################################################"
   },
   "source": [
    "Vertical Interpolation Function\n",
    "-------------------------------\n",
    "\n",
    "Function interpolates to given pressure level data to set grid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 1
   },
   "outputs": [],
   "source": [
    "def vertical_interpolate(vcoord_data, interp_var, interp_levels):\n",
    "    \"\"\"A function to interpolate sounding data from each station to\n",
    "    every millibar. Assumes a log-linear relationship.\n",
    "\n",
    "    Input\n",
    "    -----\n",
    "    vcoord_data : A 1D array of vertical level values (e.g., pressure from a radiosonde)\n",
    "    interp_var : A 1D array of the variable to be interpolated to all pressure levels\n",
    "    vcoord_interp_levels : A 1D array containing veritcal levels to interpolate to\n",
    "\n",
    "    Return\n",
    "    ------\n",
    "    interp_data : A 1D array that contains the interpolated variable on the interp_levels\n",
    "    \"\"\"\n",
    "\n",
    "    # Make veritcal coordinate data and grid level log variables\n",
    "    lnp = np.log(vcoord_data)\n",
    "    lnp_intervals = np.log(interp_levels)\n",
    "\n",
    "    # Use numpy to interpolate from observed levels to grid levels\n",
    "    interp_data = np.interp(lnp_intervals[::-1], lnp[::-1], interp_var[::-1])[::-1]\n",
    "\n",
    "    # Mask for missing data (generally only near the surface)\n",
    "    mask_low = interp_levels > vcoord_data[0]\n",
    "    mask_high = interp_levels < vcoord_data[-1]\n",
    "    interp_data[mask_low] = interp_var[0]\n",
    "    interp_data[mask_high] = interp_var[-1]\n",
    "\n",
    "    return interp_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_marker": "######################################################################"
   },
   "source": [
    "Radiosonde Observation Interpolation Function\n",
    "---------------------------------------------\n",
    "\n",
    "This function interpolates given radiosonde data into a 2D array for all\n",
    "meteorological variables given in dataframe. Returns a dictionary that\n",
    "will have requesite data for plotting a cross section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 1
   },
   "outputs": [],
   "source": [
    "def radisonde_cross_section(stns, data, start=1000, end=100, step=10):\n",
    "    \"\"\"This function takes a list of radiosonde observation sites with a\n",
    "    dictionary of Pandas Dataframes with the requesite data for each station.\n",
    "\n",
    "    Input\n",
    "    -----\n",
    "    stns : List of statition three-letter identifiers\n",
    "    data : A dictionary of Pandas Dataframes containing the radiosonde observations\n",
    "    for the stations\n",
    "    start : interpolation start value, optional (default = 1000 hPa)\n",
    "    end : Interpolation end value, optional (default = 100 hPa)\n",
    "    step : Interpolation interval, option (default = 10 hPa)\n",
    "\n",
    "    Return\n",
    "    ------\n",
    "    cross_section : A dictionary that contains the following variables\n",
    "\n",
    "        grid_data : An interpolated grid with 100 points between the first and last station,\n",
    "        with the corresponding number of vertical points based on start, end, and interval\n",
    "        (default is 90)\n",
    "        obs_distance : An array of distances between each radiosonde observation location\n",
    "        x_grid : A 2D array of horizontal direction grid points\n",
    "        p_grid : A 2D array of vertical pressure levels\n",
    "        ground_elevation : A representation of the terrain between radiosonde observation sites\n",
    "        based on the elevation of each station converted to pressure using the standard\n",
    "        atmosphere\n",
    "\n",
    "    \"\"\"\n",
    "    # Set up vertical grid, largest value first (high pressure nearest surface)\n",
    "    vertical_levels = np.arange(start, end-1, -step)\n",
    "\n",
    "    # Number of vertical levels and stations\n",
    "    plevs = len(vertical_levels)\n",
    "    nstns = len(stns)\n",
    "\n",
    "    # Create dictionary of interpolated values and include neccsary attribute data\n",
    "    # including lat, lon, and elevation of each station\n",
    "    lats = []\n",
    "    lons = []\n",
    "    elev = []\n",
    "    keys = data[stns[0]].keys()[:8]\n",
    "    tmp_grid = dict.fromkeys(keys)\n",
    "\n",
    "    # Interpolate all variables for each radiosonde observation\n",
    "    # Temperature, Dewpoint, U-wind, V-wind\n",
    "    for key in tmp_grid.keys():\n",
    "        tmp_grid[key] = np.empty((nstns, plevs))\n",
    "        for station, loc in zip(stns, range(nstns)):\n",
    "            if key == 'pressure':\n",
    "                lats.append(data[station].latitude[0])\n",
    "                lons.append(data[station].longitude[0])\n",
    "                elev.append(data[station].elevation[0])\n",
    "                tmp_grid[key][loc, :] = vertical_levels\n",
    "            else:\n",
    "                tmp_grid[key][loc, :] = vertical_interpolate(\n",
    "                    data[station]['pressure'].values, data[station][key].values,\n",
    "                    vertical_levels)\n",
    "\n",
    "    # Compute distance between each station using Pyproj\n",
    "    g = Geod(ellps='sphere')\n",
    "    _, _, dist = g.inv(nstns*[lons[0]], nstns*[lats[0]], lons[:], lats[:])\n",
    "\n",
    "    # Compute sudo ground elevation in pressure from standard atmsophere and the elevation\n",
    "    # of each station\n",
    "    ground_elevation = mpcalc.height_to_pressure_std(np.array(elev) * units('meters'))\n",
    "\n",
    "    # Set up grid for 2D interpolation\n",
    "    grid = dict.fromkeys(keys)\n",
    "    x = np.linspace(dist[0], dist[-1], 100)\n",
    "    nx = len(x)\n",
    "\n",
    "    pp, xx = np.meshgrid(vertical_levels, x)\n",
    "    pdist, ddist = np.meshgrid(vertical_levels, dist)\n",
    "\n",
    "    # Interpolate to 2D grid using scipy.interpolate.griddata\n",
    "    for key in grid.keys():\n",
    "        grid[key] = np.empty((nx, plevs))\n",
    "        grid[key][:] = griddata((ddist.flatten(), pdist.flatten()),\n",
    "                                tmp_grid[key][:].flatten(),\n",
    "                                (xx, pp),\n",
    "                                method='cubic')\n",
    "\n",
    "    # Gather needed data in dictionary for return\n",
    "    cross_section = {'grid_data': grid, 'obs_distance': dist,\n",
    "                     'x_grid': xx, 'p_grid': pp, 'elevation': ground_elevation}\n",
    "    return cross_section"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_marker": "######################################################################"
   },
   "source": [
    "Stations and Time\n",
    "-----------------\n",
    "\n",
    "Select cross section stations by creating a list of three-letter\n",
    "identifiers and choose a date by creating a datetime object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A roughly east-west cross section\n",
    "stn_list = ['DNR', 'LBF', 'OAX', 'DVN', 'DTX', 'BUF']\n",
    "\n",
    "# Set a date and hour of your choosing\n",
    "date = datetime(2019, 6, 1, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_marker": "######################################################################"
   },
   "source": [
    "Get Radiosonde Data\n",
    "-------------------\n",
    "\n",
    "This example is built around the data from the University of Wyoming\n",
    "sounding archive and using the Siphon package to remotely access that\n",
    "data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up empty dictionary to fill with Wyoming Sounding data\n",
    "df = {}\n",
    "\n",
    "# Loop over stations to get data and put into dictionary\n",
    "for station in stn_list:\n",
    "    df[station] = WyomingUpperAir.request_data(date, station)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_marker": "######################################################################"
   },
   "source": [
    "Create Interpolated fields\n",
    "--------------------------\n",
    "\n",
    "Use the function ``radisonde_cross_section`` to generate the 2D grid (x,\n",
    "p) for all radiosonde variables including, Temperature, Dewpoint,\n",
    "u-component of the wind, and v-component of the wind."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xsect = radisonde_cross_section(stn_list, df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_marker": "######################################################################"
   },
   "source": [
    "Calculate Variables for Plotting\n",
    "--------------------------------\n",
    "\n",
    "Use MetPy to calculate common variables for plotting a cross section,\n",
    "specifically potential temperature and mixing ratio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "potemp = mpcalc.potential_temperature(\n",
    "    xsect['p_grid'] * units('hPa'), xsect['grid_data']['temperature'] * units('degC'))\n",
    "\n",
    "relhum = mpcalc.relative_humidity_from_dewpoint(\n",
    "    xsect['grid_data']['temperature'] * units('degC'),\n",
    "    xsect['grid_data']['dewpoint'] * units('degC'))\n",
    "\n",
    "mixrat = mpcalc.mixing_ratio_from_relative_humidity(relhum,\n",
    "                                                    xsect['grid_data']['temperature'] *\n",
    "                                                    units('degC'),\n",
    "                                                    xsect['p_grid'] * units('hPa'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_marker": "######################################################################"
   },
   "source": [
    "Plot Cross Section\n",
    "------------------\n",
    "\n",
    "Use standard Matplotlib to plot the now 2D cross section grid using the\n",
    "data from xsect and those calculated above. Additionally, the actualy\n",
    "radiosonde wind observations are plotted as barbs on this plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Start Figure, set big size for cross section\n",
    "fig = plt.figure(figsize=(17, 11))\n",
    "\n",
    "# Specify plotting axis (single panel)\n",
    "ax = plt.subplot(111)\n",
    "\n",
    "# Set y-scale to be log since pressure decreases exponentially with height\n",
    "ax.set_yscale('log')\n",
    "\n",
    "# Set limits, tickmarks, and ticklabels for y-axis\n",
    "ax.set_ylim([1030, 101])\n",
    "ax.set_yticks(range(1000, 101, -100))\n",
    "ax.set_yticklabels(range(1000, 101, -100))\n",
    "\n",
    "# Invert the y-axis since pressure decreases with increasing height\n",
    "ax.yaxis_inverted()\n",
    "\n",
    "# Plot the sudo elevation on the cross section\n",
    "ax.fill_between(xsect['obs_distance'], xsect['elevation'].m, 1030,\n",
    "                where=xsect['elevation'].m <= 1030, facecolor='lightgrey',\n",
    "                interpolate=True, zorder=10)\n",
    "# Don't plot xticks\n",
    "plt.xticks([], [])\n",
    "\n",
    "# Plot wind barbs for each sounding location\n",
    "for stn, stn_name in zip(range(len(stn_list)), stn_list):\n",
    "    ax.axvline(xsect['obs_distance'][stn], ymin=0, ymax=1,\n",
    "               linewidth=2, color='blue', zorder=11)\n",
    "    ax.text(xsect['obs_distance'][stn], 1100, stn_name, ha='center', color='blue')\n",
    "    ax.barbs(xsect['obs_distance'][stn], df[stn_name]['pressure'][::2],\n",
    "             df[stn_name]['u_wind'][::2, None],\n",
    "             df[stn_name]['v_wind'][::2, None], zorder=15)\n",
    "\n",
    "# Plot smoothed potential temperature grid (K)\n",
    "cs = ax.contour(xsect['x_grid'], xsect['p_grid'], gaussian_filter(\n",
    "    potemp, sigma=1.0), range(0, 500, 5), colors='red')\n",
    "ax.clabel(cs, fmt='%i')\n",
    "\n",
    "# Plot smoothed mixing ratio grid (g/kg)\n",
    "cs = ax.contour(xsect['x_grid'], xsect['p_grid'], gaussian_filter(\n",
    "    mixrat*1000, sigma=2.0), range(0, 41, 2), colors='tab:green')\n",
    "ax.clabel(cs, fmt='%i')\n",
    "\n",
    "# Add some informative titles\n",
    "plt.title('Cross-Section from DNR to BUF Potential Temp. '\n",
    "          '(K; red) and Mix. Rat. (g/kg; green)', loc='left')\n",
    "plt.title(date, loc='right');"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "-all",
   "main_language": "python",
   "notebook_metadata_filter": "-all"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}