GOES_aircraft

GOES Aircraft Example

This example plots the most recent GOES-16 data with an overlay of the current position of the NCAR C130 research aircraft. It also demonstrates projecting vectors into map coordinates and matplotlib marker manipulation.

In [1]:
from datetime import datetime
import json
import urllib.request
In [2]:
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from matplotlib import patheffects
import matplotlib.pyplot as plt
import metpy  # noqa: F401
import metpy.calc as mpcalc
from metpy.plots.ctables import registry
from metpy.units import units
import numpy as np
from siphon.catalog import TDSCatalog
import xarray as xr
from xarray.backends import NetCDF4DataStore
In [3]:
def get_plane_data():
    """Get JSON data from NCAR aircraft."""
    endpoint_url = 'https://www.eol.ucar.edu/flight_data/C130/position.json'
    with urllib.request.urlopen(endpoint_url) as f:
        jstring = f.read()
    payload = json.loads(jstring.decode('utf-8'))
    data = {'latitude': float(payload['lat']),
            'longitude': float(payload['lon']),
            'altitude': float(payload['alt']),
            'heading': float(payload['head']),
            'time': payload['timestamp']}
    return data
In [4]:
def get_goes_image(date=datetime.utcnow(), channel=8, region='CONUS'):
    """Return dataset of GOES-16 data."""
    cat = TDSCatalog('https://thredds.ucar.edu/thredds/catalog/satellite/goes/east/products/'
                     'CloudAndMoistureImagery/{}/Channel{:02d}/{:%Y%m%d}/'
                     'catalog.xml'.format(region, channel, date))

    ds = cat.datasets[-1]  # Get most recent dataset
    ds = ds.remote_access(service='OPENDAP')
    ds = NetCDF4DataStore(ds)
    ds = xr.open_dataset(ds)
    return ds
In [5]:
ds = get_goes_image()
data = get_plane_data()
In [6]:
# Parse out the projection data from the satellite file
dat = ds.metpy.parse_cf('Sectorized_CMI')
proj = dat.metpy.cartopy_crs
In [7]:
# Pull out what we need from the GOES netCDF file
x = dat['x']
y = dat['y']
In [8]:
# Make the plot
fig = plt.figure(figsize=(1.375 * 40, 40))
ax = fig.add_subplot(1, 1, 1, projection=proj)
plt.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0, hspace=0)

wv_norm, wv_cmap = registry.get_with_range('WVCIMSS_r', 195, 265)

im = ax.imshow(dat, extent=(x.min(), x.max(), y.min(), y.max()),
               origin='upper')

im.set_cmap(wv_cmap)
im.set_norm(wv_norm)

ax.add_feature(cfeature.BORDERS, linewidth=8, edgecolor='black')
ax.add_feature(cfeature.STATES.with_scale('50m'), linestyle='-',
               edgecolor='black', linewidth=4)

timestamp = datetime.strptime(ds.start_date_time, '%Y%j%H%M%S')

text_time = ax.text(0.01, 0.01, timestamp.strftime('%d %B %Y %H%MZ'),
                    horizontalalignment='left', transform=ax.transAxes,
                    color='white', fontsize=100, weight='bold')

outline_effect = [patheffects.withStroke(linewidth=15, foreground='black')]
text_time.set_path_effects(outline_effect)

ax.set_extent([-124.5, -105, 38.5, 50])

# Transform plane heading to a map direction and plot a rotated marker
u, v = mpcalc.wind_components(1 * units('m/s'),
                              data['heading'] * units('degrees'))
u, v = proj.transform_vectors(ccrs.PlateCarree(), np.array([data['longitude']]),
                              np.array([data['latitude']]), np.array([u.m]),
                              np.array([v.m]))
map_direction = -mpcalc.wind_direction(u * units.knots, v * units.knots).to('degrees')
map_direction = map_direction[0].m

ax.scatter(data['longitude'], data['latitude'],
           transform=ccrs.PlateCarree(),
           marker=(3, 0, map_direction),
           color='red',
           s=4000)

ax.text(data['longitude'], data['latitude'] - 0.3,
        'Altitude: {}\nHeading: {}\nTime:{}'.format(data['altitude'],
                                                    data['heading'],
                                                    data['time']),
        transform=ccrs.PlateCarree(), fontsize=40,
        ha='center', va='top',
        bbox={'facecolor': 'white', 'edgecolor': 'black',
              'boxstyle': 'round,pad=0.5', 'alpha': 0.6})

ax.gridlines(linestyle=':', color='black', linewidth=2)
Out[8]:
<cartopy.mpl.gridliner.Gridliner at 0x7f54aa40ad50>