Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
GOES-16 RGB Animation
from datetime import datetime
from scipy import interpolate
from siphon.catalog import TDSCatalog
from netCDF4 import Dataset
import matplotlib.pyplot as plt
from matplotlib import patheffects
import cartopy.feature as cfeat
import cartopy.crs as ccrs
import numpy as np
from metpy.plots import add_logo
from matplotlib.animation import ArtistAnimation
def open_dataset(date, channel, idx):
"""
Open and return a netCDF Dataset object for a given date, channel, and image index
of GOES-16 data from THREDDS test server.
"""
cat = TDSCatalog('http://thredds-test.unidata.ucar.edu/thredds/catalog/satellite/goes16/GOES16/'
'Mesoscale-1/Channel{:02d}/{}/catalog.xml'.format(channel, date)) #Mesoscale-1
dataset = cat.datasets[idx]
ds = Dataset(dataset.access_urls['OPENDAP'])
return ds
def interpolate_red_channel(red_ds, resampled_ds):
"""
Interpolate the red data channel to the same grid as another channel.
"""
x_new = resampled_ds.variables['x'][:]
y_new = resampled_ds.variables['y'][::-1]
f = interpolate.interp2d(red_ds.variables['x'][:], red_ds.variables['y'][::-1],
red_ds.variables['Sectorized_CMI'][::-1,], fill_value=0)
red_interpolated = f(x_new, y_new[::-1])
return x_new, y_new, red_interpolated
def make_RGB_data(date, idx):
"""
Make an RGB image array, returning the time, coordinates, projection, and data.
"""
red_channel = 2
green_channel = 3
blue_channel = 1
red_ds = open_dataset(date, red_channel, idx)
blue_ds = open_dataset(date, blue_channel, idx)
green_ds = open_dataset(date, green_channel, idx)
green_data = green_ds.variables['Sectorized_CMI'][:]
blue_data = blue_ds.variables['Sectorized_CMI'][:]
x, y, red_data = interpolate_red_channel(red_ds, blue_ds)
# Clip out maxes
red_data[np.where(red_data <= 0.0001)] = 1
blue_data[np.where(blue_data <= 0.0001)] = 1
green_data[np.where(green_data <= 0.0001)] = 1
red_data = np.sqrt(red_data)
blue_data = np.sqrt(blue_data)
green_data = np.sqrt(green_data)
#red_data[np.where(red_data >= 0.8)] = 1
#blue_data[np.where(blue_data >= 0.8)] = 1
#green_data[np.where(green_data >= 0.8)] = 1
# Make better fake green channel
green_data = green_data*0.1 + blue_data*0.45 + red_data[::-1,:]*0.45
rgb_data = np.dstack([red_data[::-1,:], green_data, blue_data])
x = green_ds.variables['x'][:]
y = green_ds.variables['y'][:]
proj_var = green_ds.variables[green_ds.variables['Sectorized_CMI'].grid_mapping]
time = datetime.strptime(green_ds.start_date_time, '%Y%j%H%M%S')
return time, x, y, proj_var, rgb_data
# Get the most recent image and make the RGB data array.
date = datetime.strftime(datetime.utcnow(), '%Y%m%d')
timestamp, x, y, proj_var, rgb_data = make_RGB_data(date, -2)
# Create a Globe specifying a spherical earth with the correct radius
globe = ccrs.Globe(ellipse='sphere', semimajor_axis=proj_var.semi_major,
semiminor_axis=proj_var.semi_minor)
print(proj_var)
proj = ccrs.LambertConformal(central_longitude=proj_var.longitude_of_central_meridian,
central_latitude=proj_var.latitude_of_projection_origin,
standard_parallels=[proj_var.standard_parallel],
globe=globe)
# Set up a feature for the state/province lines. Tell cartopy not to fill in the polygons
state_boundaries = cfeat.NaturalEarthFeature(category='cultural',
name='admin_1_states_provinces_lakes',
scale='50m', facecolor='none')
artists = []
# Same as before, except we call imshow with our colormap and norm.
fig = plt.figure(figsize=(18, 15.7)) #18,20
ax = fig.add_subplot(1, 1, 1, projection=proj)
plt.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0, hspace=0)
ax.coastlines(resolution='50m', color='black')
ax.add_feature(state_boundaries, linestyle=':', edgecolor='black')
ax.add_feature(cfeat.BORDERS, linewidth='2', edgecolor='black')
add_logo(fig, size='large')
#-581
for idx in range(-291, -1):
print(idx)
# Get the most recent image and make the RGB data array.
date = datetime.strftime(datetime.utcnow(), '%Y%m%d')
timestamp, x, y, proj_var, rgb_data = make_RGB_data(date, idx)
im = ax.imshow(rgb_data, extent=(x.min(), x.max(), y.min(), y.max()), origin='upper')
# Add text (aligned to the right); save the returned object so we can manipulate it.
text_time = ax.text(0.99, 0.01, timestamp.strftime('%d %B %Y %H%MZ'),
horizontalalignment='right', transform=ax.transAxes,
color='white', fontsize='x-large', weight='bold')
text_channel = ax.text(0.5, 0.97, 'Experimental GOES-16 RGB Composite',
horizontalalignment='center', transform=ax.transAxes,
color='white', fontsize='large', weight='bold')
# Make the text stand out even better using matplotlib's path effects
outline_effect = [patheffects.withStroke(linewidth=2, foreground='black')]
text_time.set_path_effects(outline_effect)
text_channel.set_path_effects(outline_effect)
artists.append((im, text_time, text_channel))
# Create the animation--in addition to the required args, we also state that each
# frame should last 200 milliseconds
anim = ArtistAnimation(fig, artists, interval=50., blit=False)
anim.save('Irma_RGB.mp4')
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment