{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Dependencies | Version\n", "--- | ---\n", "satpy | 0.10.0\n", "pyresample | 1.10.1\n", "metpy | 0.9.1\n", "siphon | 0.8.0\n", "xarray | 0.10\n", "fiona | 1.7.13\n", "\n", "# Hurricane Florence Plot with SatPy, MetPy, and Siphon with THREDDS Server\n", "\n", "This notebook provides a complex example of how SatPy, MetPy, and Siphon can be used to load data from different sources and combine them in to one image using cartopy and matplotlib.\n", "\n", "This example is based on an [example notebook](https://gist.github.com/dopplershift/48001c3102b1583b78c2c6542618beac) originally created by Ryan May (@dopplershift)." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import metpy.calc as mpcalc\n", "from metpy.plots import StationPlot, add_metpy_logo, add_unidata_logo, add_timestamp\n", "from metpy.units import units\n", "from siphon.catalog import TDSCatalog\n", "from siphon.simplewebservice.ndbc import NDBC\n", "\n", "import cartopy.crs as ccrs\n", "import cartopy.feature as cfeature\n", "import matplotlib.patheffects as mpatheffects\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "from satpy import Scene\n", "from satpy.writers import get_enhanced_image\n", "\n", "from urllib.request import urlopen\n", "import fiona\n", "\n", "def get_zip(url):\n", " data = urlopen(url).read()\n", " return fiona.BytesCollection(data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 1: Load GOES-16 ABI from a remote THREDDS server using SatPy and Siphon" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "base_cat_url = 'http://thredds-test.unidata.ucar.edu/thredds/catalog/satellite/goes16/{platform}/{sector}/{channel}/current/catalog.xml'\n", "urls = []\n", "# Access ABI channels 1-3 to use in making a true color image (ABI has 16 total channels)\n", "for channel in ['Channel{:02d}'.format(x) for x in range(1, 4)]:\n", " cat_url = base_cat_url.format(platform='GOES16', sector='CONUS', channel=channel)\n", " cat = TDSCatalog(cat_url)\n", " # get the latest dataset from each\n", " ds = cat.datasets[-1]\n", " # get the opendap url for this dataset\n", " urls.append(ds.access_urls['OPENDAP']) \n", "\n", "# Access the files to figure out what is available and load a true_color RGB\n", "scn = Scene(reader='abi_scmi', filenames=urls)\n", "scn.load(['true_color']) \n", "new_scn = scn.resample(resampler='native')\n", "# Scale the reflectance data to look better as an RGB on a 0 to 1 scale\n", "var = get_enhanced_image(new_scn['true_color']).data\n", "# Get true color data to use later and reorder the dimensions so matplotlib can use the image\n", "# Sadly, this operation is not lazy (bad performance) in xarray at the time of writing\n", "var = var.transpose('y', 'x', 'bands')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 2: Load other data using MetPy and Siphon" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "best_track = get_zip('https://www.nhc.noaa.gov/gis/best_track/al062018_best_track.zip')\n", "forecast = get_zip('https://www.nhc.noaa.gov/gis/forecast/archive/al062018_5day_latest.zip')\n", "track_x, track_y = np.array(list(zip(*([item['geometry']['coordinates'] for _, item in best_track.items()]\n", " + list(forecast.items())[0][1]['geometry']['coordinates']))))" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "buoy = NDBC.latest_observations()\n", "buoy.dropna(subset=['pressure', 'wind_speed', 'wind_direction'], inplace=True)\n", "buoy_u, buoy_v = mpcalc.wind_components(buoy['wind_speed'].values, buoy['wind_direction'].values * units.degree)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "sst_cat = TDSCatalog('https://www.ncei.noaa.gov/thredds/catalog/OisstBase/NetCDF/AVHRR/201809/catalog.xml')\n", "sst = sst_cat.datasets[-1].remote_access(use_xarray=True)\n", "sst_data = sst.metpy.parse_cf('sst')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 3: Plot all of the data using Cartopy" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "%matplotlib inline\n", "fig = plt.figure(figsize=(20, 10), dpi=200)\n", "abi_crs = var.attrs['area'].to_cartopy_crs()\n", "ax = fig.add_subplot(1, 1, 1, projection=abi_crs)\n", "\n", "xy = abi_crs.transform_points(ccrs.PlateCarree(), buoy['longitude'].values, buoy['latitude'].values)\n", "mask = mpcalc.reduce_point_density(xy[..., 0:2], 30000, priority=buoy['pressure'].values)\n", "\n", "kwargs = dict(path_effects=[mpatheffects.withStroke(foreground='white', linewidth=2)],\n", " clip_on=True)\n", "sp = StationPlot(ax, buoy['longitude'].values[mask], buoy['latitude'].values[mask], transform=ccrs.PlateCarree())\n", "sp.plot_parameter('NW', buoy['air_temperature'].values[mask], color='red', **kwargs)\n", "sp.plot_parameter('SW', buoy['dewpoint'].values[mask], color='darkgreen', **kwargs)\n", "sp.plot_parameter('NE', buoy['pressure'].values[mask], color='black', **kwargs)\n", "sp.plot_barb(buoy_u[mask], buoy_v[mask], edgecolor='white')\n", "\n", "ax.imshow(var.data, extent=(var.x[0], var.x[-1], var.y[-1], var.y[0]), origin='upper')\n", "ax.contour(sst_data.lon, sst_data.lat, sst_data.squeeze(),\n", " [20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30], transform=sst_data.metpy.cartopy_crs)\n", "ax.plot(track_x, track_y, marker='o', color='tab:blue', transform=ccrs.Geodetic())\n", "ax.add_feature(cfeature.COASTLINE.with_scale('10m'), edgecolor='orange')\n", "ax.add_feature(cfeature.STATES.with_scale('10m'), edgecolor='orange')\n", "ax.set_extent([-85, -69, 27, 38], crs=ccrs.PlateCarree())\n", "\n", "ax.set_title('GOES-16 Visible, SST (contours), NDBC Buoy Observations, NHC Best and Forecast Track')\n", "add_unidata_logo(fig, y=1375, x=2350, size='large')\n", "add_metpy_logo(fig, y=1375, x=25, size='large')\n", "add_timestamp(ax, high_contrast=True, y=0.01)\n", "plt.savefig('florence.png', dpi=200, bbox_inches='tight')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "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.6.4" } }, "nbformat": 4, "nbformat_minor": 2 }