Python Hub

Home Page Anaconda Download Tutorial ABI CONUS AOD Tutorial VIIRS Active Fires Tutorial
Block 1 Block 2 Block 3 Block 4 Block 5 Block 6 Block 7 Block 8 Block 9 Block 10 Block 11 Block 12 Block 13 Block 14

Processing and Visualizing GOES-16/ABI CONUS AOD using Python Jupyter Notebooks
Instructions and Detailed Code Annotation

Workshop on New Generation Satellite Products for Operational Fire and Smoke Applications
3rd International Smoke Symposium
April 20, 2020


The Advanced Baseline Imager (ABI) on the GOES-16 (GOES-East) and GOES-17 GOES-West) satellites is NOAA’s next generation instrument for imaging Earth’s weather, oceans, and environment. ABI has 16 spectral bands: 2 visible bands, 4 near-infrared bands, and 10 infrared bands. Observations from ABI’s 16 bands can be used individually or in combination (in scientific algorithms) to measure different aspects of the Earth’s surface (e.g., vegetation, fires) and atmosphere (e.g., clouds, moisture, aerosols). Aerosol optical depth (AOD) is an ABI product that quantitatively measures aerosols in the atmosphere from such sources as smoke plumes, blowing dust, and urban haze.

In scan mode 6 (the operational scan mode as of April 2, 2019), ABI observes the continental U.S. (CONUS) every 5 minutes and the full disk of the hemisphere every 10 minutes. The ABI AOD product has 2 km spatial resolution in both views.

The Python Jupyter Notebook code (.ipynb) documented here was written for the Workshop on New Generation Satellite Products for Operational Fire and Smoke Applications, held at the Third International Smoke Symposium (ISS3) on April 20, 2020. The code shows users how to open and explore a netCDF-4 file containing GOES-16 CONUS view ABI AOD data; how to process AOD data and data quality flags (DQFs) for a smoke event on May 30, 2019; how to plot AOD data on maps to create professional-looking figures; and finally how to create an animation of multiple AOD figures, representing 1 hour of ABI AOD observations.

To download the 12 netCDF-4 ABI CONUS view AOD data files used in this tutorial click here. If you don’t have Jupyter Notebook on your computer, follow the instructions here to download Anaconda and install Jupyter Notebook.


Block #1: Import Libraries and Settings

To make things simple in subsequent code blocks, first we import all of the libraries and settings we will need to read, process, and visualize the ABI AOD data. We need libraries to do things like perform array operations (numpy), plot data (matplotlib), make maps (cartopy), and read netCDF-4 files (netCDF4, Dataset).

#Block 1: Import libraries and settings

#Library to perform array operations
import numpy as np 

#Libraries for making plots
import matplotlib as mpl
from matplotlib import pyplot as plt
import matplotlib.ticker as ticker

#Libaries for drawing maps
import cartopy
from cartopy import crs as ccrs
import cartopy.feature as cfeature
from cartopy.feature import NaturalEarthFeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter

#Library for creating animations
from PIL import Image

#Library for accessing files in the directory
import os

#Library to read in netCDF files
from netCDF4 import Dataset

#Library for using math functions
import math

#Library for collecting lists of files from folders
import glob

import warnings
warnings.filterwarnings('ignore')

#Sets font size to 12
plt.rcParams.update({'font.size': 12})

#Option to keep numpy from printing in scientific notation by default
np.set_printoptions(suppress = True)

Block #2: Explore an ABI AOD Data File

Let’s open a netCDF-4 data file containing ABI AOD data. The data file names are very long!

#Block 2: Explore an ABI AOD data file

#Enter file name
fname = os.getcwd() + '/OR_ABI-L2-AODC-M6_G16_s20191501801401_e20191501804174_c20191501806221.nc'

