GOES CIRA Product Writer

Notebook Python-AWIPS Tutorial Notebook


  • Use python-awips to connect to an EDEX server

  • Define and filter the data request specifically for new CIRA GOES16 data products

  • Resize the products to their native resolution

  • Write the individual bands (channels) locally

  • Combine and write the RGB product locally

Table of Contents

1 Imports

The imports below are used throughout the notebook. Note the first import is coming directly from python-awips and allows us to connect to an EDEX server. The subsequent imports are for data manipulation and visualization.

from awips.dataaccess import DataAccessLayer
import cartopy.crs as ccrs
import cartopy.feature as cfeat
import matplotlib.pyplot as plt
from datetime import datetime
import numpy as np
import os


2 Initial Setup

2.1 EDEX Connection

First we establish a connection to Unidata’s public EDEX server. With that connection made, we can create a new data request object and set the data type to satellite.

# Create an EDEX data request
request = DataAccessLayer.newDataRequest()

2.2 Parameter Definition

After establishing the python-awips specific objects, we create a few other parameters that will be used for the data query based off of known values: projection, and extent.

# Create a projection for ECONUS and WCONUS
# Set up the projection using known parameters (from the netcdf of GOES products)
globe = ccrs.Globe(semimajor_axis=6378137.0, semiminor_axis=6356752.5, ellipse=None)
sat_h = 35785830.0
proj = ccrs.Geostationary(globe=globe, central_longitude=-75.0, satellite_height=sat_h, sweep_axis='x')

# Define the extents for ECONUS and WCONUS in goes native coords
# (originally taken from netcdf GOES data)
extent = (-3626751., 1382263.5, 1583666.1, 4588674.)


3 Function: set_size()

Here we’re defining a function that will allow us to pass in the dimensions of the output file we desire in pixels. Default Python methods require the size to be set in inches, which is confusing in our case, since we know what the size of GOES images are in pixels. Also, default Python functions add a padding when creating figures, and we don’t want that.

This function allows the exact final image to be specified based in pixels, with no padding or buffers.

def set_size(w,h, plt):
    """ w, h: width, height in pixels """

    # Convert from pixels to inches
    DPI = plt.figure().get_dpi()
    w = w/float(DPI)
    h = h/float(DPI)

    # Get the axes

    # Remove the padding
    l = ax.figure.subplotpars.left
    r = ax.figure.subplotpars.right
    t = ax.figure.subplotpars.top
    b = ax.figure.subplotpars.bottom
    figw = float(w)/(r-l)
    figh = float(h)/(t-b)

    # Set the final size
    ax.figure.set_size_inches(figw, figh)

    # Return the DPI, this is used when in the
    # write_image() function
    return DPI


4 Function: write_img()

Next, we’re defining another function which takes the image data, file name, projection, extent, reference time, and whether or not to print out a footnote.

This method specifies the size of the output image and creates a plot object to draw all our data into. Then it draws the GOES data, coastlines, state boundaries, and lat/lon lines onto the image. Additionally, if we want, it writes out a short footnote describing what product we’re looking at. Finally, it writes out the figure to disk.

By default we’re specifying the output dimensions to be 5000x4000 pixels, because that is the native GOES image size, but feel free to modify these values if you wish to print out an image of another size (you may want to keep the w:h ratio the same though).

def write_img(data, name, proj, extent, reftime, footnote):

    # Specify the desired size, in pixels
    px_width = 5000.0
    px_height = 3000.0

    # Create the plot with proper projection, and set the figure size
    fig = plt.figure()
    DPI = set_size(px_width, px_height, plt)
    ax = plt.axes(projection=proj)

    # Draw GOES data
    ax.imshow(data, cmap='gray', transform=proj, extent=extent)

    # Add Coastlines and States
    ax.coastlines(resolution='50m', color='magenta', linewidth=1.0)
    ax.add_feature(cfeat.STATES, edgecolor='magenta', linewidth=1.0)
    ax.gridlines(color='cyan', linewidth=2.0, xlocs=np.arange(-180, 180, 10), linestyle=(0,(5,10)))

    # Create and draw the footnote if needed
    if footnote:
        footnoteStr = '   CIRA-'+name[7:-4]+'-'+str(reftime)
        plt.annotate(str(footnoteStr), (0,0), (0, 0), xycoords='axes fraction', textcoords='offset points', va='top')

    # Write out the figure
    plt.savefig(name, dpi=DPI, bbox_inches='tight', pad_inches=0)


5 Get the Data and Write it Out!

5.1 Filter the Data

