add nexrad plotting example

Author: kgoebber

pages/gallery/plotting_nexrad_radar.ipynb

@@ -0,0 +1,161 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from datetime import datetime, timedelta\n",
+    "\n",
+    "import cartopy.crs as ccrs\n",
+    "import cartopy.feature as cfeature\n",
+    "import matplotlib.pyplot as plt\n",
+    "from metpy.plots import colortables, USSTATES, USCOUNTIES\n",
+    "import numpy as np\n",
+    "from pyproj import Proj\n",
+    "from siphon.catalog import TDSCatalog\n",
+    "import xarray as xr\n",
+    "\n",
+    "import warnings\n",
+    "warnings.filterwarnings(\"ignore\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def get_radar_file_url(datasets, date):\n",
+    "    '''A function to help find the desired satellite data based on the time given.\n",
+    "    \n",
+    "    Input:\n",
+    "     - List of datasets from a THREDDS Catalog\n",
+    "     - Date of desired dataset (datetime object)\n",
+    "     \n",
+    "    Output:\n",
+    "     - Index value of dataset closest to desired time\n",
+    "    '''\n",
+    "    rad_date_hour = date.strftime('%Y%m%d_%H')\n",
+    "    files = []\n",
+    "    times = []\n",
+    "    for file in cat.datasets:\n",
+    "        if rad_date_hour in file:\n",
+    "            times.append(datetime.strptime(file[-18:-5], '%Y%m%d_%H%M'))\n",
+    "            files.append(file)\n",
+    "    if not files:\n",
+    "        date = date - timedelta(hours=1)\n",
+    "        rad_date_hour = date.strftime('%Y%m%d_%H')\n",
+    "        for file in cat.datasets:\n",
+    "            if rad_date_hour in file:\n",
+    "                times.append(datetime.strptime(file[-18:-5], '%Y%m%d_%H%M'))\n",
+    "                files.append(file)\n",
+    "    find_file = np.abs(np.array(times) - date)\n",
+    "    return list(cat.datasets).index(files[np.argmin(find_file)])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "date = datetime.utcnow()\n",
+    "\n",
+    "# Create variables for URL generation\n",
+    "station = 'KLOT'\n",
+    "\n",
+    "# Construct the data_url string\n",
+    "data_url = (f'https://thredds.ucar.edu/thredds/catalog/nexrad/level2/'\n",
+    "            f'{station}/{date:%Y%m%d}/catalog.html')\n",
+    "\n",
+    "# Get list of files available for particular day\n",
+    "cat = TDSCatalog(data_url)\n",
+    "\n",
+    "# Use homemade function to get dataset for desired time\n",
+    "dataset = cat.datasets[get_radar_file_url(cat.datasets, date)]\n",
+    "\n",
+    "ds = xr.open_dataset(dataset.access_urls['OPENDAP'], decode_times=False,\n",
+    "                     decode_coords=False, mask_and_scale=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "station = ds.Station\n",
+    "slat = ds.StationLatitude\n",
+    "slon = ds.StationLongitude\n",
+    "elevation = ds.StationElevationInMeters\n",
+    "vtime = datetime.strptime(ds.time_coverage_start, '%Y-%m-%dT%H:%M:%SZ')\n",
+    "\n",
+    "dataproj = Proj(f\"+proj=stere +lat_0={slat} +lat_ts={slat} +lon_0={slon} +ellps=WGS84 +units=m\")\n",
+    "\n",
+    "sweep = 0\n",
+    "rng = ds.distanceR_HI.values\n",
+    "az = ds.azimuthR_HI.values[sweep]\n",
+    "ref = ds.Reflectivity_HI.values[sweep]\n",
+    "\n",
+    "x = rng * np.sin(np.deg2rad(az))[:, None]\n",
+    "y = rng * np.cos(np.deg2rad(az))[:, None]\n",
+    "\n",
+    "lon, lat = dataproj(x, y, inverse=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "cmap = colortables.get_colortable('NWSStormClearReflectivity')\n",
+    "\n",
+    "fig, ax = plt.subplots(1, 1, figsize=(12, 10), subplot_kw=dict(projection=ccrs.PlateCarree()))\n",
+    "\n",
+    "img = ax.pcolormesh(lon, lat, ref, vmin=-30, vmax=79, cmap=cmap)\n",
+    "plt.colorbar(img, aspect=50, pad=0)\n",
+    "\n",
+    "ax.set_extent([slon-2.5, slon+2.5, slat-2.5, slat+2.5], ccrs.PlateCarree())\n",
+    "ax.set_aspect('equal', 'datalim')\n",
+    "\n",
+    "ax.add_feature(USCOUNTIES.with_scale('5m'), edgecolor='darkgrey')\n",
+    "ax.add_feature(USSTATES.with_scale('5m'))\n",
+    "\n",
+    "plt.title(f'{station}: {ds.Reflectivity_HI.name}', loc='left')\n",
+    "plt.title(f'Valid Time: {vtime}', loc='right')\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "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.8.2"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}