{
 "metadata": {
  "name": "",
  "signature": "sha256:fa99317ee94ae5086634916914e7293636eb8df9fd5ef4cbd7ca0e80f68d7871"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Using Python to get the latest NEXRAD Composite\n",
      "Objective: Visualize the latest available reflectivity data composited data\n",
      "    \n",
      "Steps involved:\n",
      "\n",
      "- Define a color map based on the one used by the National Weather Service\n",
      "- Construct the appropriate URL to get the latest data file\n",
      "- Open the URL using netCDF4-python\n",
      "- Read the basic metadata\n",
      "- Create the appropriate CartoPy projection and plot the Radar Reflectivity"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Set-up for notebook\n",
      "%matplotlib inline\n",
      "\n",
      "# Some needed imports\n",
      "import datetime as dt\n",
      "import matplotlib.pyplot as plt\n",
      "import matplotlib as mpl\n",
      "import cartopy\n",
      "import numpy as np\n",
      "from netCDF4 import Dataset\n",
      "from pyudl.tds import TDSCatalog"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Define the NWS Color Map for reflectivity\n",
      "- While matplotlib has some really nice, built-in colormaps, nothing quite matches the NWS colormap\n",
      "- Python List of HTML colors (used photoshop color picker on a screenshot of the colorbar off of the NWS radar page)\n",
      "- Converted into matplotlib colormap using ListedColormap"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Function to create a colortable that matches the NWS colortable\n",
      "def radar_colormap():\n",
      "    nws_reflectivity_colors = [\n",
      "    \"#646464\", # ND\n",
      "    \"#ccffff\", # -30\n",
      "    \"#cc99cc\", # -25\n",
      "    \"#996699\", # -20\n",
      "    \"#663366\", # -15\n",
      "    \"#cccc99\", # -10\n",
      "    \"#999966\", # -5\n",
      "    \"#646464\", # 0\n",
      "    \"#04e9e7\", # 5\n",
      "    \"#019ff4\", # 10\n",
      "    \"#0300f4\", # 15\n",
      "    \"#02fd02\", # 20\n",
      "    \"#01c501\", # 25\n",
      "    \"#008e00\", # 30\n",
      "    \"#fdf802\", # 35\n",
      "    \"#e5bc00\", # 40\n",
      "    \"#fd9500\", # 45\n",
      "    \"#fd0000\", # 50\n",
      "    \"#d40000\", # 55\n",
      "    \"#bc0000\", # 60\n",
      "    \"#f800fd\", # 65\n",
      "    \"#9854c6\", # 70\n",
      "    \"#fdfdfd\" # 75\n",
      "    ]\n",
      "\n",
      "    return mpl.colors.ListedColormap(nws_reflectivity_colors)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Get the latest data URL, grab the metadata, and request the data"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Get today's date...\n",
      "today = dt.datetime.utcnow()\n",
      "\n",
      "# ...and use that to assemble the URL and grab the catalog\n",
      "url = \"http://thredds.ucar.edu/thredds/catalog/nexrad/composite/gini/n0r/1km/{}/catalog.xml\".format(today.strftime(\"%Y%m%d\"))\n",
      "cat = TDSCatalog(url)\n",
      "\n",
      "# Get the list of names of datasets\n",
      "names = cat.datasets.keys()\n",
      "print names"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# sort, so that the last dataset is the latest\n",
      "names.sort()\n",
      "latest = names[-1]\n",
      "print latest"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Grab the dataset for the latest\n",
      "latestDs = cat.datasets[latest]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Construct a NetCDF dataset using the OPeNDAP access URL\n",
      "dataset = Dataset(latestDs.accessUrls['OPENDAP'])\n",
      "print dataset.variables.keys()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dataset.variables['LambertConformal'].ncattrs()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "##################\n",
      "# Projection fun #\n",
      "##################\n",
      "\n",
      "# get basic info from the file regarding projection attributes\n",
      "# see the following for more info on projections as described here:\n",
      "#   http://www.nco.ncep.noaa.gov/pmb/docs/on388/tableb.html#GRID218\n",
      "#   http://www.wmo.int/pages/prog/www/WDM/Guides/Guide-binary-2.html [see LAMBERT CONFORMAL SECANT OR TANGENT CONE GRIDS]\n",
      "#   http://www.unidata.ucar.edu/mailing_lists/archives/netcdf-java/2006/msg00006.html [starndard parallels in CDM]\n",
      "proj = dataset.variables['LambertConformal']\n",
      "rsphere = proj.earth_radius\n",
      "\n",
      "# lat_0\t: center of desired map domain (in degrees) [Basemap]\n",
      "# CDM : 'latitude_of_projection_origin'\n",
      "lat_0 = proj.latitude_of_projection_origin\n",
      "\n",
      "# lon_0\t: center of desired map domain (in degrees) [Basemap]\n",
      "# CDM : 'longitude_of_central_meridian'\n",
      "lon_0 = proj.longitude_of_central_meridian\n",
      "\n",
      "# lat_1, lat_2 : 1st and second parallels [Basemap]\n",
      "# CDM : 'standard_parallel' - this attr contains both 1st and 2nd\n",
      "lat_1 = proj.standard_parallel\n",
      "print lat_0, lon_0, lat_1, rsphere"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Grab the data"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Used to subset the data\n",
      "xstride = 10\n",
      "ystride = 10\n",
      "\n",
      "# download x and y coords and convert them from km to m\n",
      "x = dataset.variables['x'][::xstride] * 1000.\n",
      "y = dataset.variables['y'][::ystride] * 1000.\n",
      "\n",
      "# Grab the reflectivity data. Mask values less than -30 dBz\n",
      "data = dataset.variables[\"Reflectivity\"][0, 0::ystride, 0::xstride]\n",
      "data = np.ma.array(data, mask=data<=-30)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Create the Plot"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Set up the projection for the LambertConformal projection we know we have\n",
      "lcc = cartopy.crs.LambertConformal(central_longitude=lon_0, central_latitude=lat_0, secant_latitudes=(lat_0, lat_1))\n",
      "\n",
      "# Create a large figure and axes with this projection\n",
      "fig = plt.figure(figsize=(24, 12))\n",
      "ax = fig.add_subplot(1, 1, 1, projection=lcc)\n",
      "\n",
      "# Limit to the bounds of the data we have\n",
      "ax.set_extent([-129., -63., 22., 49.], cartopy.crs.Geodetic())\n",
      "\n",
      "# Add some map features\n",
      "ax.stock_img()\n",
      "ax.coastlines(resolution='50m')\n",
      "ax.add_feature(cartopy.feature.NaturalEarthFeature(category='cultural', name='admin_1_states_provinces_lines', scale='50m',\n",
      "                                                   facecolor='none'))\n",
      "ax.add_feature(cartopy.feature.BORDERS, linewidth='2', edgecolor='black')\n",
      "ax.gridlines()\n",
      "\n",
      "# Convert the time to text and add as the title\n",
      "time = dataset.variables[\"time\"][:][0] / 1000.\n",
      "title = dt.datetime.fromtimestamp(time).isoformat()\n",
      "ax.set_title(title)\n",
      "\n",
      "# Plot the data as an image, using the x and y values we have as the extents of the image\n",
      "# NOTE: This assumes equal-spaced points\n",
      "cmap = radar_colormap()\n",
      "norm = mpl.colors.Normalize(vmin=-35, vmax=80)\n",
      "cax = ax.imshow(data, extent=(x.min(), x.max(), y.min(), y.max()), cmap=cmap, norm=norm, origin=\"upper\",\n",
      "                transform=lcc)\n",
      "plt.colorbar(cax)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Exercise"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Using what was done above, plot the digital hybrid reflectivity (DHR):\n",
      "- Look at http://thredds.ucar.edu/thredds/catalog/nexrad/composite/gini/catalog.html\n",
      "- Instead of plotting over all of the U.S., limit to an area of interest\n",
      "- DHR was chosen to keep the colormap from the NWS the same. Can also look at:\n",
      "    - Echo Tops (EET)\n",
      "    - Digital VIL (DVL)\n",
      "    - Others in catalog\n",
      "- ***Bonus points***: Plot the data into a new coordinate system, like Orthographic"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}