======================== GOES CIRA Product Writer ======================== `Notebook `_ Python-AWIPS Tutorial Notebook -------------- Objectives ========== - 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 `__\ | `2 Initial Setup `__\ |     `2.1 EDEX Connection `__\ |     `2.2 Parameter Definition `__\ | `3 Function: set_size() `__\ | `4 Function: write_img() `__\ | `5 Get the Data and Write it Out! `__\ |     `5.1 Filter the Data `__\ |     `5.2 Define Output Location `__\ |     `5.3 Write Out GOES Images `__\ | `6 See Also `__\ |     `6.1 Related Notebooks `__\ |     `6.2 Additional Documentation `__\ -------------- 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. .. code:: ipython3 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 `Top `__ -------------- 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**. .. code:: ipython3 # Create an EDEX data request DataAccessLayer.changeEDEXHost("edex-cloud.unidata.ucar.edu") request = DataAccessLayer.newDataRequest() request.setDatatype("satellite") 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. .. code:: ipython3 # 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.) `Top `__ -------------- 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. .. code:: ipython3 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 ax=plt.gca() # 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 `Top `__ -------------- 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). .. code:: ipython3 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) `Top `__ -------------- 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((DataAccessLayer.getAvailableLocationNames(request)) ## Print Available Identifiers and Values ids = DataAccessLayer.getOptionalIdentifiers(request) print(ids) for id in ids: print(id, DataAccessLayer.getIdentifierValues(request, id)) .. code:: ipython3 # 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/' .. code:: ipython3 # 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) os.makedirs(outputDir) else: # If so, let the user know print('Output directory exists!') .. parsed-literal:: Output directory exists! 5.3 Write Out GOES Images ~~~~~~~~~~~~~~~~~~~~~~~~~ .. code:: ipython3 # First loop through the sectors (location names) for sector in sectors: # Set the location on the request request.setLocationNames(sector) # 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: request.setParameters(channel) # 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 print('writing',name) 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 print('Done!') .. parsed-literal:: selected time: [] 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: [] 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: [] 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 Done! .. parsed-literal::
.. image:: GOES_CIRA_Product_Writer_files/GOES_CIRA_Product_Writer_25_2.png .. parsed-literal::
.. image:: GOES_CIRA_Product_Writer_files/GOES_CIRA_Product_Writer_25_4.png .. parsed-literal::
.. image:: GOES_CIRA_Product_Writer_files/GOES_CIRA_Product_Writer_25_6.png `Top `__ -------------- 6 See Also ---------- 6.1 Related Notebooks ~~~~~~~~~~~~~~~~~~~~~ - `Satellite Imagery `__ 6.2 Additional Documentation ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ **CIRA Quick Guides** - `DEBRA-Dust `__ - `Cloud-Snow `__ - `GEOCOLOR `__ **python-awips** - `DataAccessLayer.changeEDEXHost() `__ - `DataAccessLayer.newDataRequest() `__ - `DataAccessLayer.getAvailableLocationNames() `__ - `DataAccessLayer.getOptionalIdentifiers() `__ - `DataAccessLayer.getIdentifierValues() `__ - `DataAccessLayer.getAvailableTimes() `__ - `IDataRequest `__ **matplotlib** - `matplotlib.pyplot() `__ - `matplotlib.pyplot.axes() `__ - `matplotlib.pyplot.figure() `__ **numpy** - `numpy.clip() `__ - `numpy.dstack() `__ `Top `__ --------------