Define exactly what data we want to be printing out. This notebook is designed to loop through and print out multiple images, so here we can pick which images we’re wanting to print out. We’re specifying ECONUS (for East CONUS), CLDSNOW, DBRDUST, and GEOCOLR (for the new CIRA products) and the three channels for the RBG composites.

Tip: More information could be gathered by looking at all the available location names (sectors), identifiers (entities), and parameters (channels). To see those run the following lines of code after the dataType has been set to satellite on the request object:

## Print Available Location Names

## Print Available Identifiers and Values
ids = DataAccessLayer.getOptionalIdentifiers(request)
for id in ids:
    print(id, DataAccessLayer.getIdentifierValues(request, id))
# Define Location names
sectors = ["ECONUS"]

# Define creatingEntity Identifiers
entities = ["CLDSNOW", "DBRDUST", "GEOCOLR"]

# Define parameters
ch1 = "CH-01-0.47um"
ch2 = "CH-02-0.64um"
ch3 = "CH-03-0.87um"
channels = [ch1, ch2, ch3]

5.2 Define Output Location

Here we define a folder for where the satellite images will be written to. The default directory is a new folder called ‘output’ that lives whereever this notebook lives.

Tip: If you specify the fully qualified path, it will no longer depend on where this notebook is located. For example (for a Mac):

outputDir = '/Users/scarter/test_dir/output/'
# Define name of the desired end directory
outputDir = 'output/'

# Check to see if this folder exists
if not os.path.exists(outputDir):
    # If not, create the directory
    print('Creating new output directory: ',outputDir)
    # If so, let the user know
    print('Output directory exists!')
Output directory exists!

5.3 Write Out GOES Images

# First loop through the sectors (location names)
for sector in sectors:

    # Set the location on the request

    # Next loop through the Products (entities)
    for entity in entities:
        # Reset the time and channel variables since we're on a new product
        time = None
        R = None
        G = None
        B = None

        # Set the product
        request.addIdentifier("creatingEntity", entity)

        # Cycle through the channels (parameters)
        for channel in channels:

            # Set the time for this product if it hasn't been set
            # If it has been set, then we proceed with that value
            # so that all bands in for the one product are pulled
            # from the same time
            if(time is None):
                times = DataAccessLayer.getAvailableTimes(request)
                time = [times[-1]]
                print("selected time:", time)

            # Request the data from EDEX
            response = DataAccessLayer.getGridData(request, time)
            # Grab the actual data from the response
            grid = response[0]
            # Get the raw data from the response
            data = grid.getRawData()
            reftime = grid.getDataTime().getRefTime()

            # Set the R,G,B channel
            if(channel == ch1):
                B = data
            elif (channel == ch2):
                R = data
            elif (channel == ch3):
                G = data

            # Create the single channel name
            name = outputDir+entity+'-'+sector+'-'+channel+'.png'

            # Write out the single channel
            write_img(data, name, proj, extent, reftime, False)

        # --- End of channel loop

        # Create the RGB product
        # Apply range limits for each channel. RGB values must be between 0 and 1
        R = np.clip(R, 0, 1)
        G = np.clip(G, 0, 1)
        B = np.clip(B, 0, 1)
        RGB = np.dstack([R, G, B])

        # Create RGB name
        rgbName = outputDir+entity+'-'+sector+'-RGB.png'
        # Write out the RGB image
        print('writing', rgbName)
        write_img(RGB, rgbName, proj, extent, time, False)

    # --- End of entity loop
#--- End of sector loop

selected time: [<DataTime instance: 2021-05-28 06:51:14 >]
writing output/CLDSNOW-ECONUS-CH-01-0.47um.png
writing output/CLDSNOW-ECONUS-CH-02-0.64um.png
writing output/CLDSNOW-ECONUS-CH-03-0.87um.png
writing output/CLDSNOW-ECONUS-RGB.png
selected time: [<DataTime instance: 2021-05-28 06:51:14 >]
writing output/DBRDUST-ECONUS-CH-01-0.47um.png
writing output/DBRDUST-ECONUS-CH-02-0.64um.png
writing output/DBRDUST-ECONUS-CH-03-0.87um.png
writing output/DBRDUST-ECONUS-RGB.png
selected time: [<DataTime instance: 2021-05-28 06:56:14 >]
writing output/GEOCOLR-ECONUS-CH-01-0.47um.png
writing output/GEOCOLR-ECONUS-CH-02-0.64um.png
writing output/GEOCOLR-ECONUS-CH-03-0.87um.png
writing output/GEOCOLR-ECONUS-RGB.png
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>


6 See Also

6.2 Additional Documentation

CIRA Quick Guides