python-awips/examples/notebooks/NEXRAD_Level_3_Plot_with_Matplotlib.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Shown here are plots for Base Reflectivity (N0Q, 94) and Base Velocity (N0U, 99) using AWIPS data rendered with Matplotlib, Cartopy, and MetPy.  This example improves upon existing Level 3 Python rendering by doing the following:\n",
    "\n",
    "* Display scaled and labeled colorbar below each figure.\n",
    "* Plot radar radial images as coordinate maps in Cartopy and label with lat/lon.\n",
    "* Colormap and data scaling copied from operational AWIPS.\n",
    "* Level 3 data are retrieved from the [Unidata EDEX Cloud server](http://unidata.github.io/awips2/docs/install/install-cave.html#how-to-run-cave) (`edex-cloud.unidata.ucar.edu`)\n",
    "* Raw HDF5 byte data are converted to product values and scaled according to (page 3-34 https://www.roc.noaa.gov/wsr88d/PublicDocs/ICDS/2620001U.pdf)\n",
    "\n",
    "        The threshold level fields are used to describe (up to) 256 levels as follows:\n",
    "         halfword 31 contains the minimum data value in m/s*10 (or dBZ*10)\n",
    "         halfword 32 contains the increment in m/s*10 (or dBZ*10)\n",
    "         halfword 33 contains the number of levels (0 - 255) \n",
    "\n",
    "According to the [ICD for the Product Specification](https://www.roc.noaa.gov/WSR88D/PublicDocs/NewTechnology/B17_2620003W_draft.pdf), *\"the 256 data levels of the digital product cover a range of reflectivity between -32.0 to +94.5 dBZ, in increments of 0.5 dBZ. Level codes 0 and 1 correspond to 'Below Threshold' and 'Range Folded', respectively, while level codes 2 through 255 correspond to the reflectivity data itself\"*.\n",
    "\n",
    "So it's really 254 color values between -32 and +94.5 dBZ.\n",
    "\n",
    "The ICD lists 16 specific color levels and directs 256-level reflectivity products to use corresponding colors, leaving it the rendering application to scale and blend between the 16 color values, and to make decisions about discrete color changes, apparently.\n",
    "![](http://i.imgur.com/cqphoe3.png)\n",
    "\n",
    "For AWIPS, the National Weather Service uses a mostly-blended color scale with a discrete jump to red at reflectivity values of 50 dBZ:\n",
    " \n",
    "![](http://i.imgur.com/o18gmio.png)\n",
    "\n",
    "50 dBZ corresponds to the 16-level color *light red* (**FF6060**). Note that `FF6060` is not used in the NWS AWIPS color scale, instead RGB value is given as `255,0,0` (or as percent `1.0,0.0,0.0`). 60 is not quite exactly where white starts, but it makes convenient sense enough I don't see a reason why not. Obviously the AWIPS D2D authors took some liberties with their 256-level rendering, not adhering to \"dark red\" for dBZ values between 60-65 (white was for 70 dBZ and above on the 16-level colormap).  Still the goal is to follow operational AWIPS as closely as possible, and that means recreating the rendering of their DHR/N0Q plots to the colored pixel, so we will assume 50 dBZ should be red and 60 dBZ white, and 75 dBZ cyan."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtIAAAJ7CAYAAADOeXaMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXtcVGX+x99nhuFqCWooqRla3qpNsaJct5S2tUz3t4st\ngmRt2qqpuFvUmpVWZqtdbFu1MlNrUxDYsotmua1KZSpbom7mrZRMjSQVTIGRYc7z++M558wMDHIR\nufm8Xy9ezJw5l2ceBvjMdz7P56sJIVAoFAqFQqFQKBS1w9bYA1AoFAqFQqFQKJojSkgrFAqFQqFQ\nKBR1QAlphUKhUCgUCoWiDighrVAoFAqFQqFQ1AElpBUKhUKhUCgUijqghLRCoVAoFAqFQlEHAhp7\nAHVB0zSV2adQKBQKhUKhaBCEEJq/7c1SSAOo/GvIzs5m4MCBjT2MeiU/P5/evXvz/fffc8EFF9Tp\nHC1xXuoDNS/+UfNSNWpuICFwYaVtBfoeIm09qj02q2ys3/P5236m65nnqviYeZ5TWrDP9lbCWe3Y\nzgXq9VI1am7801zmRdP8amigGQtpRcskKiqKm266iaysLMaMGdPYw1EoFOc5plitSuDWhZqcq6Jw\n9neMuS2rkYSzQqEArTlWdjVNE81x3IqasWrVKp5++mk2bdrU2ENRKBQKoPZCuqpq9Jn2r801zlSh\nrsh+LYiu4nSNz61QKHzRNK1Ka4cS0oomR3l5OZdeeilr1qzhiiuuaOzhKBQKBfxnrXUzYci+anev\nStT6E8u1Fd3+jq9KVO/Xgny2Vyeovc9zJguKQnE+cSYhrVI7mjHZ2dmNPYRzQkBAAIMGDWLLli11\nOr6lzsvZoubFP2peqkbNjQf9o0Lra8Ksk2St7lan82SVjbW+TBICF/p8ASxxTa7xOauyfSQELqSr\nOG19gRTW5ldV46t4jppS8fVS2+NbMup3yT8tYV6UR1rRJGnfvj1Hjhxp7GEoFAoFO8f+EoDeCz+X\nG7KzYUskICvTWW/3JGH4bmv/s63kJgQuBMdc61xHtCDC0BhtbKsrNbV3VPSFW17sOj4vVeVWtGSU\ntUPRJHnuuec4cuQIzz//fGMPRaFQnOeYQtqb3j1eA0CUFgIwYsbXVR5fG5sHQPpgKcpHrunps/0F\nVwoPOOb5bEsbsJnkDddXee2ajKMhUGJa0ZxRHmlFs+PNN9/k3//+N8uWLWvsoSgUCkVlVn4KwO7f\n/pruD6/A1q4rCVM3VLl7bXzQGZPboOdthbJSz8bAEEau7GLd9RbVVYnpqhJHGltQKyGtaG4oj3QL\npSV4i6ribKwdLXlezgY1L/5R81I1am78k/3a36BtCbQtoefn72PrfyF0P0rW2z2rP9gLf4Iyc2Y/\nCItg5MoujFzTE1uP/jhXvewjotOHHbBE9DTXJB8RnTZgs3Xb25JR0Zd9LqjJ6+V8FdHqd8k/LWFe\nlEda0SRRHmmFQtFk+TkIDraG7keh1CFF9aZL4Op8zyLEX98MVF+Brbi4r2jKL3EAaQPiCRiWiigu\nIjR1GUvm3ElIyiKSFpT4iOqnHPOt21JUz/c+fSVLRW2E7Nl6oxWK8wFl7VA0SfLz8+nTp48S0wqF\noumx8lMocUCoS94PccH2KHn76nxrt7W33ArAzcJVq9N7i99XXZNxACGjnsJ+Vwz6R4WUGKJaCwsn\ncc4PdXoKNRHHjW0JUTYQRVNBeaQVzY7y8nJCQkJwOp3Y7fbGHo5CoVBIzDzptiVwqLVnuymsQ1yy\nSr23HWsfvNrvKWojrP01U3Hf9Zz0UHthHzCyTqK6qQrqxhbxCoU359wjrWlasKZpOZqmbdM0baem\nabOM7VdrmrZJ07T/aZr2vqZpF3gds8TY/3bj/qWapumapk3y2me+pml318cYWyItwVtUFQEBAYSH\nh3P06NFaH9uS5+VsUPPiHzUvVaPmxg+hLrK3b5MiOsQQzsdC5dfB1nL73nbQ6QSvOl6yvhxnccmK\n2c9JGa3B5SQ5ZyDJOQMBcG9IJy02u0aCs7aitKaWkPp8vVSVtd1cUb9L/mkJ81IvHmkhhFPTtEFC\niBJN0wKADZqmDQD+DjwghPhM07R7gIeA6ZqmXQl8D/wJSAc+ME5VAEzWNO1VIYQLUGXn85gOHTpw\n5MgR2rdv39hDUSgUjcHu+ZW3RYRX3tb+znM/FpCWjhAHBLllRdob877pmT4W6vPwjRmfkjDqKwBu\nrsUlvavRppieDRzMARzSOoLLSXKuvJ1FZd91ReoqSKuL8SvQ95BdPrBO567umt6Z1qo6rWhK1Lu1\nQ9O0UOAT4I/A50KIcGN7Z+AjIcQVmqb1BMYAjwOvCyFGaJp2KbAS2ABsEUIs0jRtHvClEOKfFa6h\nrB3nAb/+9a+ZMmUKt9xyS2MPRaFQNAbfLvC9X14OF7Ty3A/wqgXVREzvWeF7v0d87cZjRN75cCzU\nI5xNW4c/QlyWkD4bIei63jPmgzkfWLcfdsyr8XnP1Ka8rkK1pq3PK+6rRLGiOXAma0e9pXZommYD\ncoFuwCtCiK81Tfta07T/E0K8B/wB6AwghNhtVK4/AVIrnOpZ4ENN05bU19gUzROV3KFQnOdcNr7y\ntu8zPWI1+Dg420gB+32m736XjKh8rCmcTUFtfq+poC4xrtu2xBqDOHQajVBPRTrUhdgVhhYpfIW1\nl8CuriKcVTa2SrF5Omc1QbFDAOgce7tnn801F6QVz19xTHURulVlVleHSgZRNHfqLUdaCKELIfoA\nnYAbNU0bCIwGJmia9iXQCijz2v9+IcS1QohPK5wnD8gBRtbX2FoqLcFbdCbat2/Pjz/+WOvjWvq8\n1BU1L/5R81I1TXJutkdJD3KJA463l2LaXxX4yBmaOfWI9/2qCXN2Sv/zodZkv31QjmFvO7Sr3XIs\nhk9abLej2UMQBRr6Ti/7R4nDJ+P5TPgTuebXaMdcknNvxT54IqdzVtds7H6oKlu64v3aeJOzs7Nr\nlFnt7/Hm7H+uCU3yd6kJ0BLmpd5zpIUQJzRN+wC4RgjxPDAYQNO07sDtZzzYw9+At5AVa0UVbNu2\nDYCBAwcCnhdkS7l/6tQpvvjiC+v51vT42u5/vtxv6a+Xut43aSrjaUr3t23b1qTGAzCwRC64y/5q\nKwSVM7BPHwgvIfvzXXAqiIGxV3ru8ygDf9lL7r9Mfro1MDW1btf/9r/y/mXX+d7nOuh0guxNskX4\nwJu7o2/cz6elPxqP94FjoWT/9Bni8liy7ogl4S85TJhg44n5u4i09QCkvxio8f24p1eDfSKRuT3I\nKhtb4+fz8m/2Wud74j/DGDhwoM/xgM/9l3+zt8bnr+3rZcK/uzNw4EASAhdKf3V2duO/vtTf3wa9\n3xKoF4+0pmntgHIhRJGmaSHAGuBJ4H9CiJ8M28cbwDohxBtVnONSYKUQ4irjfiZwPTBNCPFmhX2V\nR/o84I033mDdunW8+eab1e+sUChaPhU9yiXGwr7OJ3y3RR7z9U8DfBXl/5y/rs3SPy9e3u65bVbE\nO52otMhQP3gQW+fOEOpCzzuBFhaOFikQBRojHttS5ekrLhhMi9tG8ro+1e5bE1S0nEJROxrCIx0F\n/NMQzDZgqRBiraZpf9Y0bYKxz9tViWgvvNXx08DWqnZUtHyUR1qhUPgw7MbK2zI3SfFqium2JbBB\nVm2tzoMgM567VYjT/LAnvPEl/PGa2o8l1OXxTHc/6rF3mGMwbpsimhIHWlg4orgIraQd+p5PgSC/\np/YnbL1FtLeorosIVkkYCkX9YauPkwghvhJCxAgh+gghfiGEeM7Y/g8hRA/j65FqzvGdEOIXXvf/\nJ4SwV6xGKzy0pI9G/FFXId3S56WuqHnxj5qXqmkWczNgjxSqZo7ztiiPqN7bTrbuLnHIajHIyrT5\nZfLGl75f1WDNS6ixmPBYqPxudjqs6Nk2Fh1q9hDrOFt0DMtHnWb5mOr/DVcUuaaITovbhuvWCf4O\nqRE18TPXhsZ+vTTlrOn
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10d493f90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtIAAAJ7CAYAAADOeXaMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXt8VNW99r97MrekM5MMJGBM8JIANuAFrYDHtngpqLzi\nDeu1tl7aarVvW209R61Va1urVvvWczxtT2tPFW+oWOu1iKKlVq2AIlUBFRIrEANJIGRmTGYyk9nv\nH3vW2mvvmUCA3LO+n8/+MLNnZu81K0Py7N886/kZpmmi0Wg0Go1Go9Fodg/PYA9Ao9FoNBqNRqMZ\njmghrdFoNBqNRqPR7AFaSGs0Go1Go9FoNHuAFtIajUaj0Wg0Gs0eoIW0RqPRaDQajUazB2ghrdFo\nNBqNRqPR7AHewR7AnmAYhs7s02g0Go1Go9EMCKZpGoX2D0shDaDzr2HZsmUce+yxgz2MPqWpqYkp\nU6awceNGwuHwHh1jJM5LX6DnpTB6XnpmtM7N3CsXFdwfnTYBgK3vr2L8Z4+Q+x++6KgBGddQZ7R+\nXnqDnpvCDJd5MYyCGhrQ1g7NEKOyspJjjjmGxx57bLCHotFoRimL7zorb58nMGzrThqNph8xhmNl\n1zAMcziOW9M7nn32WW655Rb+8Y9/DPZQNBqNBoBTb3yG7lQmb39pXWW/VaTP+f2rADx66Rf65fga\njaZ3GIbRo7VDC2nNkCOTyXDAAQewZMkSpk6dOtjD0Wg0Gk698Zm8fd2pDM/dfka/nE+I6EJoYa3R\nDCw7E9La2jGMWbZs2WAPoV/wer0cd9xxvPXWW3v0+pE6L3uLnpfC6HnpGT03FvPvXEqopkJun3Zt\nIlRTQaAizPw7lzq2vuLRS7+QJ5g7N7fRubmNU298Rm5DCf156Rk9N4UZCfOiTV+aIcn48ePZunXr\nYA9Do9FoeOLq2Q6R3NXeSbI1UfC58+9cyhNXz+6zc6ti+oKFbxJb1+R4/NQbn8EbKZbj1Gg0A4u2\ndmiGJHfccQdbt27lzjvvHOyhaDSaUc7uVpr3VtCefM2f5W1hHTnr7mUEykOO58XWNbHPnCm0Lv+o\nz8eg0WhstLVDM+wYP348W7ZsGexhaDQaTUF8kaBj6ytUES3un3zNn+lOZUi1JuQGEKmrBKB85oF5\nx+lLm4lGo+kZLaSHMSPBW9QTe2PtGMnzsjfoeSmMnpee0XNj8cTVsx1bl38rRX6vYxOCem8rwUUB\nr9wEJdVRwFrcKLZUawJfpJh0LAlYYlps2VSabCrN6bcs3qux7C7689Izem4KMxLmRXukNUMS7ZHW\naDRDmVBthbydqG+Rt/cmss69eFCI6VRLnECF3aBKiOnmVz4EYNysyfIxd3VcFdNPXj93t8ek0Wh2\njvZIa4YkTU1NTJs2TYtpjUYz5PjWio8d97OpDO2uRYAqfSGqn/7JKfL2WXcvo2NzW8HXie6LAImG\nFsdjWkhrNHuGzpHWDDsymQzFxcUkk0mKiooGezgajUYDWMkZbkK1FY6qdDpu2S1ia21xXahbYl/g\n9lSrlXJf2FmdfvC8I/tlDBrNSKffFxsahhE0DGO5YRirDcNYaxjGrbn9hxmG8Q/DMN4xDONpwzDC\nymv+mHv+ybn7BxiGkTUM4/8qz/lvwzAu7IsxjkRGgreoJ7xeL2VlZbS2tu72a0fyvOwNel4Ko+el\nZ/TcFGbrWmfGfTqWJFARlpsvHMwTsXOvXCS3vuS528+QGzjzp9PxpNwGQkTrz0vP6LkpzEiYlz7x\nSJummTQM4zjTNDsMw/ACrxqG8QXgV8D3TdP8u2EYFwP/DtxoGMbBwEbgm8DDwHO5QzUD3zUM43em\naaYBXXYexeyzzz5s3bqV8ePHD/ZQNBrNIHDC5Qvz9lXOqQNgwfxpAz0cTr3xGcbOrKGoxI8vl92c\njnUSrLBj6ZItCQIVYdpWb6K4OkpnDxaM/kDtsqi7H2o0A0OfLTY0TbMjd9MPFAFtwCTTNP+e278U\neB64EcgAnwECrsO0AK8CFwJ/6KuxjVSOPfbYwR5Cv7KnCw5H+rzsKXpeCqPnpWeG8txc+MRqx/2B\nEtbbljfgp5xtyxsYO7OGUG0F7YqFo3SKFUknF/3lEjcgf9HgnnDW3cvk7UXfOXavjtXXDOXPy2Cj\n56YwI2Fe+kxIG4bhAVYBtcBvTdNcYxjGGsMwTjNN8yngLGACgGma7+cq138DfuA61C+AxYZh/LGv\nxqYZnujkDo1mdKNGwAHsf/4MgLwqryfg4+Ln1gBw78lTezxeoZba6iK+XVGsiGKwRPW+8w6VvuRE\nfYtDVIMzRSPVEifVEi/YcGVPGMqiWqMZLfRZjrRpmlnTNKcB1cAswzCOBS4BrjAM400gBHQpz7/K\nNM3ppmm+4jrOR8By4Py+GttIZSR4i3bGnjZlGenzsqfoeSmMnpeeGey5WXzXWY5NUFwddWwqFz+3\nhoufW8MJly90bD1RSFwX4vRbFkvv8/ZNa+Q4VEK1FVJUuz3SRf7CdSv3YsFdseg7xw5Z0TzYn5eh\njJ6bwoyEeenzHGnTNNsNw3gOONI0zTuBEwEMw5gMnNzLw/wceByrYq3pgdWrra82xVcj4gM5Uu4n\nEglWrlwp329vX7+7zx8t90f652VP7wuGyniG0v3Vq1cPqfH86x8fccCXjgfgk7f+AcA+B0/HFwny\nydtvkPm0i8pDZ5BqibO9cS0AY6qmAPDFy38JQEXtYQC01P+T7k9TlB9wSK/O39rwDq0NUF5zKJ2f\n7GA7a/CVFROKWcK5Zf3b1vEnHU46nqT5g7fxlvgZP+Vz1uv/9S5po5Pxnz0CgK3vryLVEqe85lBO\nv2UxrQ3v8LOvzuzx/Eed91MAxu53MAD/Pjc64PO/q/tD7fMylO7r37+F748E+iT+zjCMciBjmuYO\nwzCKgSXAzcA7pmm25Gwf9wEvm6Z5Xw/HOAB4xjTNQ3L3HwWOAm4wTfN+13N1/N0o4L777uPll1/m\n/vvv3/WTNRrNiEfYN1RKqssAZIc/cDZIAWhd/pEjFk6lt50I3Z7s9nVNhGrsY7q9z2qGc6imgmwq\n43g81RKno7HwQsRCec89Va6zqUy/RetpNBqLncXf9VVFuhJYkBPMHuAB0zRfMgzje4ZhXJF7zp96\nEtEKqjq+BXi7j8anGYZoj7RGoxHMvXIRnkD+n6wJZ1oVXtWLrFouOje3UT7zQNpWb3K8zm3L2Bmq\nFxmgpKqM0rpKxzGEeBeCWhXZ2VQmb+yegNchrgu9NxVVqKda4vK4gCNST4tqjWZg8fTFQUzTfNc0\nzSNM05xmmuahpmnekdv/n6ZpHpTbfriLY/zLNM1DlfvvmKZZ5K5Ga2xG0lcjhdhTIT3S52VP0fNS\nGD0vPTPU5iabyuRtohIt2mZ3uyq/bj+12Do3t/U6mk74ksX2r9y8qMcI1VZQUh3NWyAJVkSee9wA\nkbpKInWVee+tEE9cPVtWz0VedaGLgb7Oqd4dhtrnZSih56YwI2Fe+twjrdH0FboirdFoBIUqrRcs\nfFNWZwXBihBFuVznZEsCgJZXPnQ8xxsplhVsdbFhoQQPdzV60XeO5eb503i4yK4QCzEtRL1qIxHC\nPh3rlPs6Gnc4jinEdPMr6/PO76YnK8pgCmiNZjSjW4RrhizpdJqSkhJSqRQeT598eaLRaEYQ59/3\nRp4lQuQ4q3RsdgrXbcsbejymW0y7hTRA1ew6+9itllhXrSVtb9s2El8kiDfXvEU+7rKZQL6vW/DC\nb8/rcawajWZgGAiPtEbT5/h8PkpLS9m2bRsVFYUXCmk0mtFLR2MbJVVOe4OoCqvCVl3UV1IVZezM\nGsCZOd1TDJ6oFgPE1lkZ0eku235RUh7K2xc9fIIcS6olTkapRscbWsmm0vK+J+BzvEYV4WB1d9Ri\nWqMZuugy3zBmJHiLdsWeZEmPhnnZE/S8FEbPS8/0x9ycfstiue0tT14/l47GNrmpnuF0LEk6lqTt\n7U0y/9kXDtLR2MaWF9e
      "text/plain": [
       "<matplotlib.figure.Figure at 0x113d958d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "from awips.dataaccess import DataAccessLayer\n",
    "from awips import ThriftClient, RadarCommon\n",
    "from dynamicserialize.dstypes.com.raytheon.uf.common.time import TimeRange\n",
    "from dynamicserialize.dstypes.com.raytheon.uf.common.dataplugin.radar.request import GetRadarDataRecordRequest\n",
    "from datetime import datetime\n",
    "from datetime import timedelta\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from numpy import ma\n",
    "from metpy.plots import ctables\n",
    "import cartopy.crs as ccrs\n",
    "from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER\n",
    "\n",
    "# set EDEX server and radar site definitions\n",
    "site = 'kmux'\n",
    "DataAccessLayer.changeEDEXHost('edex-cloud.unidata.ucar.edu')\n",
    "request = DataAccessLayer.newDataRequest()\n",
    "request.setDatatype('radar')\n",
    "request.setLocationNames(site)\n",
    "\n",
    "# Get latest time for site\n",
    "datatimes = DataAccessLayer.getAvailableTimes(request)\n",
    "#dateTimeStr = str(datatimes[-110])\n",
    "dateTimeStr = \"2017-02-02 03:53:03\"\n",
    "buffer = 60 # seconds\n",
    "dateTime = datetime.strptime(dateTimeStr, '%Y-%m-%d %H:%M:%S')\n",
    "# Build timerange +/- buffer\n",
    "beginRange = dateTime - timedelta(0, buffer)\n",
    "endRange = dateTime + timedelta(0, buffer)\n",
    "timerange = TimeRange(beginRange, endRange)\n",
    "\n",
    "# GetRadarDataRecordRequest to query site with timerange\n",
    "client = ThriftClient.ThriftClient('edex-cloud.unidata.ucar.edu')\n",
    "request = GetRadarDataRecordRequest()\n",
    "request.setTimeRange(timerange)\n",
    "request.setRadarId(site)\n",
    "\n",
    "# Map config\n",
    "def make_map(bbox, projection=ccrs.PlateCarree()):\n",
    "    fig, ax = plt.subplots(figsize=(12, 12),\n",
    "            subplot_kw=dict(projection=projection))\n",
    "    ax.set_extent(bbox)\n",
    "    ax.coastlines(resolution='50m')\n",
    "    gl = ax.gridlines(draw_labels=True)\n",
    "    gl.xlabels_top = gl.ylabels_right = False\n",
    "    gl.xformatter = LONGITUDE_FORMATTER\n",
    "    gl.yformatter = LATITUDE_FORMATTER\n",
    "    return fig, ax\n",
    "\n",
    "# Products\n",
    "nexrad = {}\n",
    "nexrad[\"N0Q\"] = {\n",
    "    'id': 94, \n",
    "    'unit':'dBZ', \n",
    "    'name':'0.5 deg Base Reflectivity', \n",
    "    'ctable': ['NWSStormClearReflectivity',-20., 0.5], \n",
    "    'res': 1000.,\n",
    "    'elev': '0.5'\n",
    "}\n",
    "nexrad[\"N0U\"] = {\n",
    "    'id': 99, \n",
    "    'unit':'kts', \n",
    "    'name':'0.5 deg Base Velocity', \n",
    "    'ctable': ['NWS8bitVel',-100.,1.], \n",
    "    'res': 250.,\n",
    "    'elev': '0.5'\n",
    "}\n",
    "\n",
    "for code in nexrad:\n",
    "    # Colortable filename, beginning value, increment\n",
    "    ctable = nexrad[code]['ctable'][0]\n",
    "    beg = nexrad[code]['ctable'][1]\n",
    "    inc = nexrad[code]['ctable'][2]\n",
    "    # Set product code and elevation angle and send request\n",
    "    request.setProductCode(nexrad[code]['id'])\n",
    "    request.setPrimaryElevationAngle(nexrad[code]['elev'])\n",
    "    response = client.sendRequest(request)\n",
    "    \n",
    "    if response.getData():\n",
    "        for record in response.getData():\n",
    "            \n",
    "            # Get record hdf5 data\n",
    "            idra = record.getHdf5Data()\n",
    "            rdat,azdat,depVals,threshVals = RadarCommon.get_hdf5_data(idra)\n",
    "            dim = rdat.getDimension()\n",
    "            lat,lon = float(record.getLatitude()),float(record.getLongitude())\n",
    "            radials,rangeGates = rdat.getSizes()\n",
    "            \n",
    "            # Convert raw byte to pixel value\n",
    "            rawValue=np.array(rdat.getByteData())\n",
    "            array = []\n",
    "            for rec in rawValue:\n",
    "                if rec<0:\n",
    "                    rec+=256\n",
    "                array.append(rec)        \n",
    "                \n",
    "            if azdat:\n",
    "                azVals = azdat.getFloatData()\n",
    "                az = np.array(RadarCommon.encode_radial(azVals))\n",
    "                dattyp = RadarCommon.get_data_type(azdat)\n",
    "                az = np.append(az,az[-1])\n",
    "\n",
    "            header = RadarCommon.get_header(record, format, rangeGates, radials, azdat, 'description')\n",
    "            rng = np.linspace(0, rangeGates, rangeGates + 1)\n",
    "\n",
    "            # Convert az/range to a lat/lon\n",
    "            from pyproj import Geod\n",
    "            g = Geod(ellps='clrk66')\n",
    "            center_lat = np.ones([len(az),len(rng)])*lat   \n",
    "            center_lon = np.ones([len(az),len(rng)])*lon\n",
    "            az2D = np.ones_like(center_lat)*az[:,None]\n",
    "            rng2D = np.ones_like(center_lat)*np.transpose(rng[:,None])*nexrad[code]['res']\n",
    "            llon,llat,back=g.fwd(center_lon,center_lat,az2D,rng2D)\n",
    "            bbox = [llon.min(), llon.max(), llat.min(), llat.max()]\n",
    "            \n",
    "            # Create 2d array\n",
    "            multiArray = np.reshape(array, (-1, rangeGates))\n",
    "            data = ma.array(multiArray)\n",
    "            \n",
    "            # threshVals[0:2] contains halfwords 31,32,33 (min value, increment, num levels)\n",
    "            data = ma.array(threshVals[0]/10. + (multiArray)*threshVals[1]/10.)\n",
    "            \n",
    "            if nexrad[code]['unit'] == 'kts':\n",
    "                data[data<-63] = ma.masked\n",
    "                data *= 1.94384 # Convert to knots\n",
    "            else:\n",
    "                data[data<=((threshVals[0]/10.)+threshVals[1]/10.)] = ma.masked\n",
    "                \n",
    "            # Create figure\n",
    "            fig, ax = make_map(bbox=bbox)\n",
    "            norm, cmap = ctables.registry.get_with_steps(ctable, beg, inc)\n",
    "            cs = ax.pcolormesh(llon, llat, data, norm=norm, cmap=cmap)\n",
    "            ax.set_aspect('equal', 'datalim')\n",
    "            cbar = plt.colorbar(cs, extend='both', shrink=0.75, orientation='horizontal')\n",
    "            cbar.set_label(site.upper()+\" \"+ str(nexrad[code]['res']/1000.) +\"km \" \\\n",
    "                           +nexrad[code]['name']+\" (\"+code+\") \" \\\n",
    "                           +nexrad[code]['unit']+\" \" \\\n",
    "                           +str(record.getDataTime()))\n",
    "            \n",
    "            # Zoom to within +-2 deg of center\n",
    "            ax.set_xlim(lon-2., lon+2.)\n",
    "            ax.set_ylim(lat-2., lat+2.)\n",
    "            \n",
    "            plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "compare with the same product scan rendered in AWIPS CAVE (slightly different projections and still some color mapping differences, most noticeable in ground clutter).\n",
    "\n",
    "![](http://i.imgur.com/q7zPRod.gif)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}