#Set the file name to read
file_id = Dataset(fname)
The graphic below shows how to decipher the data file name:
  • “OR” represents the system environment that generated the file; in this example, the file was created by the Operational system using Real-time data
  • “ABI” represents the sensor (instrument); in this example, the Advanced Baseline Imager
  • “L2” represents the data processing level; in this example, it indicates Level 2 data
  • “AODC” refers to the product name; in this example, AOD data, CONUS view
  • “M6” represents the satellite scan mode; in this example, ABI scan Mode 6: CONUS view observations every 5 min and full disk view observations every 10 min
  • “G16” represents the satellite; in this example, the GOES-16 satellite
  • For the observation Start (“s”) and End (“e”) times, and the file Creation (“c”) time:
    • “YYYY” is the year (e.g., 2019)
    • “DDD” is the Julian day (e.g., 150)
    • “HH” and “MM” are the hour and minute in UTC time (e.g., 18:01 UTC)
    • “SSS” are the seconds to the tenth of a second
  • “.nc” indicates the file is a netCDF file

ABI File Names
Some Notes about netCDF-4 Files
NOAA satellite data are distributed in netCDF-4 file format. netCDF files have a common organizational structure, illustrated in the graphic below. Files contain one top-level group (“root group”); content is organized by optional sub-groups. Groups have attributes (descriptive information) and dimensions. Data are organized as variables, which also have attributes and dimensions. For example, there are >70 variables in an ABI AOD netCDF-4 file. Attributes and variables each have a data type. Examples of common data types you will encounter when working with satellite netCDF-4 files include: string (str), byte, integer (int), and floating point number (float).

nc4 model

Image courtesy of University Corporation for Atmospheric Research: https://www.unidata.ucar.edu/software/netcdf/docs/netcdf_data_model.html



In Block #2, uncomment the “print” statements (remove the leading ##) one by one and re-run the block to explore the contents of the data file. The output of each “print” statement is too long to replicate here, so look at the output in your Jupyter Notebook. We give some more information and key details about the data variables below.

The metadata for the entire file shows the file contents and structure, including the root group, attributes, dimensions, and variables. The variables contain the data we are interested in: AOD, DQF, x (latitude in radians), y (longitude in radians).
#Explore the contents of the file by UNCOMMENTING the 'print' statements one by one to see various aspects- 
#of the file
#Check the contents of the entire file print(file_id)

The metadata for the AOD variable shows us the attributes of the data. Because the AOD data are 2-dimensional (2D), we add [: , :] to the end of the variable file path name to display a snippet of the AOD array values. The symbol “ - -“ in the array indicates a pixel for which no AOD was retrieved. Before we go any further, there are a couple of very important things to notice about the AOD data.
  1. AOD data are stored in the netCDF-4 file as 16-bit unsigned integers (int16 data type) to save file space. To convert the data from integers to floating point numbers, a scale factor (“scale_factor”) and offset (“add_offset”) are included in the metadata; multiply AOD by the scale and then add the offset. Python automatically applies this conversion to the AOD data! We will check the data types in Block #5 to confirm the automatic conversion.
  2. The fact that Python automatically converts the AOD data to floats results in an error in the displayed valid AOD data range; the actual range is [-0.05, 5] but it is shown as [0, -6] in the metadata.

#Check the AOD variable metadata
##print(file_id.variables['AOD'])

#Check the AOD array values
##print(file_id.variables['AOD'][:,:])

The metadata for the DQF variable show us that the variable is an integer with four possible values:

  • 0: high quality AOD
  • 1: medium quality AOD
  • 2: low quality AOD
  • 3: no AOD retrieval

The graphic below shows an example of the data quality flags for ABI AOD data, full disk view, on 31 July 2017 at 15:45 UTC, highlighting the areas of high quality AOD, medium quality AOD, low quality AOD, and areas where AOD was not retrieved. The reasons for designation of low quality and not retrieved AOD are indicated. Notice that AOD is not retrieved for areas where there are clouds, snow, or bright land surfaces. The most striking area where AOD is not retrieved is the sun glint area over the ocean (dark brown circular area in the graphic). Sun glint is a measurement artifact that occurs when sunlight reflects off the surface of the ocean at the same angle that a satellite sensor is viewing the surface. ABI Quality Flags

Similar to the AOD variable, the DQF variable is 2D, so we add [: , :] to the end of the variable file path name to display a snippet of the DQF array values. As with the AOD data array, the symbol “ - -“ in the array indicates a pixel for which no AOD was retrieved and thus, for which there is no DQF.

#Check the DQF variable metadata
##print(file_id.variables['DQF'])

#Check the DQF array values
##print(file_id.variables['DQF'][:,:])

Block #3: Check the Spatial Resolution of the Product

