{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# NCL_skewt_3_2.py\nThis script illustrates the following concepts:\n    - Drawing Skew-T plots\n    - Thinning the wind barbs in a Skew-T plot\n    - Customizing the background of a Skew_T plot\n    \nSee following URLs to see the reproduced NCL plot & script:\n    - Original NCL script: https://www.ncl.ucar.edu/Applications/Scripts/skewt_3.ncl\n    - Original NCL plot: https://www.ncl.ucar.edu/Applications/Images/skewt_3_2_lg.png\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Import packages:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport matplotlib.pyplot as plt\nimport matplotlib.lines as mlines\nimport pandas as pd\nimport metpy.calc as mpcalc\nfrom metpy.plots import  SkewT\nfrom metpy.units import units\n\nimport geocat.datafiles as gdf\nfrom geocat.viz import util as gvutil"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Read in data:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Open a netCDF data file using xarray default engine and load the data into xarrays\nds = pd.read_csv(gdf.get('ascii_files/sounding_ATS.csv'),  header=None)\n\n# Extract the data\np = ds[0].values*units.hPa   # Pressure [mb/hPa]\ntc = ds[1].values*units.degC  # Temperature [C]\ntdc = ds[2].values*units.degC  # Dew pt temp  [C]\nwspd = ds[5].values*units.knots    # Wind speed   [knots or m/s]\nwdir = ds[6].values*units.degrees    # Meteorological wind dir\nu, v = mpcalc.wind_components(wspd, wdir)   # Calculate wind components"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig = plt.figure(figsize=(12,12))\n# Adding the \"rotation\" kwarg will over-ride the default MetPy rotation of \n# 30 degrees for the 45 degree default found in NCL Skew-T plots\nskew = SkewT(fig, rotation=45)\nax = skew.ax\n\n# Shade every other section between isotherms\nx1 = np.linspace(-100, 40, 8)  # The starting x values for the shaded regions\nx2 = np.linspace(-90, 50, 8)  # The ending x values for the shaded regions\ny = [1050, 100]  # The range of y values that the shaded regions should cover\n\nfor i in range(0, 8):\n    skew.shade_area(y=y,\n                    x1=x1[i],\n                    x2=x2[i],\n                    color='limegreen',\n                    alpha=0.25,\n                    zorder=1)\n\nskew.plot(p, tc, 'black')\nskew.plot(p, tdc, 'blue')\n# Plot only every third windbarb\nskew.plot_barbs(p[::3], u[::3], v[::3],\n                xloc=1.05,\n                fill_empty=True,\n                sizes=dict(emptybarb=0.075, width=0.1, height=0.2))\n\n# Draw line underneath wind barbs\nline = mlines.Line2D([1.05, 1.05], [0, 1],\n                     color='gray',\n                     linewidth=0.5,\n                     transform=ax.transAxes,\n                     dash_joinstyle='round',\n                     clip_on=False,\n                     zorder=0)\nax.add_line(line)\n\n# Add relevant special lines\n# Choose starting temperatures in Kelvin for the dry adiabats\nt0 = units.K * np.arange(243.15, 473.15, 10)\nskew.plot_dry_adiabats(t0=t0,\n                         linestyles='solid',\n                         colors='gray',\n                         linewidth=1.5)\n\n# Choose temperatures for moist adiabats\nt0 = units.K * np.arange(281.15, 306.15, 4)\nmsa = skew.plot_moist_adiabats(t0=t0, linestyles='solid', colors='lime', linewidths=1.5)\n\n# Choose mixing ratios\nw = np.array([0.001, 0.002, 0.003, 0.005, 0.008, 0.012, 0.020]).reshape(-1, 1)\n\n# Choose the range of pressures that the mixing ratio lines are drawn over\np_levs = units.hPa * np.linspace(1000, 400, 7)\nskew.plot_mixing_lines(w=w, p=p_levs, colors='lime')\n\nskew.ax.set_ylim(1000, 100)\n\ngvutil.set_titles_and_labels(ax, maintitle=\"ATS Rawinsonde: degC + Thin wind\")\n\n# Set axes limits and ticks\ngvutil.set_axes_limits_and_ticks(\n    ax=ax,\n    xlim=[-30, 50],\n    yticks=[1000, 850, 700, 500, 400, 300, 250, 200, 150, 100])\n\n# Change the style of the gridlines\nplt.grid(True,\n         which='major',\n         axis='both',\n         color='tan',\n         linewidth=1.5,\n         alpha=0.5)\nplt.xlabel(\"Temperature (C)\")\nplt.ylabel(\"P (hPa)\")\nplt.show()"
      ]
    }
  ],
  "metadata": {
    "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.9"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}