diff --git a/pages/gallery/plotting_nexrad_radar.ipynb b/pages/gallery/plotting_nexrad_radar.ipynb new file mode 100644 index 00000000..061e6e1e --- /dev/null +++ b/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 +}