It’s always useful to check the spatial resolution of the data with which you are working. As mentioned previously, ABI AOD data have 2 km resolution at nadir (when the satellite view is looking straight down (0°)).
#Block 3: Check the spatial resolution of the data

print((file_id.__getattr__('title')),'spatial resolution is', (file_id.__getattr__('spatial_resolution')))

Block #4: Check the Units for the Variables of Interest

We strongly recommend checking the file for the units of the variables with which you are working. Occasionally, the units will change (e.g., km2 to m2) when a product is updated and the version changes; if you miss the change, your processed data will be off by a corresponding amount. In this case, AOD and DQF are unitless. ABI files always contain latitude (“x”) and longitude (“y”) in radians, in order to save file space. We will convert from radians to degrees in Block #6.
#Block 4: Check the units for the variables of interest (note: "1" means unitless)

print('AOD unit is', (file_id.variables['AOD'].__getattr__('units')))
print('DQF unit is', (file_id.variables['DQF'].__getattr__('units')))
print('Latitude unit is', (file_id.variables['x'].__getattr__('units')))
print('Longitude unit is', (file_id.variables['y'].__getattr__('units')))

Block #5: Check the Data Types for the Variables of Interest

We strongly recommend checking the data types for the variables with which you will be working. As mentioned previously, Python automatically applies the scale and offset to the AOD variable to convert from an integer (int) to floating point number (float). Be aware that if you use a program written in a different language, such as IDL, you may have to manually convert AOD using the scale and offset! In this case, we see that the AOD, latitude, and longitude are floats, and DQF is an integer (which we know it should be from looking at its metadata in Block #2), so we can proceed.
#Block 5: Check the data types for the variables of interest

print('AOD data type is', (file_id.variables['AOD'][:,:].dtype))
print('DQF data type is', (file_id.variables['DQF'][:,:].dtype))
print('Latitude data type is', (file_id.variables['x'][:].dtype))
print('Longitude data type is', (file_id.variables['y'][:].dtype))

Block #6: Algorithm to Convert Latitude and Longitude Radian Values to Degrees

Because latitude (“x”) and longitude (“y”) are stored in ABI files in radians, to save file space, we need to convert them to degrees before we can plot the AOD data on a map. More information about the theory and mathematics behind this conversion is available here.

We use the “def” command to create a function called “Degrees(file_id)” that reads in the latitude and longitude variables in radians, converts from radians to degrees, and returns latitude (“Lat”) and longitude (“Lon”) in degrees.

Functions allow us to complete a specific action like converting latitude and longitude in radians to degrees over and over again, without having to re-write the full code each time. We will create more functions in Blocks #7, #9, #10, and #11.
#Block 6: Algorithm to convert latitude and longitude radian values to degrees (as a callable function)

def Degrees(file_id):
    proj_info = file_id.variables['goes_imager_projection']
    lon_origin = proj_info.longitude_of_projection_origin
    H = proj_info.perspective_point_height+proj_info.semi_major_axis
    r_eq = proj_info.semi_major_axis
    r_pol = proj_info.semi_minor_axis
    
    #Data info
    lat_rad_1d = file_id.variables['x'][:]
    lon_rad_1d = file_id.variables['y'][:]
    
    #Create meshgrid filled with radian angles
    lat_rad,lon_rad = np.meshgrid(lat_rad_1d,lon_rad_1d)
    
    #lat/lon calculus routine from satellite radian angle vectors
    lambda_0 = (lon_origin*np.pi)/180.0
    
    a_var = np.power(np.sin(lat_rad),2.0) + (np.power(np.cos(lat_rad),2.0)*(np.power(np.cos(lon_rad),2.0)
+
(((r_eq*r_eq)/(r_pol*r_pol))*np.power(np.sin(lon_rad),2.0)))) b_var = -2.0*H*np.cos(lat_rad)*np.cos(lon_rad) c_var = (H**2.0)-(r_eq**2.0) r_s = (-1.0*b_var - np.sqrt((b_var**2)-(4.0*a_var*c_var)))/(2.0*a_var) s_x = r_s*np.cos(lat_rad)*np.cos(lon_rad) s_y = - r_s*np.sin(lat_rad) s_z = r_s*np.cos(lat_rad)*np.sin(lon_rad) Lat = (180.0/np.pi)*(np.arctan(((r_eq*r_eq)/(r_pol*r_pol))*((s_z/np.sqrt(((H-s_x)*(H-s_x))+(s_y*s_y)))))) Lon = (lambda_0 - np.arctan(s_y/(H-s_x)))*(180.0/np.pi) return Lat, Lon

Block #7: Select and Process AOD Data from a Single File

Here, we create a function (“AOD_Data(file_id)”) to read in the AOD variable and use the DQF variable to select the desired AOD data quality. Recall from Block #2 that the DQF variable designates AOD pixels as high quality (DQF = 0), medium quality (DQF = 1), low quality (DQF = 2), and not retrieved (DQF = 3). You may chose different AOD quality depending on your application; the NOAA Satellite Aerosol Team recommends:
  • For operational applications, use high and medium quality (“top 2 qualities”)
  • For science applications, like data assimilation, use high quality only
In this case, we select high and medium quality AOD by masking the low quality pixels and pixels for which AOD was not retrieved (DQF >1). The function returns only pixels with high or medium quality AOD.
#Block 7: Select and process AOD data from a single file (as a callable function)

def AOD_Data(file_id):
    #Read in AOD data
    AOD_data = file_id.variables['AOD'][:,:]

    #Select quality of AOD data pixels using the "DQF" variable
    #High quality: DQF = 0, Medium quality: DQF = 1, Low quality: DQF = 2, not retrieved (NR): DQF = 3
    #Science team recommends using High and Medium qualities for operational applications- 
#(e.g.,mask low quality and NR pixels)
DQF = file_id.variables['DQF'][:,:] Quality_Mask = (DQF > 1) AOD = np.ma.masked_where(Quality_Mask, AOD_data) return AOD

Block #8: Review Processed Data

Before attempting to plot satellite data, we recommend a quick check to see if you have data (i.e., if the processing steps in Block #7 ran correctly) and if the data look reasonable. Note how we use our functions from Blocks #7-8 to read in and process the AOD, latitude, and longitude with only two lines of code.

In this case, the min/max AOD fit within the specified range seen in the metadata in Block #2 ([-0.05, 5], and the latitude/longitude values are in degrees in a range that makes sense for the CONUS.
#Block 8: Review processed data (sanity check)

AOD = AOD_Data(file_id)
Lat, Lon = Degrees(file_id)

print('AOD: minimum value is ' + str(np.min(AOD)) + ';' + ' maximum value is ' + str(np.max(AOD)))
print('Latitude: minimum value is ' + str(np.min(Lat)) + ' degrees;' + ' maximum value is ' 
+ str(np.max(Lat)) + ' degrees') print('Longitude: minimum value is ' + str(np.min(Lon)) + ' degrees;' + ' maximum value is '
+ str(np.max(Lon)) + ' degrees')

Block #9: Plotting Settings for AOD Data

We create a function (“AOD_Data_Settings( )”) to return a couple of essential settings for plotting AOD data.

We create a custom continuous color bar for the AOD data using 11 colors from the matplotlib (abbreviated as “mpl”) library. A list of all of the available colors is found here. Typically, AOD data are plotted on a scale from 0 to 1, but as we know from the metadata, ABI AOD can be as high as 5. So we pick a different color (‘darkred’) for AOD data that are > 1.

We also need to set the range for plotting the AOD data using the numpy (abbreviated as “np”) “arrange” function. The data min is set as 0, the data max is set as 1.1 (in order to plot data up to a value of 1), and the contour interval is set as 0.05. Note that the smaller the contour interval, the higher the resolution of the plotted AOD data, but the longer it will take for the plotting code to run. We select 0.05 here as a compromise for relatively high resolution plotted data which doesn’t take too long to run. If you are plotting data from a large file, like Full Disk view ABI AOD, then select a coarser interval (e.g., 0.1). If you don’t care as much about how long it takes the code to run, but you want a really detailed figure, then select a finer interval (e.g., 0.01).
#Block 9: Plotting settings for AOD data (as a callable function)

def AOD_Data_Settings():
    #Create custom continuous colormap for AOD data
    #.set_over sets color for plotting data > max
    color_map = mpl.colors.LinearSegmentedColormap.from_list('custom_AOD', [(0, 'indigo'),(0.1, 'mediumblue'), 
(0.2, 'blue'), (0.3, 'royalblue'), (0.4, 'skyblue'), (0.5, 'cyan'), (0.6, 'yellow'), (0.7, 'orange'),
(0.8, 'darkorange'), (0.9, 'red'), (1, 'firebrick')], N = 150) color_map.set_over('darkred') #Set range for plotting AOD data (data min, data max, contour interval) (MODIFY contour interval) #interval: 0.1 = runs faster/coarser resolution, 0.01 = runs slower/higher resolution data_range = np.arange(0, 1.1, 0.05) return color_map, data_range

Block #10: Create AOD Colorbar, Independent of Plotted Data

We create a function (“AOD_Colorbar( )”) to display the AOD colorbar on the maps we will create of AOD data. In Block #9, we created a custom color map for the AOD data; in this block, we put the color map in a colorbar with labels and tick marks. The custom AOD color bar is shown below. This colorbar is special in a couple of ways. First, it is independent of the plotted data, meaning the colorbar itself won’t vary, even slightly, from figure to figure. Second, this colorbar scales automatically when the domain of the plotted area changes or the figure is re-sized. This saves us a lot of tedious adjustments when we get to the plotting steps. AOD Colorbar
#Block 10: Create AOD colorbar, independent of plotted data (as a callable function)
##cbar_ax are dummy variables
#Location/dimensions of colorbar set by .set_position (x0, y0, width, height) to scale automatically with plot

def AOD_Colorbar():
    last_axes = plt.gca()
    cbar_ax = fig.add_axes([0, 0, 0, 0])
    plt.draw()
    posn = ax.get_position()
    cbar_ax.set_position([0.35, posn.y0 - 0.07, 0.3, 0.02])
    color_map = mpl.colors.LinearSegmentedColormap.from_list('custom_AOD', [(0, 'indigo'),(0.1, 'mediumblue'), 
(0.2, 'blue'), (0.3, 'royalblue'), (0.4, 'skyblue'), (0.5, 'cyan'), (0.6, 'yellow'), (0.7, 'orange'),
(0.8, 'darkorange'), (0.9, 'red'), (1, 'firebrick')], N = 150) color_map.set_over('darkred') norm = mpl.colors.Normalize(vmin = 0, vmax = 1) cb = mpl.colorbar.ColorbarBase(cbar_ax, cmap = color_map, norm = norm, orientation = 'horizontal',
ticks = [0, 0.25, 0.5, 0.75, 1], extend = 'max') cb.set_label(label = 'AOD', size = 'medium', weight = 'bold') cb.ax.set_xticklabels(['0', '0.25', '0.50', '0.75', '1.0']) cb.ax.tick_params(labelsize = 'medium') plt.sca(last_axes)

Block #11: Format Map with Plate Carrée Projection

There are many steps in Python to create a professional quality map. Here, we create a function (“ABI_Map_Settings_PC(ax)”) to capture these formatting settings for a map with a Plate Carrée (PC) equidistant rectangular projection. To set up and label a latitude and longitude grid for our map, we use cartopy’s “LatitudeFormatter( )” and “LongitudeFormatter( )” functions. We’ve set latitude labels every 10 degrees and longitude labels every 20 degrees (for the Western Hemisphere only, since the GOES are geostationary satellites) as a general default; you can modify the labels by changing the values inside the “ax.set_xticks” and “ax.set_yticks” commands.
#Block 11: Format map with Plate Carree projection (as a callable function)

def ABI_Map_Settings_PC(ax):
    #Set up and label the lat/lon grid
    lon_formatter = LongitudeFormatter()
    lat_formatter = LatitudeFormatter()
    ax.xaxis.set_major_formatter(lon_formatter)
    ax.yaxis.set_major_formatter(lat_formatter)
    ax.set_xticks([-160, -140, -120, -100, -80, -60, -40, -20], crs = ccrs.PlateCarree())
    ax.set_yticks([-80,-70,-60,-50,-40,-30,-20,-10,0,10,20,30,40,50,60,70,80], crs = ccrs.PlateCarree())

We also need gridlines for our lat/lon grid. We add them with “ax.grid( )”. The default gridline in matplotlib is a light grey solid line, which we like; if you want a different color, linestyle, or any other formatting, consult the arguments listed here. The “zorder” argument sets the drawing order for layers; layers with larger zorder numbers plot on top of layers with smaller zorder numbers. We like a clean plot with no lat/lon ticks, so we set “length = 0” for “ax.tick_params( )”. Modify the formatting using the arguments listed here.
#Set lat/lon ticks and gridlines
ax.tick_params(length = 0)
ax.grid(linewidth = 0.5, zorder = 3)

We use cartopy and the Natural Earth Feature shapefile database to add mapping features like coastlines, international borders, and state borders. We also add polygons for land, oceans, and lakes, using the “facecolor” argument to shade the land ‘grey’ and the oceans and lakes ‘lightgrey’. The zorder arguments ensure coastlines and borders plot on top of the land, ocean, and lakes. The “scale” parameter sets the spatial resolution of the Natural Earth Features.
#Draw coastlines/borders using Cartopy; zorder sets drawing order for layers
    ax.coastlines(resolution = '50m', zorder = 3)
    ax.add_feature(cfeature.BORDERS, zorder = 3)
    ax.add_feature(cfeature.NaturalEarthFeature(category = 'cultural', name = 'admin_1_states_provinces', 
scale = '50m'), facecolor = 'none', lw = 0.5, edgecolor = 'black', zorder = 2) ax.add_feature(cfeature.NaturalEarthFeature(category = 'physical', name = 'ocean', scale = '50m'),
facecolor = 'lightgrey') ax.add_feature(cfeature.NaturalEarthFeature(category = 'physical', name = 'land', scale = '50m'),
facecolor = 'grey') ax.add_feature(cfeature.NaturalEarthFeature(category = 'physical', name = 'lakes', scale = '50m'),
facecolor = 'lightgrey', edgecolor = 'black', zorder = 2)

We zoom into an area of interest using “ax.set_extent([x0, x1, y0, y1])”, where (x0, x1) are the longitude values and (y0, y1) are the latitude values, as shown in the graphic below. Note that negative values must be used to designate degrees West longitude when working with ABI data. It is not necessary to specify a coordinate system in “ax.set_extent([ ])” if the extent values are given in the same coordinate system as the axes’ projection; in this case, Plate Carrée. Note that if you comment out (add leading ##) the “ax.set_extent([ ]) line, the map domain will automatically extend to the limits of the data.
#Set domain for map [x0, x1, y0, y1] 
    #Default longitude extent (x0, x1) values: G16 = (-135, -65); G17 = (-170, -100)
    #Use 180 degrees for longitude coordinates (i.e, -100 = 100 degrees W)
    #NOTE: comment out (add leading ##) the line below to automatically set domain to extend to limits of data
    ax.set_extent([-135, -65, 15, 55], crs = ccrs.PlateCarree())
Lat Lon Extent

Block #12: Plot AOD Data from a Single File CONUS View

Now, we put together all the functions we created to make a plot of high and medium quality AOD from a single file (5 min of observations). We read in our processed AOD data and our latitude and longitude in degrees. We set up a figure and axes with the Plate Carrée projection; this is also the source projection for the ABI data. There are many other map projection options; examples are listed here. We use the functions we made previously to add the map settings, AOD colorbar, and plotting settings for the AOD data. We also add a title to our figure, using reverse indexing (from right to left) of the file name to automatically add the satellite number, time in UTC, and year to the title.

We create a filled contour plot of the AOD data using “ax.contourf( )”. The “extend = both” argument plots any AOD values that are both above and below the set range of [0, 1]. Recall that we set a color (‘darkred’) for AOD values > 1. Values < 0 are plotted the same color (‘indigo’) as values = 0, in order to de-emphasize them. Negative AOD values represent the uncertainty in the AOD measurement and can be considered physically as very small positive AOD.

In rare cases, there may be no high or medium quality AOD data in a file. If we try to plot AOD in this case, we will get an error message. To avoid this problem, we put an “if/else” statement around the “Plot” and use the “.count( )” function to check if the processed AOD array contains any numerical values. If it doesn’t, we skip to the next data file (“pass”).

Notice that in “Plot” we set the “transform” parameter to the Plate Carrée projection (“transform = ccrs.PlateCarree( )”). This parameter transforms the data being plotted from its source coordinates to the map projection. Since the ABI data we are plotting and the map we are plotting the data onto are both in the Plate Carrée projection, specifying the “transform” parameter in this case is not strictly necessary (if you delete it, the data will still plot correctly). Nevertheless, we recommend getting into the habit of always specifying the “transform” parameter in plotting functions. The reason is because if you ever chose a different map projection such as Orthographic, for example (ABI full disk view) and you don’t specify “transform = ccrs.PlateCarree( )”, the satellite data will not plot correctly.

We "show" the map of plotted data so we can see the output but wait until the next block to save it. We recommend plotting one data file first, like this, to check the map settings before moving onto processing and plotting AOD from multiple data files.
#Block 12: Plot AOD data from a single file - CONUS View

#Select and process AOD data
AOD = AOD_Data(file_id)

#Read in latitutude and longitude values in degrees
Lat, Lon = Degrees(file_id)

#Set up figure and map projection: PlateCarree(central_longitude)
#Plate Carree: equidistant cylindrical projection w/equator as the standard parallel;- 
#default central_longitude = 0
fig = plt.figure(figsize=(8, 10)) ax = fig.add_subplot(1,1,1, projection = ccrs.PlateCarree()) #Format map with Plate Carree projection ABI_Map_Settings_PC(ax) #Add and format title #Reverse indexing (from right to left) of file name automatically adds satellite, time, and year to title plt.title('GOES-' + fname[-53:-51] + '/ABI\nHigh + Medium Quality AOD\n' + fname[-42:-40] + ':' + fname[-40:-38]
+ ' UTC, 30 May ' +fname[-49:-45], y = 1.025, ma = 'center', size = 15, weight = 'bold') #Add AOD colorbar AOD_Colorbar() #Plotting settings for AOD data color_map, data_range = AOD_Data_Settings() if AOD.count() > 0: #Create filled contour plot of AOD data Plot = ax.contourf(Lon, Lat, AOD, data_range, cmap = color_map, extend = 'both', zorder = 3,
transform = ccrs.PlateCarree()) else: pass #Show figure plt.show() #Don't save this figure now - we will save it in the next step! But the code showing how to save is here-
#for reference:
#Save figure as a .png file #dpi sets the resolution of the digital image in dots per inch filename = 'G16_CONUS_ABI_AOD_20190530_' + fname[-42:-38] ##fig.savefig(filename, bbox_inches = 'tight', dpi = 150)
G16 CONUS ABI AOD 20190530 1801

Block #13: Make multiple individual figures of AOD data (one plot from each file) CONUS View

Given the high temporal resolution of the ABI, you will likely want to make multiple images at once, to see the evolution of an event. In this case, we will make and save images from 1 hour of ABI CONUS view AOD (12 data files). The content of this block is very similar to that of Block #12. The difference is that we add the “glob( )” function to retrieve all of the ABI AOD netCDF-4 (.nc) files in our “current working directory” and a “for” loop to cycle through the 12 data files. In this way, we can process and plot the data from each file.

We “show” each map of plotted data (so we can see the output) and save each graphic file using a unique filename that includes the observation start time in UTC, taken by reverse indexing (from right to left) the netCDF-4 file name.

In the “fig.savefig” function, setting “bbox_inches = ‘tight’” ensures the entire plot is saved, including the title. The “dpi” parameter sets the resolution of the image in dots per inch. We have set “dpi = 150”, which is moderate resolution. Here are some suggested guidelines for setting the dpi parameter, depending on the application of your graphic:
  • 100 (default): relatively low resolution, suitable for webpages
  • 150: moderate resolution, suitable for presentations
  • 300: high resolution, suitable for written reports
  • 600: very high resolution, suitable for scientific journal articles

Note that the higher the dpi value, the larger the resulting file size will be, and the longer the processing time. So we recommend setting dpi to ≤ 150 for routine graphics.

Find the saved figures in your “current working directory” in this tutorial, that is your folder on Jupyter Notebook where you saved the Python Notebook .ipynb file and netCDF-4 data files.
#Block 13: Make multiple individual figures of AOD data (one plot from each data file) - CONUS View

#Collect all of the AOD CONUS view .nc files in given subdirectory
file_list = sorted(glob.glob(os.getcwd() + '/*AODC*.nc'))

#Plotting settings for AOD data
color_map, data_range = AOD_Data_Settings()

#Loop through data files, making/saving a figure for each data file
for x in file_list:
    file_id = Dataset(x)  
    
    #Select and process AOD data
    AOD = AOD_Data(file_id)
    
    #Read in latitutude and longitude values in degrees
    Lat, Lon = Degrees(file_id)
    
    #Set up figure and map projection: PlateCarree(central_longitude)
    #Plate Carree: equidistant cylindrical projection w/equator as the standard parallel;- 
#default central_longitude = 0
fig = plt.figure(figsize=(8, 10)) ax = fig.add_subplot(1,1,1, projection = ccrs.PlateCarree()) #Format map with Plate Carree projection ABI_Map_Settings_PC(ax) #Add and format title #Reverse indexing (from right to left) of file name automatically adds satellite, time, and year to title plt.title('GOES-16/ABI\nHigh + Medium Quality AOD\n' + x[-42:-40] + ':' + x[-40:-38] + ' UTC, 30 May '
+ x[-49:-45], y = 1.025, ma = 'center', size = 15, weight = 'bold') #Add AOD colorbar AOD_Colorbar() if AOD.count() > 0: #Create filled contour plot of AOD data Plot = ax.contourf(Lon, Lat, AOD, data_range, cmap = color_map, extend = 'both', zorder = 3,
transform = ccrs.PlateCarree()) else: pass #Show figure plt.show() #Save figure as a .png file #dpi sets the resolution of the digital image in dots per inch #Find the saved figures in your "current working directory" (folder containing Python Notebook file,-
#netCDF-4 data files)
filename = 'G16_CONUS_ABI_AOD_20190530_' + x[-42:-38] fig.savefig(filename, bbox_inches = 'tight', dpi = 150) #Erase plot so we can build the next one plt.close()
ABI AOD 12 Figs Output

