Color images consist of red, green, and blue components. We can utilize the red and blue channels of the GOES-16 satellite, along with the “veggie” band to create a color image of the Earth. GOES-16 does not have a sensor that strictly covers the green visible range, so the image would be too green. Currently lookup tables are being created to simulate a true green sensor on the platform. Until those are completed, we will use an approximation that does a good job of producing a near real-color image. We'll also implement a square-root brightness enhancement.
We’ll need to accomplish a few tasks - making an outline like this is a good way to start thinking of how your program can be broken up into functions that encapsulate different functionality.
In [7]:
from datetime import datetime
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
from matplotlib import patheffects
import metpy
from metpy.plots import add_metpy_logo
import numpy as np
from scipy import interpolate
from siphon.catalog import TDSCatalog
%matplotlib inline
In [8]:
def open_dataset(date, channel, region, 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.ucar.edu/thredds/catalog/satellite/goes16/GOES16/'
'{}/Channel{:02d}/{:%Y%m%d}/catalog.xml'.format(region, channel, date))
ds = cat.datasets[idx]
return ds.remote_access(service='OPENDAP', use_xarray=True)
In [9]:
def interpolate_red_channel(red_ds, resampled_ds):
"""
Interpolate the red data channel to the same grid as another channel.
"""
x_new = resampled_ds.x
y_new = resampled_ds.y
f = interpolate.interp2d(red_ds.x, red_ds.y,
red_ds, fill_value=0)
red_interpolated = f(x_new, y_new)
return x_new, y_new, red_interpolated[::-1]
In [10]:
def make_RGB_data(date, region, 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, region, idx)
blue_ds = open_dataset(date, blue_channel, region, idx)
green_ds = open_dataset(date, green_channel, region, idx)
red_dat = red_ds.metpy.parse_cf('Sectorized_CMI')
blue_dat = blue_ds.metpy.parse_cf('Sectorized_CMI')
green_dat = green_ds.metpy.parse_cf('Sectorized_CMI')
proj = green_dat.metpy.cartopy_crs
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
# Brightness Enhancement
red_data = np.sqrt(red_data)
blue_data = np.sqrt(blue_data)
green_data = np.sqrt(green_data)
# Make better fake green channel
green_data = green_data * 0.1 + blue_data * 0.45 + red_data * 0.45
rgb_data = np.dstack([red_data, green_data, blue_data])
x = green_ds.x
y = green_ds.y
time = datetime.strptime(green_ds.start_date_time, '%Y%j%H%M%S')
return time, x, y, proj_var, rgb_data
In [11]:
# Get the next to most recent image and make the RGB data array.
date = datetime.utcnow()
timestamp, x, y, proj_var, rgb_data = make_RGB_data(date, 'Mesoscale-1', -2)
In [ ]:
# Same as before, except we call imshow with our colormap and norm.
fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(1, 1, 1, projection=proj)
im = ax.imshow(rgb_data, extent=(x.min(), x.max(), y.min(), y.max()), origin='upper')
ax.coastlines(resolution='50m', color='black')
ax.add_feature(cfeature.STATES, linestyle=':', edgecolor='black')
ax.add_feature(cfeature.BORDERS, linewidth=2, edgecolor='black')
# 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)
# Add the MetPy Logo
fig = add_metpy_logo(fig)