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