Block #14: Make an animation of AOD figures using python image library (Pillow)

Now that you have 12 ABI AOD figures, representing 1 hour of observations, let’s animate them to see how the AOD evolved. There are several Python different animation libraries; we use Pillow here because it retains the features of continuous color bars relatively well.

Note that the “glob( )” function in this block collects all of the .png files that contain “ABI_AOD_20190530_18” in the filename from your current working directory, which should capture only the files we made in Block #13. If for some reason you have any additional.png files in your current working directory with the same string of text in the filename, then they will be added to the animation. If this is the case, either move those files or rename them.

The “for” loop cycles through the 12 graphics files, opens each one, and adds it to the “frames” list.

There are several parameters in “frames[0].save( )” you can use to customize the animation. The time delay in ms between each frame (how fast or slow the animation runs) is set by “duration”; “duration = 1000” sets 1 second between frames. The time delay before the animation restarts is set by “loop”; “loop = 0” means the animation runs continuously with no delay.

Find the saved animation in your “current working directory” in this tutorial, that is your folder on Jupyter Notebook where you saved the Python Notebook .ipynb file and netCDF-4 data files.
#Block 14: Make an animation of AOD figures using python image library (Pillow)
#Pillow is preferred for AOD animations because it retains the features of continuous colorbars relatively well

#Collect all of the AOD graphics files (figures) in given subdirectory
file_list = sorted(glob.glob(os.getcwd() + '/*ABI_AOD_20190530_18*.png'))

#Create an empty list to store figures
frames = []

#Loop through graphics files and append
for x in file_list:
    new_frame = Image.open(x)
    frames.append(new_frame)

#Save animation
#Find the saved animation in your "current working directory" (folder containing Python Notebook file,- 
#netCDF-4 data files)
#Duration is speed of frame animation in ms (e.g., 1000 ms = 1 second between frames) #Loop sets time before animation restarts (e.g., loop = 0 means animation loops continuously with no delay) frames[0].save('G16_CONUS_ABI_AOD-Animation_20190530_18-19.gif', format = 'GIF', append_images = frames[1:],
save_all = True, duration = 1000, loop = 0) #Close the graphics files we opened for x in file_list: new_frame.close() print('Animation done!')
G16 CONUS ABI AOD Animation 20190530 1819