{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"top\"></a>\n",
    "<div style=\"width:1000 px\">\n",
    "\n",
    "<div style=\"float:right; width:98 px; height:98px;\">\n",
    "<img src=\"https://raw.githubusercontent.com/Unidata/MetPy/master/metpy/plots/_static/unidata_150x150.png\" alt=\"Unidata Logo\" style=\"height: 98px;\">\n",
    "</div>\n",
    "\n",
    "<h1>Downloading model fields using netCDF Subset Service (NCSS)</h1>\n",
    "<h3>Unidata Python Workshop</h3>\n",
    "\n",
    "<div style=\"clear:both\"></div>\n",
    "</div>\n",
    "\n",
    "<hr style=\"height:2px;\">\n",
    "\n",
    "<div style=\"float:right; width:250 px\"><img src=\"https://unidata.github.io/siphon/latest/_static/siphon_150x150.png\" alt=\"TDS\" style=\"height: 200px;\"></div>\n",
    "\n",
    "## Overview:\n",
    "\n",
    "* **Teaching:** 15 minutes\n",
    "* **Exercises:** 15 minutes\n",
    "\n",
    "### Questions\n",
    "1. What is the netCDF Subset Service (NCSS)?\n",
    "1. How can I use Siphon to make an NCSS request?\n",
    "1. How do I plot gridded fields using CartoPy?\n",
    "\n",
    "### Objectives\n",
    "1. <a href=\"#ncss\">Use siphon to make a request using NCSS</a>\n",
    "1. <a href=\"#projection\">Making sense of netCDF</a>\n",
    "1. <a href=\"#plotting\">Plot on a map</a>\n",
    "1. <a href=\"#pointrequest\">Requesting for a single point</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"ncss\"></a>\n",
    "## 1. What is NCSS?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Resolve the latest GFS dataset\n",
    "import metpy\n",
    "from siphon.catalog import TDSCatalog\n",
    "\n",
    "# Set up access via NCSS\n",
    "gfs_catalog = ('http://thredds.ucar.edu/thredds/catalog/grib/NCEP/GFS/'\n",
    "               'Global_0p5deg/catalog.xml?dataset=grib/NCEP/GFS/Global_0p5deg/Best')\n",
    "cat = TDSCatalog(gfs_catalog)\n",
    "ncss = cat.datasets[0].subset()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see what variables are available from ncss as well:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "ncss.variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From here, we can build a query to ask for the data we want from the server."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datetime import datetime, timedelta\n",
    "\n",
    "# Create a new NCSS query\n",
    "query = ncss.query()\n",
    "\n",
    "# Request data in netCDF format\n",
    "query.accept('netcdf')\n",
    "\n",
    "# Ask for our variable\n",
    "query.variables('Temperature_isobaric')\n",
    "\n",
    "# Ask for the 500 hPa surface\n",
    "query.vertical_level(50000)\n",
    "\n",
    "# Set the time range of data we want\n",
    "now = datetime.utcnow()\n",
    "query.time_range(now, now + timedelta(days=1))\n",
    "\n",
    "# Set the spatial limits\n",
    "query.lonlat_box(west=-110, east=-45, north=50, south=10)\n",
    "\n",
    "# get the data!\n",
    "data = ncss.get_data(query)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"#top\">Top</a>\n",
    "<hr style=\"height:2px;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"projection\"></a>\n",
    "## 2. Making sense of netCDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use a library called XArray to make working with this a little simpler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from xarray.backends import NetCDF4DataStore\n",
    "import xarray as xr\n",
    "\n",
    "# We need the datastore so that we can open the existing netcdf dataset we downloaded\n",
    "ds = xr.open_dataset(NetCDF4DataStore(data))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var = ds.metpy.parse_cf('Temperature_isobaric')\n",
    "var"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "XArray handles parsing things like date and times for us"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "latitude = var.metpy.y\n",
    "longitude = var.metpy.x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"#top\">Top</a>\n",
    "<hr style=\"height:2px;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"plotting\"></a>\n",
    "## Visualize the grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature as cfeature\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# GFS uses lon/lat grid\n",
    "data_projection = ccrs.PlateCarree()\n",
    "\n",
    "# Make it easy to change what time step we look at\n",
    "t_step = 0\n",
    "\n",
    "fig = plt.figure(figsize=(12, 8))\n",
    "ax = fig.add_subplot(1, 1, 1, projection=ccrs.LambertConformal())\n",
    "mesh = ax.pcolormesh(longitude, latitude, var[t_step].squeeze(),\n",
    "                     transform=data_projection, zorder=0)\n",
    "\n",
    "# add some common geographic features\n",
    "ax.coastlines(resolution='10m', color='black')\n",
    "ax.add_feature(cfeature.STATES, edgecolor='black')\n",
    "ax.add_feature(cfeature.BORDERS)\n",
    "\n",
    "# add some lat/lon gridlines\n",
    "ax.gridlines()\n",
    "\n",
    "# add a colorbar\n",
    "cax = fig.colorbar(mesh)\n",
    "cax.set_label(var.attrs['units'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <b>EXERCISE</b>:\n",
    "     <ul>\n",
    "        <li>Extend the plot above by plotting contours of 500 hPa geopotential heights</li>\n",
    "        <li>Add a title to the plot with the correct time</li>\n",
    "    </ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up an NCSS query from thredds using siphon\n",
    "query = ncss.query()\n",
    "\n",
    "#\n",
    "# SET UP QUERY\n",
    "#\n",
    "\n",
    "# Download data using NCSS\n",
    "#ncss.get_data(query)\n",
    "\n",
    "data_projection = ccrs.PlateCarree()\n",
    "\n",
    "# Plot using CartoPy and Matplotlib\n",
    "fig = plt.figure(figsize=(12, 8))\n",
    "ax = fig.add_subplot(1, 1, 1, projection=ccrs.LambertConformal())\n",
    "\n",
    "#\n",
    "# YOUR PLOT HERE\n",
    "#\n",
    "\n",
    "# add some common geographic features\n",
    "ax.coastlines(resolution='10m', color='black')\n",
    "ax.add_feature(cfeature.STATES, edgecolor='black')\n",
    "ax.add_feature(cfeature.BORDERS)\n",
    "\n",
    "# add some lat/lon gridlines\n",
    "ax.gridlines()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "scrolled": false
   },
   "source": [
    "<button data-toggle=\"collapse\" data-target=\"#sol1\" class='btn btn-primary'>View Solution</button>\n",
    "<div id=\"sol1\" class=\"collapse\">\n",
    "<code><pre>\n",
    "import numpy as np\n",
    "\n",
    "\\# Set up an NCSS query from thredds using siphon\n",
    "query = ncss.query()\n",
    "query.accept('netcdf4')\n",
    "query.variables('Temperature_isobaric', 'Geopotential_height_isobaric')\n",
    "query.vertical_level(50000)\n",
    "now = datetime.utcnow()\n",
    "query.time_range(now, now + timedelta(days=1))\n",
    "query.lonlat_box(west=-110, east=-45, north=50, south=10)\n",
    "\n",
    "\\# Download data using NCSS\n",
    "data = ncss.get_data(query)\n",
    "ds = xr.open_dataset(NetCDF4DataStore(data))\n",
    "\n",
    "temp_var = ds.metpy.parse_cf('Temperature_isobaric')\n",
    "height_var = ds.metpy.parse_cf('Geopotential_height_isobaric')\n",
    "longitude = temp_var.metpy.x\n",
    "latitude = temp_var.metpy.y\n",
    "time_index = 0\n",
    "\n",
    "\\# Plot using CartoPy and Matplotlib\n",
    "fig = plt.figure(figsize=(12, 8))\n",
    "ax = fig.add_subplot(1, 1, 1, projection=ccrs.LambertConformal())\n",
    "\n",
    "contours = np.arange(5000, 6000, 80)\n",
    "ax.pcolormesh(longitude, latitude, temp_var[time_index].squeeze(),\n",
    "              transform=data_projection, zorder=0)\n",
    "ax.contour(longitude, latitude, height_var[time_index].squeeze(), contours, colors='k',\n",
    "           transform=data_projection, linewidths=2, zorder=1)\n",
    "ax.set_title(temp_var.metpy.time[time_index].values)\n",
    "\n",
    "\\# add some common geographic features\n",
    "ax.coastlines(resolution='10m', color='black')\n",
    "ax.add_feature(cfeature.STATES, edgecolor='black')\n",
    "ax.add_feature(cfeature.BORDERS)\n",
    "\n",
    "\\# add some lat/lon gridlines\n",
    "ax.gridlines()\n",
    "</pre></code>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"#top\">Top</a>\n",
    "<hr style=\"height:2px;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"pointrequest\"></a>\n",
    "## 4. NCSS Point Request\n",
    "We can also request data for a specfic lon/lat point, across vertical coordinates or times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "cat = TDSCatalog('http://thredds.ucar.edu/thredds/catalog/grib/NCEP/GFS/'\n",
    "                 'Global_0p5deg/catalog.xml?dataset=grib/NCEP/GFS/Global_0p5deg/Best')\n",
    "ncss = cat.datasets[0].subset()\n",
    "\n",
    "point_query = ncss.query()\n",
    "point_query.time(datetime.utcnow())\n",
    "point_query.accept('netcdf4')\n",
    "point_query.variables('Temperature_isobaric', 'Relative_humidity_isobaric')\n",
    "point_query.variables('u-component_of_wind_isobaric', 'v-component_of_wind_isobaric')\n",
    "point_query.lonlat_point(-101.877, 33.583)\n",
    "\n",
    "# get the data! Unfortunately, xarray does not quite like what comes out of thredds\n",
    "point_data = ncss.get_data(point_query)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Skew-T diagrams typical use specific units. First, let's assign units to the variables we requested:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from metpy.units import units\n",
    "\n",
    "# get netCDF variables\n",
    "pressure = point_data.variables[\"isobaric\"]\n",
    "temp = point_data.variables[\"Temperature_isobaric\"]\n",
    "u_cmp = point_data.variables[\"u-component_of_wind_isobaric\"]\n",
    "v_cmp = point_data.variables[\"v-component_of_wind_isobaric\"]\n",
    "relh = point_data.variables[\"Relative_humidity_isobaric\"]\n",
    "\n",
    "# download data and assign the units based on the variables metadata\n",
    "# Need to put units on the left to assure things work properly with masked arrays\n",
    "p = units(pressure.units) * pressure[:].squeeze()\n",
    "T = units(temp.units) * temp[:].squeeze()\n",
    "u = units(u_cmp.units) * u_cmp[:].squeeze()\n",
    "v = units(v_cmp.units) * v_cmp[:].squeeze()\n",
    "relh = units('percent') * relh[:].squeeze()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also need to calculate dewpoint from our relative humidity data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import metpy.calc as mpcalc\n",
    "\n",
    "Td = mpcalc.dewpoint_rh(T, relh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's change those units to what we typically see used in Skew-T diagrams. We use `ito` to do this in-place rather than manually reassigning to the same variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p.ito(units.millibar)\n",
    "T.ito(units.degC)\n",
    "Td.ito(units.degC)\n",
    "u.ito(units.knot)\n",
    "v.ito(units.knot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from metpy.plots import SkewT\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=45)\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, 'tab:red')\n",
    "skew.plot(p, Td, 'blue')\n",
    "skew.plot_barbs(p, u, v)\n",
    "skew.ax.set_ylim(1000, 100)\n",
    "skew.ax.set_xlim(-40, 60)\n",
    "\n",
    "# Add the relevant special lines\n",
    "skew.plot_dry_adiabats()\n",
    "skew.plot_moist_adiabats()\n",
    "skew.plot_mixing_lines()"
   ]
  }
 ],
 "metadata": {
  "gist_id": "310b790a3a2cfdc8d06e",
  "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}