plot_20210625-MattG.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The following notebook will read and plot XTRRA data from the 25 June 2021 Carroll Co. tornado case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "## You are using the Python ARM Radar Toolkit (Py-ART), an open source\n",
      "## library for working with weather radar data. Py-ART is partly\n",
      "## supported by the U.S. Department of Energy as part of the Atmospheric\n",
      "## Radiation Measurement (ARM) Climate Research Facility, an Office of\n",
      "## Science user facility.\n",
      "##\n",
      "## If you use this software to prepare a publication, please cite:\n",
      "##\n",
      "##     JJ Helmus and SM Collis, JORS 2016, doi: 10.5334/jors.119\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import sys\n",
    "sys.path.append('../scripts')\n",
    "import pytz\n",
    "from pyXTRRA import *\n",
    "from copy import deepcopy\n",
    "\n",
    "#Disable a couple of warnings:\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore', category=UserWarning, append=True)\n",
    "warnings.filterwarnings('ignore', category=RuntimeWarning, append=True)\n",
    "warnings.filterwarnings('ignore', category=FutureWarning, append=True)\n",
    "warnings.filterwarnings('ignore', category=DeprecationWarning)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Read the data files into a list. \n",
    "Each file contains a single elevation scan (aka sweep).\n",
    "\n",
    "### List of tilts: 0.5 0.9 1.3 1.8 2.4 4.0 6.4 10.0 90.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['/depot/rtanama/share/XTRRA/cases/20210625/purdue.in-20210626-001902.netcdf',\n",
       " '/depot/rtanama/share/XTRRA/cases/20210625/purdue.in-20210626-001921.netcdf',\n",
       " '/depot/rtanama/share/XTRRA/cases/20210625/purdue.in-20210626-001940.netcdf']"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Directory where the XTRRA files live:\n",
    "fileDir = '/depot/rtanama/share/XTRRA/cases/20210625/'\n",
    "\n",
    "# Note that I've used a wildcard to limit the plotted files to the time close to the tornado.\n",
    "theFiles = glob.glob(fileDir + \"purdue.in-20210626-0019*.netcdf\")\n",
    "theFiles.sort()\n",
    "theFiles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot reflectivity\n",
    "Note: These plots may not display in the notebook, but you'll find the PNG files in the directory where this notebook lives."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Plotting file  purdue.in-20210626-001902.netcdf\n",
      "Saving file  purdue.in-20210626-001902.Z.png\n",
      "Plotting file  purdue.in-20210626-001921.netcdf\n",
      "Saving file  purdue.in-20210626-001921.Z.png\n",
      "Plotting file  purdue.in-20210626-001940.netcdf\n",
      "Saving file  purdue.in-20210626-001940.Z.png\n"
     ]
    }
   ],
   "source": [
    "for f in theFiles:\n",
    "    print(\"Plotting file \",f[-32:])\n",
    "    fig = plt.figure(figsize=(8,6))\n",
    "    plt.clf()\n",
    "    radar = CreateRadarFromRXM25netCDF(f,'XTRRA')\n",
    "    display = pyart.graph.RadarMapDisplay(radar)\n",
    "    display.plot_ppi_map('reflectivity', 0, vmin=-10, vmax=60, \n",
    "                         min_lon = -86.95, max_lon = -86.5, \n",
    "                         min_lat = 40.3, max_lat = 40.6,\n",
    "                        cmap = pyart.graph.cm.LangRainbow12)\n",
    "    # If you have a county shapefile, read it into a shapely feature \n",
    "    # called \"counties\" (defined in pyXTRRA.py).\n",
    "    #cy = display.ax.add_feature(counties, edgecolor='brown', alpha=0.7, lw=0.75)\n",
    "    for r in np.arange(10,60,10):\n",
    "        display.plot_range_ring(r, line_style='k--')\n",
    "    display.plot_point(radar.longitude['data'][0],radar.latitude['data'][0],'k*',label_text=\"XTRRA\", label_offset =[0.01,-0.01])\n",
    "    plt.show()\n",
    "    figString = f[-32:-7]+\".Z.png\"\n",
    "    print(\"Saving file \",figString)\n",
    "    plt.savefig(figString,dpi=600)\n",
    "    # You can uncomment the following line if you have ImageMagick installed in your path.\n",
    "    # It removes white space from around the edges of your figure.\n",
    "    # os.system('convert -trim %s %s'%(figString,figString))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot Doppler (radial) velocity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Plotting file  purdue.in-20210626-001902.netcdf\n",
      "Saving file  purdue.in-20210626-001902.Vr.png\n",
      "Plotting file  purdue.in-20210626-001921.netcdf\n",
      "Saving file  purdue.in-20210626-001921.Vr.png\n",
      "Plotting file  purdue.in-20210626-001940.netcdf\n",
      "Saving file  purdue.in-20210626-001940.Vr.png\n"
     ]
    }
   ],
   "source": [
    "for f in theFiles:\n",
    "    print(\"Plotting file \",f[-32:])\n",
    "    fig = plt.figure(figsize=(8,6))\n",
    "    radar = CreateRadarFromRXM25netCDF(f,'XTRRA')\n",
    "    display = pyart.graph.RadarMapDisplay(radar)\n",
    "    display.plot_ppi_map('velocity', 0, vmin=-30, vmax=30, \n",
    "                         min_lon = -87.0, max_lon = -86.55, \n",
    "                         min_lat = 40.3, max_lat = 40.6)\n",
    "    #cy = display.ax.add_feature(counties, edgecolor='brown', alpha=0.7, lw=0.75)\n",
    "    for r in np.arange(10,60,10):\n",
    "        display.plot_range_ring(r, line_style='k--')\n",
    "    display.plot_point(radar.longitude['data'][0],radar.latitude['data'][0],'k*',label_text=\"XTRRA\", label_offset =[0.01,-0.01])\n",
    "    figString = f[-32:-7]+\".Vr.png\"\n",
    "    print(\"Saving file \",figString)\n",
    "    plt.savefig(figString,dpi=600)\n",
    "    plt.close(fig)\n",
    "    #os.system('convert -trim %s %s'%(figString,figString))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot ZDR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Plotting file  purdue.in-20210626-001902.netcdf\n",
      "Saving file  purdue.in-20210626-001902.Zdr.png\n",
      "Plotting file  purdue.in-20210626-001921.netcdf\n",
      "Saving file  purdue.in-20210626-001921.Zdr.png\n",
      "Plotting file  purdue.in-20210626-001940.netcdf\n",
      "Saving file  purdue.in-20210626-001940.Zdr.png\n"
     ]
    }
   ],
   "source": [
    "for f in theFiles:\n",
    "    print(\"Plotting file \",f[-32:])\n",
    "    fig = plt.figure(figsize=(8,6))\n",
    "    radar = CreateRadarFromRXM25netCDF(f,'XTRRA')\n",
    "    display = pyart.graph.RadarMapDisplay(radar)\n",
    "    display.plot_ppi_map('differential_reflectivity', 0, vmin=-2, vmax=6, \n",
    "                         min_lon = -87.0, max_lon = -86.55, \n",
    "                         min_lat = 40.3, max_lat = 40.6,\n",
    "                        cmap = pyart.graph.cm.LangRainbow12)\n",
    "    #cy = display.ax.add_feature(counties, edgecolor='brown', alpha=0.7, lw=0.75)\n",
    "    for r in np.arange(10,60,10):\n",
    "        display.plot_range_ring(r, line_style='k--')\n",
    "    display.plot_point(radar.longitude['data'][0],radar.latitude['data'][0],'k*',label_text=\"XTRRA\", label_offset =[0.01,-0.01])\n",
    "    figString = f[-32:-7]+\".Zdr.png\"\n",
    "    print(\"Saving file \",figString)\n",
    "    plt.savefig(figString,dpi=600)\n",
    "    plt.close(fig)\n",
    "    #os.system('convert -trim %s %s'%(figString,figString))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot KDP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Plotting file  purdue.in-20210626-001902.netcdf\n",
      "Saving file  purdue.in-20210626-001902.Kdp.png\n",
      "Plotting file  purdue.in-20210626-001921.netcdf\n",
      "Saving file  purdue.in-20210626-001921.Kdp.png\n",
      "Plotting file  purdue.in-20210626-001940.netcdf\n",
      "Saving file  purdue.in-20210626-001940.Kdp.png\n"
     ]
    }
   ],
   "source": [
    "for f in theFiles:\n",
    "    print(\"Plotting file \",f[-32:])\n",
    "    fig = plt.figure(figsize=(8,6))\n",
    "    radar = CreateRadarFromRXM25netCDF(f,'XTRRA')\n",
    "    display = pyart.graph.RadarMapDisplay(radar)\n",
    "    display.plot_ppi_map('specific_differential_phase', 0, vmin=0, vmax=6, \n",
    "                         min_lon = -87.0, max_lon = -86.55, \n",
    "                         min_lat = 40.3, max_lat = 40.6,\n",
    "                        cmap = pyart.graph.cm.LangRainbow12)\n",
    "    #cy = display.ax.add_feature(counties, edgecolor='brown', alpha=0.7, lw=0.75)\n",
    "    for r in np.arange(10,60,10):\n",
    "        display.plot_range_ring(r, line_style='k--')\n",
    "    display.plot_point(radar.longitude['data'][0],radar.latitude['data'][0],'k*',label_text=\"XTRRA\", label_offset =[0.01,-0.01])\n",
    "    #display.ax.add_feature(states, edgecolor='black', alpha=0.7, lw=1.5)\n",
    "    figString = f[-32:-7]+\".Kdp.png\"\n",
    "    print(\"Saving file \",figString)\n",
    "    plt.savefig(figString,dpi=600)\n",
    "    plt.close(fig)\n",
    "    #os.system('convert -trim %s %s'%(figString,figString))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot CC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Plotting file  purdue.in-20210626-001902.netcdf\n",
      "Saving file  purdue.in-20210626-001902.rhv.png\n",
      "Plotting file  purdue.in-20210626-001921.netcdf\n",
      "Saving file  purdue.in-20210626-001921.rhv.png\n",
      "Plotting file  purdue.in-20210626-001940.netcdf\n",
      "Saving file  purdue.in-20210626-001940.rhv.png\n"
     ]
    }
   ],
   "source": [
    "for f in theFiles:\n",
    "    print(\"Plotting file \",f[-32:])\n",
    "    fig = plt.figure(figsize=(8,6))\n",
    "    radar = CreateRadarFromRXM25netCDF(f,'XTRRA')\n",
    "    display = pyart.graph.RadarMapDisplay(radar)\n",
    "    display.plot_ppi_map('cross_correlation_ratio', 0, vmin=0, vmax=1, \n",
    "                         min_lon = -87.0, max_lon = -86.55, \n",
    "                         min_lat = 40.3, max_lat = 40.6,\n",
    "                        cmap = 'cividis')\n",
    "    #cy = display.ax.add_feature(counties, edgecolor='brown', alpha=0.7, lw=0.75)\n",
    "    for r in np.arange(10,60,10):\n",
    "        display.plot_range_ring(r, line_style='k--')\n",
    "    display.plot_point(radar.longitude['data'][0],radar.latitude['data'][0],'k*',label_text=\"XTRRA\", label_offset =[0.01,-0.01])\n",
    "    #display.ax.add_feature(states, edgecolor='black', alpha=0.7, lw=1.5)\n",
    "    figString = f[-32:-7]+\".rhv.png\"\n",
    "    print(\"Saving file \",figString)\n",
    "    plt.savefig(figString,dpi=600)\n",
    "    plt.close(fig)\n",
    "    #os.system('convert -trim %s %s'%(figString,figString))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Calculate & plot azimuthal shear\n",
    "The azimuthal shear routine is adopted from PyMeso by Jordan Brook: https://github.com/jordanbrook/PyMeso\n",
    "It uses a parallelization module called \"numba\" to speed up the calculation of shear.\n",
    "\n",
    "## Looking at the Doppler velocity fields (above), we find that they are really noisy!\n",
    "I've already done a cursory filtering of the data where SNR < 3 dB, but more filtering is required to get rid of the clutter near the radar. We can accomplish that using the \"gatefilter\" class in Py-ART."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Plotting file  purdue.in-20210626-001902.netcdf\n",
      "Saving file  purdue.in-20210626-001902.text_Vr.png\n"
     ]
    }
   ],
   "source": [
    "f = theFiles[0] # Just do the first file for now.\n",
    "print(\"Plotting file \",f[-32:])\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "radar = CreateRadarFromRXM25netCDF(f,'XTRRA')\n",
    "# Restore the entire velocity field (i.e., remove the mask)\n",
    "radar.fields['velocity']['data'] = radar.fields['velocity']['data'].data\n",
    "# Calculate texture of velocity\n",
    "texture = pyart.retrieve.calculate_velocity_texture(radar, vel_field='velocity',\n",
    "                                                    wind_size=3, nyq = 33.44)\n",
    "radar.add_field('velocity_texture',texture,replace_existing=True)\n",
    "display = pyart.graph.RadarMapDisplay(radar)\n",
    "display.plot_ppi_map('velocity_texture', 0, vmin=0, vmax=15, \n",
    "                     min_lon = -87.0, max_lon = -86.55, \n",
    "                     min_lat = 40.3, max_lat = 40.6)\n",
    "#cy = display.ax.add_feature(counties, edgecolor='brown', alpha=0.7, lw=0.75)\n",
    "for r in np.arange(10,60,10):\n",
    "    display.plot_range_ring(r, line_style='k--')\n",
    "display.plot_point(radar.longitude['data'][0],radar.latitude['data'][0],'k*',label_text=\"XTRRA\", label_offset =[0.01,-0.01])\n",
    "figString = f[-32:-7]+\".text_Vr.png\"\n",
    "print(\"Saving file \",figString)\n",
    "plt.savefig(figString,dpi=600)\n",
    "plt.close(fig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Plotting file  purdue.in-20210626-001902.netcdf\n",
      "Saving file  purdue.in-20210626-001902.corr_Vr.png\n"
     ]
    }
   ],
   "source": [
    "# Create an instance of class GateFilter, tied to the current radar object.\n",
    "gf = pyart.filters.GateFilter(radar)\n",
    "# Add exclusion conditions to the GateFilter.\n",
    "gf.exclude_below('cross_correlation_ratio',0.7)\n",
    "gf.exclude_above('velocity_texture',8)  #Value found by trial and error.\n",
    "# Add despeckling (removal of isolated gates) to the GateFilter.\n",
    "gf_despeckle = pyart.correct.despeckle_field(radar, 'velocity', gatefilter=gf)\n",
    "# Create a true (deep) copy of the velocity field (a dictionary), \n",
    "# to hold the filtered, despeckled data, and apply the GateFilter to it.\n",
    "# A 'shallow' copy just creates a pointer back to the original data.\n",
    "corr_velocity = deepcopy(radar.fields['velocity'])\n",
    "corr_velocity_data = np.ma.masked_where(gf_despeckle.gate_included == False, \n",
    "                                        radar.fields['velocity']['data'])\n",
    "corr_velocity['data'] = corr_velocity_data\n",
    "corr_velocity['long_name'] = 'Mean dopper velocity (corrected)'\n",
    "radar.add_field('corr_velocity',corr_velocity,replace_existing=True)\n",
    "\n",
    "# Plot the resulting azimuthal shear field\n",
    "print(\"Plotting file \",f[-32:])\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "display = pyart.graph.RadarMapDisplay(radar)\n",
    "display.plot_ppi_map('corr_velocity', 0, vmin=-30, vmax=30, \n",
    "                     min_lon = -87.0, max_lon = -86.55, \n",
    "                     min_lat = 40.3, max_lat = 40.6, cmap = 'seismic')\n",
    "#cy = display.ax.add_feature(counties, edgecolor='brown', alpha=0.7, lw=0.75)\n",
    "for r in np.arange(10,60,10):\n",
    "    display.plot_range_ring(r, line_style='k--')\n",
    "display.plot_point(radar.longitude['data'][0],radar.latitude['data'][0],'k*',label_text=\"XTRRA\", label_offset =[0.01,-0.01])\n",
    "figString = f[-32:-7]+\".corr_Vr.png\"\n",
    "print(\"Saving file \",figString)\n",
    "plt.savefig(figString,dpi=600)\n",
    "plt.close(fig)\n",
    "#os.system('convert -trim %s %s'%(figString,figString))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Now, finally, we can apply the azimuthal shear routine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'data': masked_array(\n",
       "   data=[[0.0, nan, 0.0028114644810557365, ..., --, --, --],\n",
       "         [0.0, nan, 0.002844721544533968, ..., --, --, --],\n",
       "         [--, nan, 0.002992850262671709, ..., --, --, --],\n",
       "         ...,\n",
       "         [--, nan, --, ..., --, --, --],\n",
       "         [0.0, nan, --, ..., --, --, --],\n",
       "         [0.0, nan, --, ..., --, --, --]],\n",
       "   mask=[[False, False, False, ...,  True,  True,  True],\n",
       "         [False, False, False, ...,  True,  True,  True],\n",
       "         [ True, False, False, ...,  True,  True,  True],\n",
       "         ...,\n",
       "         [ True, False,  True, ...,  True,  True,  True],\n",
       "         [False, False,  True, ...,  True,  True,  True],\n",
       "         [False, False,  True, ...,  True,  True,  True]],\n",
       "   fill_value=1e+20,\n",
       "   dtype=float32),\n",
       " 'long_name': 'LLSD Azimuthal Shear',\n",
       " 'standard_name': 'azimuthal_shear',\n",
       " '_FillValue': nan,\n",
       " '_Least_significant_digit': 2,\n",
       " 'comment': 'LLSD azimuthal shear calculation from Miller, M. L., Lakshmanan, V., and Smith, T. M. (2013). An Automated Method for Depicting Mesocyclone Paths and Intensities. Weather and Forecasting, 28(3): 570-585.',\n",
       " 'units': 'Hz'}"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "azi_shear_meta = llsd_main(radar, 'reflectivity', 'corr_velocity')\n",
    "azi_shear_meta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This just returns a dictionary; we have to add the field back to the radar object and then plot it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Plotting file  purdue.in-20210626-001902.netcdf\n",
      "Saving file  purdue.in-20210626-001902.az_shear.png\n"
     ]
    }
   ],
   "source": [
    "radar.add_field('az_shear',azi_shear_meta,replace_existing=True)\n",
    "print(\"Plotting file \",f[-32:])\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "display = pyart.graph.RadarMapDisplay(radar)\n",
    "display.plot_ppi_map('az_shear', 0, vmin=-0.003, vmax=0.003, \n",
    "                     min_lon = -87.0, max_lon = -86.55, \n",
    "                     min_lat = 40.3, max_lat = 40.6, cmap = 'seismic')\n",
    "#cy = display.ax.add_feature(counties, edgecolor='brown', alpha=0.7, lw=0.75)\n",
    "for r in np.arange(10,60,10):\n",
    "    display.plot_range_ring(r, line_style='k--')\n",
    "display.plot_point(radar.longitude['data'][0],radar.latitude['data'][0],'k*',label_text=\"XTRRA\", label_offset =[0.01,-0.01])\n",
    "figString = f[-32:-7]+\".az_shear.png\"\n",
    "print(\"Saving file \",figString)\n",
    "plt.savefig(figString,dpi=600)\n",
    "plt.close(fig)\n",
    "#os.system('convert -trim %s %s'%(figString,figString))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (XTRRA)",
   "language": "python",
   "name": "xtrra"
  },
  "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}