example_notebooks/pyart_nexrad_maps_reports.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get NEXRAD data from AWS, get severe reports from SPC archive, and plot a loop of them! Here we use the August 2020 derecho as an example\n",
    "------------\n",
    "\n",
    "#### Pulling radar data follows from: https://nexradaws.readthedocs.io/en/latest/Tutorial.html\n",
    "#### Radar plotting follows from:  https://github.com/ARM-DOE/pyart/blob/master/examples/plotting/plot_ppi_with_rings.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors as mcolors\n",
    "import matplotlib.dates as mdates\n",
    "\n",
    "import cartopy.crs as ccrs\n",
    "import pyart\n",
    "import pandas as pd\n",
    "\n",
    "import nexradaws\n",
    "import tempfile\n",
    "import pytz\n",
    "\n",
    "templocation = tempfile.mkdtemp()\n",
    "\n",
    "import cartopy.feature as cfeature\n",
    "from metpy.plots import USCOUNTIES\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Get NEXRAD data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Define the radar, start time and end time\n",
    "radar_id = 'KDVN'\n",
    "start = pd.Timestamp(2020,8,10,16,30).tz_localize(tz='UTC')\n",
    "end = pd.Timestamp(2020,8,10,21,0).tz_localize(tz='UTC')\n",
    "\n",
    "### Bounds of map we want to plot\n",
    "min_lon = -93.25\n",
    "max_lon = -88.\n",
    "min_lat = 40.35\n",
    "max_lat = 43.35"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#### and get the data\n",
    "conn = nexradaws.NexradAwsInterface()\n",
    "scans = conn.get_avail_scans_in_range(start, end, radar_id)\n",
    "print(\"There are {} scans available between {} and {}\\n\".format(len(scans), start, end))\n",
    "print(scans[0:4])\n",
    "\n",
    "## download these files\n",
    "#results = conn.download(scans[0:2], templocation)\n",
    "results = conn.download(scans, templocation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now get the severe reports from the SPC site. This assumes you're plotting a year far enough in the past that SPC has official records available. If plotting a more recent time period, then the local storm reports archive at IEM is a good source"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### wind reports\n",
    "wind_rpts = pd.read_csv(\"https://www.spc.noaa.gov/wcm/data/\"+str(start.year)+\"_wind.csv\")\n",
    "\n",
    "wind_rpts['datetime'] = pd.to_datetime(wind_rpts.date + ' ' + wind_rpts.time) ## convert to datetime\n",
    "wind_rpts.set_index(\"datetime\",inplace=True)\n",
    "\n",
    "### times in the file are given in central standard time (UTC+6). Localize, and convert to UTC\n",
    "wind_rpts.index = wind_rpts.index.tz_localize(\"Etc/GMT+6\",ambiguous='NaT',nonexistent='shift_forward').tz_convert(\"UTC\")\n",
    "\n",
    "## subset down to 30 minutes before/after the radar times we're plotting\n",
    "wind_rpts = wind_rpts[((start-pd.Timedelta(minutes=30)).strftime(\"%Y-%m-%d %H:%M\")):((end+pd.Timedelta(minutes=30)).strftime(\"%Y-%m-%d %H:%M\"))]\n",
    "wind_rpts\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### repeat for tornado reports\n",
    "tor_rpts = pd.read_csv(\"https://www.spc.noaa.gov/wcm/data/\"+str(start.year)+\"_torn.csv\")\n",
    "\n",
    "tor_rpts['datetime'] = pd.to_datetime(tor_rpts.date + ' ' + tor_rpts.time) ## convert to datetime\n",
    "tor_rpts.set_index(\"datetime\",inplace=True)\n",
    "\n",
    "### times in the file are given in central standard time (UTC+6). Localize, and convert to UTC\n",
    "tor_rpts.index = tor_rpts.index.tz_localize(\"Etc/GMT+6\",ambiguous='NaT',nonexistent='shift_forward').tz_convert(\"UTC\")\n",
    "\n",
    "## subset down to 30 minutes before/after the radar times we're plotting\n",
    "tor_rpts = tor_rpts[((start-pd.Timedelta(minutes=30)).strftime(\"%Y-%m-%d %H:%M\")):((end+pd.Timedelta(minutes=30)).strftime(\"%Y-%m-%d %H:%M\"))]\n",
    "tor_rpts\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### repeat for hail\n",
    "hail_rpts = pd.read_csv(\"https://www.spc.noaa.gov/wcm/data/\"+str(start.year)+\"_hail.csv\")\n",
    "\n",
    "hail_rpts['datetime'] = pd.to_datetime(hail_rpts.date + ' ' + hail_rpts.time) ## convert to datetime\n",
    "hail_rpts.set_index(\"datetime\",inplace=True)\n",
    "\n",
    "### times in the file are given in central standard time (UTC+6). Localize, and convert to UTC\n",
    "hail_rpts.index = hail_rpts.index.tz_localize(\"Etc/GMT+6\",ambiguous='NaT',nonexistent='shift_forward').tz_convert(\"UTC\")\n",
    "\n",
    "## subset down to 30 minutes before/after the radar times we're plotting\n",
    "hail_rpts = hail_rpts[((start-pd.Timedelta(minutes=30)).strftime(\"%Y-%m-%d %H:%M\")):((end+pd.Timedelta(minutes=30)).strftime(\"%Y-%m-%d %H:%M\"))]\n",
    "hail_rpts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now, plot the maps (reflectivity with severe reports)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "### loop over the radar images that have been downloaded\n",
    "\n",
    "for i,scan in enumerate(results.iter_success(),start=1):\n",
    "#for i in range(0,1):\n",
    "    ## skip the files ending in \"MDM\"\n",
    "    if scan.filename[-3:] != \"MDM\":\n",
    "        \n",
    "        print(str(i))\n",
    "        print(\"working on \"+scan.filename)\n",
    "        \n",
    "        this_time = pd.to_datetime(scan.filename[4:17], format=\"%Y%m%d_%H%M\").tz_localize(\"UTC\")\n",
    "    \n",
    "        radar = scan.open_pyart()\n",
    "        #display = pyart.graph.RadarDisplay(radar)\n",
    "\n",
    "        fig = plt.figure(figsize=[15, 7])\n",
    "\n",
    "        map_panel_axes = [0.05, 0.05, .4, .80]\n",
    "        x_cut_panel_axes = [0.55, 0.10, .4, .25]\n",
    "        y_cut_panel_axes = [0.55, 0.50, .4, .25]\n",
    "\n",
    "        projection = ccrs.PlateCarree()\n",
    "\n",
    "        ## apply gatefilter (see here: https://arm-doe.github.io/pyart/notebooks/masking_data_with_gatefilters.html)\n",
    "        #gatefilter = pyart.correct.moment_based_gate_filter(radar)\n",
    "        gatefilter = pyart.filters.GateFilter(radar)\n",
    "\n",
    "        # Lets remove reflectivity values below a threshold.\n",
    "        gatefilter.exclude_below('reflectivity', -2.5)\n",
    "\n",
    "        display = pyart.graph.RadarMapDisplay(radar)\n",
    "       \n",
    "        ### set up plot\n",
    "        ax1 = fig.add_axes(map_panel_axes, projection=projection)\n",
    "        \n",
    "        # Add some various map elements to the plot to make it recognizable.\n",
    "        ax1.add_feature(USCOUNTIES.with_scale('500k'), edgecolor=\"gray\", linewidth=0.4)\n",
    "        #ax1.coastlines('50m', edgecolor='black', linewidth=0.75)\n",
    "        ax1.add_feature(cfeature.STATES.with_scale('10m'), linewidth=1.0)\n",
    "    \n",
    "        cf = display.plot_ppi_map('reflectivity', 0, vmin=-7.5, vmax=65,\n",
    "                          min_lon=min_lon, max_lon=max_lon, min_lat=min_lat, max_lat=max_lat,\n",
    "                          title=radar_id+\" reflectivity and severe weather reports, \"+this_time.strftime(\"%H%M UTC %d %b %Y\"),\n",
    "                          projection=projection, resolution='10m',\n",
    "                             gatefilter=gatefilter,\n",
    "                          cmap='pyart_HomeyerRainbow', \n",
    "                          colorbar_flag=False,\n",
    "                          lat_lines=[0,0], lon_lines=[0,0]) ## turns off lat/lon grid lines\n",
    "        #display.plot_crosshairs(lon=lon, lat=lat)\n",
    "        \n",
    "        ## plot horizontal colorbar \n",
    "        display.plot_colorbar(cf,orient='horizontal', pad=0.07)\n",
    "        \n",
    "        # Plot range rings if desired\n",
    "        #display.plot_range_ring(25., color='gray', linestyle='dashed')\n",
    "        #display.plot_range_ring(50., color='gray', linestyle='dashed')\n",
    "        #display.plot_range_ring(100., color='gray', linestyle='dashed')\n",
    "\n",
    "        ax1.set_xticks(np.arange(min_lon, max_lon, .5), crs=ccrs.PlateCarree())\n",
    "        ax1.set_yticks(np.arange(min_lat, max_lat, .5), crs=ccrs.PlateCarree())     \n",
    "       \n",
    "        ## add marker points for severe reports\n",
    "        wind_rpts_now = wind_rpts[((start-pd.Timedelta(minutes=30)).strftime(\"%Y-%m-%d %H:%M\")):this_time.strftime(\"%Y-%m-%d %H:%M\")]\n",
    "        ax1.scatter(wind_rpts_now.slon.values.tolist(), wind_rpts_now.slat.values.tolist(), s=20, facecolors='none', edgecolors='mediumblue', linewidths=1.8)\n",
    "        tor_rpts_now = tor_rpts[((start-pd.Timedelta(minutes=30)).strftime(\"%Y-%m-%d %H:%M\")):this_time.strftime(\"%Y-%m-%d %H:%M\")]\n",
    "        ax1.scatter(tor_rpts_now.slon.values.tolist(), tor_rpts_now.slat.values.tolist(), s=20, facecolors='red', edgecolors='black', marker=\"v\",linewidths=1.5)\n",
    "        hail_rpts_now = hail_rpts[((start-pd.Timedelta(minutes=30)).strftime(\"%Y-%m-%d %H:%M\")):this_time.strftime(\"%Y-%m-%d %H:%M\")]\n",
    "        ax1.scatter(hail_rpts_now.slon.values.tolist(), hail_rpts_now.slat.values.tolist(), s=20, facecolors='none', edgecolors='green', linewidths=1.8)\n",
    "\n",
    "        \n",
    "        plt.savefig(scan.radar_id+\"_\"+scan.filename[4:17]+\"_dz_rpts.png\",bbox_inches='tight',dpi=300,\n",
    "                   facecolor='white', transparent=False)\n",
    "        #plt.show()\n",
    "        plt.close('all')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}