{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# NCL_skewt_2_2.py\nThis script illustrates the following concepts:\n - Customizing the background of a Skew-T plot\n - Plotting temperature, dewpoint, and wind data on 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_2.ncl\n - Original NCL plots: https://www.ncl.ucar.edu/Applications/Images/skewt_2_2_lg.png\n\nNote:\n Currently functions to calculate CAPE, precipitable water, the showalter\n index, the pressure of the LCL, and the temperature of the LCL do not\n exist in ``geocat-comp``. An `issue `_\n has been opened on the ``geocat-comp`` GitHub. The subtitle with those\n values will be added at a later date once that issue has been closed.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Import packages:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\nimport matplotlib.lines as mlines\nimport numpy as np\nimport pandas as pd\nfrom metpy.plots import SkewT\nfrom metpy.units import units\nimport metpy.calc as mpcalc\n\nimport geocat.viz.util as gvutil\nimport geocat.datafiles as gdf" ] }, { "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.testdata'), delimiter='\\\\s+', header=None)\n\n# Extract the data\np = ds[1].values*units.hPa # Pressure [mb/hPa]\ntc = (ds[5].values + 2)*units.degC # Temperature [C]\ntdc = ds[9].values*units.degC # Dew pt temp [C]\n\n# Create dummy wind data\nwspd = np.linspace(0, 150, len(p))*units.knots # Wind speed [knots or m/s]\nwdir = np.linspace(0, 360, len(p))*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": [ "# Note that MetPy forces the x axis scale to be in Celsius and the y axis\n# scale to be in hectoPascals. Once data is plotted, then the axes labels are\n# automatically added\nfig = plt.figure(figsize=(12, 12))\n\n# The rotation keyword changes how skewed the temperature lines are. MetPy has\n# a default skew of 30 degrees\nskew = SkewT(fig, rotation=45)\nax = skew.ax\n\n# Plot temperature and dew point\nskew.plot(p, tc, color='black')\nskew.plot(p, tdc, color='blue')\n\n# Draw parcel path\nparcel_prof = mpcalc.parcel_profile(p, tc[0], tdc[0]).to('degC')\nskew.plot(p, parcel_prof, color='red', linestyle='--')\nu = np.where(p>=100*units.hPa, u, np.nan)\nv = np.where(p>=100*units.hPa, v, np.nan)\np = np.where(p>=100*units.hPa, p, np.nan)\n\n# Add wind barbs\nskew.plot_barbs(p=p[::2],\n u=u[::2],\n v=v[::2],\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 clip_on=False,\n zorder=1)\nax.add_line(line)\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 shades regions should cover\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\n# Choose starting temperatures in Kelvin for the dry adiabats\nt0 = units.K * np.arange(243.15, 444.15, 10)\nskew.plot_dry_adiabats(t0=t0, linestyles='solid', colors='tan', linewidths=1.5)\n\n# Choose starting temperatures in Kelvin for the moist adiabats\nt0 = units.K * np.arange(281.15, 306.15, 4)\nskew.plot_moist_adiabats(t0=t0,\n linestyles='solid',\n colors='lime',\n linewidth=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)\n\n# Plot mixing ratio lines\nskew.plot_mixing_lines(w=w,\n p=p_levs,\n linestyle='dashed',\n colors='lime',\n linewidths=1)\n\n# Use geocat.viz utility functions to set axes limits and ticks\ngvutil.set_axes_limits_and_ticks(\n ax=ax,\n xlim=[-32, 38],\n yticks=[1000, 850, 700, 500, 400, 300, 250, 200, 150, 100])\n\n# Use geocat.viz utility function to change the look of ticks and ticklabels\ngvutil.add_major_minor_ticks(ax=ax,\n x_minor_per_major=1,\n y_minor_per_major=1,\n labelsize=14)\n# The utility function draws tickmarks all around the plot. We only need ticks\n# on the left and bottom edges\nax.tick_params('both', which='both', top=False, right=False)\n\n# Use geocat.viz utility functions to add a main title\ngvutil.set_titles_and_labels(ax=ax,\n maintitle=\"Raob; [Wind Reports]\",\n maintitlefontsize=22,\n xlabel='Temperature (C)',\n ylabel='P (hPa)',\n labelfontsize=14)\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)\n\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 }