GOES-16 RGB Animation

rgb_animation.py

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')