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 Block 15 Block 16

Tutorial: Processing and Visualizing VIIRS Active Fires Satellite Data

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


The Visible Infrared Imaging Radiometer Suite (VIIRS) instrument on the Joint Polar Satellite Series (JPSS) satellites, S-NPP and NOAA-20, is a scanning radiometer used to measure properties of the Earth’s atmosphere, ocean, and land. VIIRS extends the heritage of MODIS on the Terra and Aqua satellites and of AVHRR on the NOAA and Metop series of satellites. VIIRS has 22 spectral bands ranging from the visible to longwave infrared: 16 moderate resolution band (M-bands) with 750 m spatial resolution at nadir, 5 imaging resolution bands (I-bands) with 375 m spatial resolution at nadir, and one panchromatic day-night band (DNB) with 750 m spatial resolution. VIIRS observations can be used individually or in combination (in scientific algorithms) to measure such properties as fires, clouds, aerosols, ocean color, sea surface temperature, land surface temperature, ice motion, and albedo.

S-NPP and NOAA-20 are in low altitude sun-synchronous orbits, approximately 50 minutes apart, with a 1:30 PM ascending node. This means the satellites have global coverage, with daytime observations occurring in the early afternoon local time. VIIRS has a swath width of 3060 km, providing full global coverage daily in the day and night sides of Earth.

The VIIRS Active Fires product demonstrated in this tutorial has 750 m resolution and includes a variety of fire datasets. The most basic is fire detections, which simply show the locations of fires sensed by VIIRS with high, nominal, and low confidence. A new dataset, available since December 12, 2019, is the fire persistent anomalies. These are the locations of potential fire “false alarms” due to sources such as oil/gas industry flaring and solar panels. These locations are flagged as persistent anomalies because they are suspected false alarms, so the user must interpret the satellite information for these fires with caution. In some cases, the fire persistent anomalies may indicate actual fires, such as fires next to a solar farm in California south of the Salton Sea or lava from a volcano igniting vegetation in Hawaii. Also included is fire radiative power (FRP), a measure of fire intensity.

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 VIIRS Active Fires data; how to process fire datasets for an event on March 8, 2020; and how to plot nominal and high confidence fire detections, fire persistent anomalies, and FRP on maps to create professional-looking figures.

NOAA satellite data are distributed in netCDF-4 file format. You can download the 4 tar files (containing netCDF-4 VIIRS Active Fires data files) used in this tutorial by clicking 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 VIIRS Active Fires data. We need libraries to do things like perform array operations (numpy), plot data (matplotlib), make maps (cartopy), read netCDF-4 files (netCDF4, Dataset), and open .tar files (tarfile).
#Block 1: Import libraries and settings

#To perform array operations
import numpy as np 

#Main plotting libraries
import matplotlib as mpl
from matplotlib import pyplot as plt
import matplotlib.ticker as ticker
from matplotlib.ticker import AutoMinorLocator

#Libaries for drawing figures and map projections
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 accessing files in the directory
import os

#To read in netCDF files
from netCDF4 import Dataset

#Library for working with tar files
import tarfile

#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: Extract tar files

Because the JPSS satellites have global coverage, there are over 1000 VIIRS netCDF-4 files produced each day per satellite. To make file distribution easier, the netCDF-4 files are bundled in tar files (.tar). If you download VIIRS data from NOAA’s CLASS, for example, the data will be delivered as tar files. Therefore, before we can start working with the VIIRS Active Fires data, we have to unpack these tar files.
#Block 2: Extract tar files

file_list = glob.glob(os.getcwd() + '/*.tar')

for x in file_list:
    tar = tarfile.open(x)
    tar.extractall(path = os.getcwd())
    tar.close

The “for” loop in Block #2 will progress through the 4 VIIRS Active Fires tar files in your “current working directory” and extract the netCDF-4 files (.nc ) using the tarfile module. The extracted netCDF-4 files will appear in your “current working directory” in this tutorial, that is the same folder on Jupyter Notebook where you saved the Python Notebook .ipynb file and the tar files.


Block #3: Explore a VIIRS Active Fires Data File

Let’s open a netCDF-4 data file containing VIIRS Active Fires data. Each VIIRS netCDF-4 file contains a granule of VIIRS data. VIIRS scans a swath (cross-track) as it orbits on S-NPP or NOAA-20. A granule corresponds to 48 scans, so each VIIRS netCDF-4 file contains a small subset of the global VIIRS observations. This is different from ABI data files, which contain all of the ABI observations for a given view (e.g., Full Disk, CONUS, mesoscale) at a particular time step.
#Block 3: Explore a VIIRS Active Fires data file

#Enter file name
fname = os.getcwd() + '/AF_v1r2_j01_s202003081844102_e202003081845348_c202003081911360.nc'

#Set the file name to read
file_id = Dataset(fname)

#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) #Check the "groups" metadata ##print(file_id.groups) #Check the "Fire Pixels" metadata ##print(file_id.groups['Fire Pixels']) #Check the "FP_power" metadata ##print(file_id.groups['Fire Pixels'].variables['FP_power']) #Check the "FP_power" array values ##print(file_id.groups['Fire Pixels'].variables['FP_power'][:]) #Check the "FP_PersistentAnomalyCategory" metadata #ONLY for data files beginning 12 Dec 2019 ##print(file_id.groups['Fire Pixels'].variables['FP_PersistentAnomalyCategory']) #Check the "FP_PersistentAnomalyCategory" array values #ONLY for data files beginning 12 Dec 2019 ##print(file_id.groups['Fire Pixels'].variables['FP_PersistentAnomalyCategory'][:]) #Check the "fire_mask" metadata ##print(file_id.groups['Fire Mask'].variables['fire_mask']) #Check the "fire_mask" array values ##print(file_id.groups['Fire Mask'].variables['fire_mask'][:,:])
The VIIRS data file names are very long! The graphic below shows how to decipher the data file name:
VIIRS File Names
Some Notes about netCDF-4 Files
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 >20 variables in the Fire Pixels group in an Active Fires 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 #3, 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. We need to explore the “groups” sub-group to find the data in which we are interested. In the “groups” metadata, the “fire_mask” variable in the “Fire Mask” group, and the “FP_latitude”, “FP_longitude”, FP_power”, and “FP_PersistentAnomalyCategory” variables in the “Fire Pixels” group are what we need.

The metadata for “FP_power” indicate that this is the fire radiative power (FRP) variable. Because the data are 1-dimensional (1D), we add [:] to the end of the variable file path name to display the FRP array values.

The metadata for “FP_ PersistentAnomalyCategory” indicate that this is the new fire persistent anomalies dataset indicating potential fire “false alarms”. The metadata shows us that the variable is an integer with six possible values:

Similar to the FRP data, the fire persistent anomalies data are 1D, so we add [:] to the end of the variable file path name to display the array values. In this granule of data, most of the values are “0”, indicating the pixels correspond to no anomaly, but there are two values of “5” indicating “unclassified” anomalies.

The metadata for the “fire_mask” variable indicate this variable is an integer with ten possible values, which classify each pixel in the swath as non-fire (0-6), low confidence fire (7), nominal confidence fire (8), or high confidence fire (9). Since this variable is 2D, we add [: , :] to the end of the variable path name to display a snippet of the array values.


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.

#Block 4: Check the units for the variables of interest

print('Fire Radiative Power unit is', (file_id.groups['Fire Pixels'].variables['FP_power'].__getattr__('units')))
print('Latitude unit is', (file_id.groups['Fire Pixels'].variables['FP_latitude'].__getattr__('units')))
print('Longitude unit is', (file_id.groups['Fire Pixels'].variables['FP_longitude'].__getattr__('units')))
In this case, FRP is units of MW. Latitude and longitude are in degrees. This is different from ABI data files, which contain latitude and longitude in radians, to save file space. Don’t get confused when working with different satellite datasets!


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. Although it is not the case for VIIRS Active Fires data, some satellite data files store variables as integers, which require conversion to floating point numbers using a scale and offset. It is a good habit to check the data types for satellite variables to ensure you know what you’re working with!

#Block 5: Check the data type for the variables of interest

print('Fire Radiative Power data type is', (file_id.groups['Fire Pixels'].variables['FP_power'][:].dtype))
print('Latitude data type is', (file_id.groups['Fire Pixels'].variables['FP_latitude'][:].dtype))
print('Longitude data type is', (file_id.groups['Fire Pixels'].variables['FP_longitude'][:].dtype))
print('Fire Persistent Anomalies data type is', 
(file_id.groups['Fire Pixels'].variables['FP_PersistentAnomalyCategory'][:].dtype)) print('Fire Mask data type is', (file_id.groups['Fire Mask'].variables['fire_mask'][:].dtype))
In this case, we see that the FRP, latitude, and longitude are floats, and the fire persistent anomalies and fire mask variables are integers (which we expect from looking at the metadata in Block #3), so we can proceed.


Block #6: Select and Process Fire Detections (nominal and high confidence fires)

We use the “def” command to create a function called “Fire_Detect_Data(file_id)” that captures the steps needed to read in and process nominal and high confidence fires and fire persistent anomalies. Functions allow us to complete a specific action like selecting and processing data over and over again, without having to re-write the full code each time. We will create more functions in Blocks #7, #8, #9, #12, #13, and #14.

To plot fire detections, we need to start with the “Fire Mask” variable. Recall from Block #3, this variable is equal to 7 for low confidence fires, 8 for nominal confidence fires, and 9 for high confidence fires. We want to remove all of the non-fire pixels, so after reading in the “Fire Mask” data, we mask all pixels with values < 7 (“Fire_Mask < 7”). We use “.compressed( )” to retain only the unmasked (e.g. fire) pixels. In order to plot fire pixels on a map, we will need the corresponding latitude and longitude data. Recall from Block #3 that the “Fire Mask” variable is 2D, while the latitude and longitude variables are 1D. So we use the numpy (abbreviated as “np”) “ravel” function to convert the 2D array (“Fire_2D”) to 1D (“Fire_1D”).

#Block 6: Select and process fire detections (nominal and high confidence fires) from a single file -
#(as a callable function)
def Fire_Detect_Data(file_id): #Read in "Fire Mask" variable, mask all non-fire pixels (Fire_Mask < 7), and convert from 2D to 1D array #.compressed() retains only unmasked pixels in array Fire_Mask = file_id.groups['Fire Mask'].variables['fire_mask'][:,:] Fire_Conf_Mask = (Fire_Mask < 7) Fire_Data = np.ma.masked_where(Fire_Conf_Mask, Fire_Mask) Fire_2D = Fire_Data.compressed() Fire_1D = np.ravel(Fire_2D)

Next we read in the latitude and longitude variables (“FP_latitude”, “FP_longitude”). The longitude variable designates degrees West longitude using negative values (e.g., -100 = 100 degrees West). In this tutorial, we are focusing on the eastern U.S., but since VIIRS data are global, there may be situations when you want to plot VIIRS data near the International Date Line (IDL; 180° longitude). Python has a display error when data cross the IDL and the map domain is set using negative longitude values. We circumvent this problem by converting the range of longitude degrees from [0, 180] to [0, 360]. This means all longitude values will be positive (e.g., 100 degrees West = 260).

    #Read in latitude and longitude
    #Convert longitude from [0,180] to [0,360] to avoid plotting errors near International Date Line
    Lat = file_id.groups['Fire Pixels'].variables['FP_latitude'][:]
    Lon_180 = file_id.groups['Fire Pixels'].variables['FP_longitude'][:]
    Lon = (Lon_180 % 360)

We want to plot only the nominal and high confidence fires. So next, we select the pixels corresponding to nominal and high confidence fires (“Fire_Final”) and apply them to the latitude (“Lat”) and longitude (“Lon”). This gives us 1D arrays (“Lat_Detect” and “Lon_Detect”) containing only the latitude/longitude corresponding to nominal and high confidence fire pixels; we will plot these locations in Blocks #10 and #11.

    #Select nominal confidence (fire_mask = 8) and high confidence (fire_mask = 9) fire pixels
    Fire_Final = (Fire_1D == 8) | (Fire_1D == 9)
    Lat_Detect = Lat[Fire_Final]
    Lon_Detect = Lon[Fire_Final]

In the plots we are making, we also want to show the locations of the fire potential “false alarms” due to sources of potential fire persistent anomalies. Recall from Block #3, the fire persistent anomalies variable (“FP_ PersistentAnomalyCategory”) is equal to 0 for no anomaly and has values of 1-5 for the various anomaly sources. Here, we only want to plot the anomalies, so we read in the dataset and mask all of the pixels corresponding to no anomaly (“PA_Flag == 0”). Again, we use “.compressed( )” to retain only the unmasked pixels; in this case, that’s all the pixels that correspond to fire persistent anomalies (“PA_Final”). The last step is to apply the “PA_Mask” to the original latitude (“Lat”) and longitude (“Lon”) and use “.compressed( )” to retain only the latitude/longitude corresponding to the fire persistent anomalies (“Lat_PA” and “Lon_PA”). We will plot the fire persistent anomalies along with the fire detections in Blocks #10 and #11.

 #Read in fire persistent anomalies and mask fire pixels with no anomalies (FP_PersistentAnomalyCategory = 0) 
    #.compressed() retains only unmasked pixels in array
    PA_Flag = file_id.groups['Fire Pixels'].variables['FP_PersistentAnomalyCategory'][:]
    PA_Mask = PA_Flag == 0
    PA_Pixels = np.ma.masked_where(PA_Mask, PA_Flag)
    PA_Final = PA_Pixels.compressed()

    #Mask latitude and longitude that correspond to pixels with no anomalies (FP_PersistentAnomalyCategory = 0)
    #.compressed() retains only unmasked pixels in array
    Lat_int = np.ma.masked_where(PA_Mask, Lat)
    Lon_int = np.ma.masked_where(PA_Mask, Lon)
    Lat_PA = Lat_int.compressed()
    Lon_PA = Lon_int.compressed()
    
    return Lon_Detect, Lat_Detect, PA_Final, Lon_PA, Lat_PA

Block #7: Create Colormap for Plotting Fire Persistent Anomalies

Here we create a function to make a colormap for the fire persistent anomalies data we will be plotting. This colormap (“PA_Colormap( )”) is a discrete colormap with a different color for each of the five anomaly categories. We chose colors that match the recommended colors from the VIIRS fire team, as shown on the NOAA JSTAR Mapper website. A list of all of the available colors in the matplotlib (abbreviated “mpl”) library is found here.

#Block 7: Create colormap for plotting fire persistent anomalies (as a callable function)

def PA_Colormap():
    #Discrete colormap for fire persistent anomalies
    anomaly_color_map = mpl.colors.ListedColormap(['lime', 'fuchsia', 'cyan', 'sienna', 'plum'])
    return anomaly_color_map

Block #8: Create Fire Persistent Anomalies Colorbar, Independent of Plotted Data

We create a function (“PA_Colorbar( )”) to display the fire persistent anomalies colorbar on the maps we will create. In Block #7, we created a colormap for the fire persistent anomalies; in this block, we put the colormap in a colorbar with labels and tick marks. The colorbar is shown below.

#Block 8: Create fire persistent anomalies colorbar, independent of plotted data (as a callable function)
#cbar_ax are dummy variables
#Location/attributes of colorbar set by .set_position (x0, y0, width, height) to scale automatically with plot

def PA_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.25, posn.y0 - 0.06, 0.5, 0.02])
    anomaly_color_map = mpl.colors.ListedColormap(['lime', 'fuchsia', 'cyan', 'sienna', 'plum'])
    bounds = [1, 2, 3, 4, 5, 6]
    norm_PA = mpl.colors.BoundaryNorm(bounds, anomaly_color_map.N)
    cb = mpl.colorbar.ColorbarBase(cbar_ax, cmap = anomaly_color_map, norm = norm_PA, 
orientation = 'horizontal', ticks = bounds, spacing = 'uniform') cb.set_label(label = 'Fire Persistent Anomalies', size = 'small', weight = 'bold') cb.ax.xaxis.set_major_formatter(ticker.NullFormatter()) cb.ax.xaxis.set_minor_locator(AutoMinorLocator(2)) cb.ax.set_xticklabels(['Oil/Gas Flare', 'Volcano', 'Solar Panel', 'Urban', 'Unclassified'],
minor = True, ha = 'center') cb.ax.tick_params(which = 'both', length = 0, labelsize = 'x-small') plt.sca(last_axes)
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, the colorbar position 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. PA Colorbar

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

There are many steps in Python to create a professional quality map. Here, we create a function (“VIIRS_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 as a general default; you can modify the labels by changing the values inside the “ax.set_xticks” and “ax.set_yticks” commands. Remember to specify longitude values in a range from 0 to 360; recall we changed the default range for the longitude variable back in Block #6, in order to avoid potential plotting errors near the IDL.

#Block 9: Format map with Plate Carree projection (as a callable function)

def VIIRS_Map_Settings_PC(ax):
    #Set up and label the lat/lon grid
    lon_formatter = LongitudeFormatter(zero_direction_label = True)
    lat_formatter = LatitudeFormatter()
    ax.xaxis.set_major_formatter(lon_formatter)
    ax.yaxis.set_major_formatter(lat_formatter)
    ax.set_xticks([0, 20, 40, 60, 80, 100, 120, 140, 160, 180, 200, 220, 240, 260, 280, 300, 320, 340, 360], 
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( )”. You can 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. Remember to specify longitude values in a range from 0 to 360. 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 [xmin, xmax, ymin, ymax]
    #Use 360 degrees for longitude coordinates (i.e, -100 (100 deg W) = 260)
    #NOTE: comment out (add leading ##) the line below to automatically set domain to extend to limits of data
    ax.set_extent([260, 290, 25, 50], crs = ccrs.PlateCarree())
Lat Lon Extent

Block #10: Plot Fire Detections, Including Fire Persistent Anomalies, from a Single Data File (granule)

Now, we put together all of the functions we created to make a plot of fire detections and fire persistent anomalies from one granule of VIIRS data, contained in the file we opened back in Block #3. We set up a figure and axes with the Plate Carrée projection; this is also the source projection for the VIIRS data. There are many other map projection options to explore; examples are listed here. We use the function (“VIIRS_Map_Settings_PC(ax)”) we made in Block #9 to add the map settings. We also add a title to our figure, using reverse indexing (from right to left) of the file name to automatically add the time in UTC and year to the title.

#Block 10: Plot fire detections, including fire persistent anomalies, from a single data file (granule)

#Set up figure and map projection: PlateCarree(central_longitude)
#Plate Carree: equidistant cylindrical projection w/equator as the standard parallel;-
#default central_longitude = 0
#Use central_longitude = 180 to avoid plotting errors near International Date Line fig = plt.figure(figsize=(8, 10)) ax = fig.add_subplot(1,1,1, projection = ccrs.PlateCarree(central_longitude = 180)) #Format map with Plate Carree projection VIIRS_Map_Settings_PC(ax) #Add and format title #Reverse indexing (from right to left) of file name automatically adds time and year to title plt.title('NOAA-20/VIIRS\nFire Detections\n' + fname[-44:-42] + ':' + fname[-42:-40] + ' UTC, 8 March ' +
fname[-52:-48], y = 1.025, ma = 'center', size = 15, weight = 'bold')

We read in and process the fire detections and fire persistent anomalies using the function (”Fire_Detect_Data(file_id)”) we made in Block #6. We create two scatter plots, one for each fire dataset. “Plot1” uses the simpler “plot” function to plot fire detections (fire locations only) using red filled circles (‘ro’) size 3 (“ms = 3’) with a fine black border (“mec = ‘black’, mew = 0.5”) as markers. All of the parameters and arguments for the “plot” function are listed here. We use the “plot” function for the fire detections because all of our markers are the same (i.e., red dots) - we don’t have to color the markers using a legend or colorbar to denote additional information.

“Plot2” uses the “scatter” function to plot fire persistent anomalies as small (“s = 8”) filled circles (default marker) colored using the colormap we made in Block #7 (“cmap = PA_Colormap( )”) with a thin black border (“linewidths = 0.5, edgecolors = ‘black’”); “vmin” and “vmax” set the min and max, respectively, of the plotted data. All of the parameters and arguments for the “scatter” function are listed here. We use the “scatter” function for the anomalies because we need to color the markers different colors corresponding to the different fire persistent anomaly categories. We also add the fire persistent anomalies colorbar (“PA_Colorbar( )”) with the function we made in Block #8. Because the fire persistent anomalies are relatively rare, they don’t necessarily occur in every granule. So we add an “if/else” statement around “Plot2” and the colorbar function to check if there is anything in the “PA_Final” array: if there is, we plot it with “Plot2” and add the colorbar to the map; if “PA_Final” is empty, we do nothing (“pass”).

#Read in and process fire detections (nominal and high confidence fires)
Lon_Detect, Lat_Detect, PA_Final, Lon_PA, Lat_PA = Fire_Detect_Data(file_id)

#Create scatter plot of fire locations, including fire persistent anomalies
Plot1 = ax.plot(Lon_Detect, Lat_Detect, 'ro', mec = 'black', mew = 0.5, ms = 3, zorder = 4, 
transform = ccrs.PlateCarree()) if len(PA_Final) > 0: Plot2 = ax.scatter(Lon_PA, Lat_PA, c = PA_Final, s = 8, linewidths = 0.5, edgecolors = 'black',
cmap = PA_Colormap(), vmin = 1, vmax = 5, zorder = 5, transform = ccrs.PlateCarree()) #Add fire persistent anomalies colorbar PA_Colorbar() else: pass

Notice that in both “Plot1” and “Plot2” 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 VIIRS 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 (full disk view) and you don’t specify “transform = ccrs.PlateCarree( )”, the satellite data will not plot correctly.

Finally, we “show” the map of plotted data (so we can see the output) and save the 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:

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 figure in your “current working directory” in this tutorial, that is your folder on Jupyter Notebook where you saved the Python Notebook .ipynb file and tar files.

#Show figure
plt.show()

#Save figure as a .png file
filename = 'N20_VIIRS_FireDet_' + fname[-52:-44] + '_' + fname[-44:-40]
fig.savefig(filename, bbox_inches = 'tight', dpi = 150)
N20 VIIRS FireDet 20200308 1844

Block #11: Plot Fire Detections, Including Fire Persistent Anomalies, from Multiple Data Files (granules) in One Figure

Since each VIIRS granule only covers a relatively small area, you will likely want to plot data from multiple files together in the same figure. We recommend plotting one granule of data first, like in Block #10, to check the data processing and map settings before moving on to multiple data files.

#Block 11: Plot fire detections, including fire persistent anomalies, from multiple data files (granules)-
#in one figure
#Collect all of the .nc files in the specified subdirectory file_list = glob.glob(os.getcwd() + '/AF*.nc') #Set up figure and map projection: PlateCarree(central_longitude) #Plate Carree: equidistant cylindrical projection w/equator as the standard parallel;-
#default central_longitude = 0
#Use central_longitude = 180 to avoid plotting errors near International Date Line fig = plt.figure(figsize=(8, 10)) ax = fig.add_subplot(1,1,1, projection = ccrs.PlateCarree(central_longitude = 180)) #Format map with Plate Carree projection VIIRS_Map_Settings_PC(ax) #Add and format title plt.title('NOAA-20/VIIRS\nFire Detections\n8 Mar 2020', y = 1.025, ma = 'center', size = 15, weight = 'bold')

Plotting data from multiple files is very similar to plotting data from one file, so the main content of Block #11 is the same as in Block #10. We add the “glob( )” function to retrieve all of the VIIRS netCDF-4 (.nc) files in our “current working directory” and a “for” loop to cycle through those files and process and plot the fire detections and fire persistent anomalies from each file.

#Set up a loop to plot the data from multiple files
for x in file_list:
    file_id = Dataset(x)
    
    #Check if data file contains fire data; if file is empty, skip it
    if len(file_id.groups['Fire Pixels'].variables['FP_power'][:]) == 0:
        continue
    
    #Read in and process fire detections and fire persistent anomalies data
    Lon_Detect, Lat_Detect, PA_Final, Lon_PA, Lat_PA = Fire_Detect_Data(file_id)
    
    #Create scatter plot of fire locations, including fire persistent anomalies
    if len(Lat_Detect) > 0:
        Plot1 = ax.plot(Lon_Detect, Lat_Detect, 'ro', mec = 'black', mew = 0.5, ms = 3, zorder = 4, 
transform = ccrs.PlateCarree()) else: continue if len(PA_Final) > 0: Plot2 = ax.scatter(Lon_PA, Lat_PA, c = PA_Final, s = 8, linewidths = 0.5, edgecolors = 'black',
cmap = PA_Colormap(), vmin = 1, vmax = 5, zorder = 5, transform = ccrs.PlateCarree()) #Add fire persistent anomalies colorbar PA_Colorbar() else: continue

Inside the loop, we add an additional “if” statement to check each file to see if it contains any fire data, using “FP_power” as a representative variable. If there are fire data in the file, the length of the “FP_power” array will be > 0, and we proceed with processing and plotting the fire detections and fire persistent anomalies in that file. If there are no fire data in the file, the “FP_power” array will be empty, and we “continue” to the next data file.

#Show figure
plt.show()

#Save figure as a .png file
fig.savefig('N20_VIIRS_FireDet_20200308', bbox_inches = 'tight', dpi = 150)   
N20 VIIRS FireDet 20200308

Block #12: Select and Process Fire Radiative Power (FRP) Data

Fire detections show where fires are located, but convey no information about the fires themselves. In some cases, you will want to display the intensity of fires, represented by the FRP. This block is designed to work independently of Block #6, for cases when you want to plot fire detections or FRP, but not both. We start by re-entering the file name of the same data file from Block #3. Then we create a function (“FRP_Data(file_id)”) that captures the steps needed to read in and process the FRP data and remove any fire persistent anomalies. We read in the FRP data (“FP_power”) and latitude and longitude (“FP_latitude”, “FP_longitude”) and convert the range of longitude degrees from [0, 180] to [0, 360], just like in Block #6.

#Block 12: Select and process fire radiative power (FRP) data from a single file (as a callable function)

#Enter file name
fname = os.getcwd() + '/AF_v1r2_j01_s202003081844102_e202003081845348_c202003081911360.nc'

#Set the file name to read
file_id = Dataset(fname)

def FRP_Data(file_id):
    #Read in FRP data and convert to log10
    FRP = file_id.groups['Fire Pixels'].variables['FP_power'][:]

    #Read in latitude and longitude
    #Convert longitude from [0,180] to [0,360] to avoid plotting errors near International Date Line
    Lat = file_id.groups['Fire Pixels'].variables['FP_latitude'][:]
    Lon_180 = file_id.groups['Fire Pixels'].variables['FP_longitude'][:]
    Lon = (Lon_180 % 360)

For the plot of FRP data, we want to use the fire persistent anomalies (“FP_ PersistentAnomalyCategory”) dataset a bit differently from what we did in Block #6. Instead of plotting the anomalies on a map with the FRP values, instead we want to use the anomalies variable to remove any FRP values that correspond to anomalies; since the anomalies indicate potential fire “false alarms”, - we don’t want to plot FRPs for the fire anomalies. To do this, we read in the anomalies variable (just like in Block #6) but here, we make a mask (“FRP_Mask”) for all of the pixels corresponding to fire persistent anomalies (“PA_Flag != 0”) and then apply that mask to the FRP dataset. Then we use “.compressed( )” to retain only the unmasked pixels; in this case, that’s all the FRP pixels the correspond to pixels with “no anomaly” (“FRP_Final”). Observed FRP values cover a wide range, varying from 0 to 5000 MW in VIIRS Active Fires data files. For any variable that has a large dynamic range like this, it is customary to plot the data using a log scale. As a result, we will map FRP as log10(FRP) in order to display the range in FRP observations most realistically. So as a final step, we take the log10 (“FRP_log10_Final”).

#Mask fire pixels corresponding to fire persistent anomalies (FP_PersistentAnomalyCategory != 0)
    #.compressed() retains only unmasked pixels in array
    PA_Flag = file_id.groups['Fire Pixels'].variables['FP_PersistentAnomalyCategory'][:]
    FRP_Mask = PA_Flag != 0
    FRP_Pixels = np.ma.masked_where(FRP_Mask, FRP)
    FRP_Final = FRP_Pixels.compressed()
    FRP_log10_Final = np.log10(FRP_Final)

Then we need to apply the “FRP_Mask” to the original latitude (“Lat”) and longitude (“Lon”) and use “.compressed( )” to retain only the latitude/longitude corresponding to the FRP pixels with the fire persistent anomalies removed (“Lat_FRP” and “Lon_FRP”). In Blocks #15 and #16, we will plot the FRP, corrected for the fire persistent anomalies.

    #Mask latitude and longitude that correspond to fire persistent anomalies (FP_PersistentAnomalyCategory != 0)
    #.compressed() retains only unmasked pixels in array
    Lat_int2 = np.ma.masked_where(FRP_Mask, Lat)
    Lon_int2 = np.ma.masked_where(FRP_Mask, Lon)
    Lat_FRP = Lat_int2.compressed()
    Lon_FRP = Lon_int2.compressed()
    
    return FRP_log10_Final, Lon_FRP, Lat_FRP

Block #13: Create Colormap for Plotting FRP Data

We create a function (“FRP_Colormap( )”) to make a custom continuous colormap for the FRP data we will be plotting, using 11 colors from the matplotlib library. Although FRP values range from 0 to 5000 MW in VIIRS Active Fires data files, as mentioned previously, most fires have FRP values < 1000 MW. To capture the range in FRP observations most realistically, we will set a data maximum corresponding to 1000 MW on our plots, but we still want to represent fires that have FRP > 1000 MW. So we pick a different color (‘darkred’) to plot FRP data that are > 1000 MW.
#Block 13: Create colormap for plotting FRP data (as a callable function)

def FRP_Colormap():
    #Custom continuous colormap for FRP
    #.set_over sets color for plotting data > max
    frp_color_map = mpl.colors.LinearSegmentedColormap.from_list('custom_FRP', [(0, 'yellow'),(0.2, 'gold'), 
(0.4, 'orange'), (0.6, 'darkorange'), (0.8, 'red'), (1, 'firebrick')], N = 150) frp_color_map.set_over('darkred') return frp_color_map


Block #14: Create FRP Colorbar, Independent of Plotted Data

We create a function (“FRP_Colorbar( )”) to display the FRP colorbar on the maps we will create of FRP data. In Block #13, we created a colormap for the FRP data; in this block, we put the colormap in a colorbar with labels and tick marks. The colorbar is shown below.

#Block 14: Create FRP colorbar, independent of plotted data (as a callable function)
#cbar_ax are dummy variables
#Location/attributes of colorbar set by .set_position (x0, y0, width, height) to scale automatically with plot

def FRP_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.06, 0.3, 0.02])
    frp_color_map = mpl.colors.LinearSegmentedColormap.from_list('custom_FRP', [(0, 'yellow'),(0.2, 'gold'), 
(0.4, 'orange'), (0.6, 'darkorange'), (0.8, 'red'), (1, 'firebrick')], N = 150) frp_color_map.set_over('darkred') norm_FA = mpl.colors.Normalize(vmin = 0, vmax = 3) cb = mpl.colorbar.ColorbarBase(cbar_ax, cmap = frp_color_map, norm = norm_FA, orientation = 'horizontal',
ticks = [0, 1, 2, 3], extend = 'max') cb.set_label(label = 'FRP (log$_{10}$, MW)', size = 'small', weight = 'bold') cb.ax.set_xticklabels(['0', '1.0', '2.0', '3.0']) cb.ax.tick_params(labelsize = 'small') plt.sca(last_axes)
The code in this block is similar to that in Block #8; this colorbar is also independent of the plotted data and scales automatically. FRP Colorbar


Block #15: Plot FRP from a Single Data File (granule)

This block is similar to Block #10 in that we put together all of the functions we created to plot the FRP data (corrected for fire persistent anomalies) from one granule of VIIRS data contained in the file we opened in Block #12. We set up the figure and axes with the Plate Carrée projection; add the map settings, title, and FRP colorbar; and read in and process the FRP data.

#Block 15: Plot FRP from a single data file (granule)

#Set up figure and map projection: PlateCarree(central_longitude)
#Plate Carree: equidistant cylindrical projection w/equator as the standard parallel;- 
#default central_longitude = 0
#Use central_longitude = 180 to avoid plotting errors near International Date Line fig = plt.figure(figsize=(8, 10)) ax = fig.add_subplot(1,1,1, projection = ccrs.PlateCarree(central_longitude = 180)) #Format map with Plate Carree projection VIIRS_Map_Settings_PC(ax) #Add and format title #Reverse indexing (from right to left) of file name automatically adds time and year to title plt.title('NOAA-20/VIIRS\nFire Radiative Power\n' + fname[-44:-42] + ':' + fname[-42:-40] + ' UTC, 8 March ' +
fname[-52:-48], y = 1.025, ma = 'center', size = 15, weight = 'bold') #Add FRP colorbar FRP_Colorbar() #Read in and process FRP data FRP_log10_Final, Lon_FRP, Lat_FRP = FRP_Data(file_id)

The main difference for this plot is that we have only one set of data, so we only need one “scatter” function to plot FRP as small (“s = 8”) filled circles (default marker) colored using the colormap we made in Block #13 (“cmap = FRP_Colormap( )”) with a thin black border (“linewidths = 0.5, edgecolors = ‘black’”); “vmin” and “vmax” are 0 and 3, respectively, corresponding to a range of 0-1000 MW (in log10 MW).

#Create scatter plot of FRP data, with persistent fire anomalies removed
Plot = ax.scatter(Lon_FRP, Lat_FRP, c = FRP_log10_Final, s = 8, linewidths = 0.5, edgecolors = 'black', 
cmap = FRP_Colormap(), vmin = 0, vmax = 3, zorder = 4, transform = ccrs.PlateCarree()) #Show figure plt.show() #Save figure as a .png file filename = 'N20_VIIRS_FRP_' + fname[-52:-44] + '_' + fname[-44:-40] fig.savefig(filename, bbox_inches = 'tight', dpi = 150)
N20 VIIRS FRP 20200308 1844
Find the saved figure in your “current working directory” in this tutorial, that is your folder on Jupyter Notebook where you saved the Python Notebook .ipynb file and tar files.


Block #16: Plot FRP from Multiple Data Files (granules) in One Figure

As with the fire detections data, plotting FRP from multiple data files is very similar to plotting data from one file, so the main content of Block #16 is the same as in Block #15. We add the “glob( )” function to retrieve all of the VIIRS netCDF-4 (.nc) files in our “current working directory” and a “for” loop to cycle through those files and process and plot the FRP data (corrected for fire persistent anomalies) from each file.

#Block 16: Plot FRP from multiple data files (granules) in one figure

#Collect all of the .nc files in the specified subdirectory
file_list = glob.glob(os.getcwd() + '/AF*.nc')

#Set up figure and map projection: PlateCarree(central_longitude)
#Plate Carree: equidistant cylindrical projection w/equator as the standard parallel;- 
#default central_longitude = 0
#Use central_longitude = 180 to avoid plotting errors near International Date Line fig = plt.figure(figsize=(8, 10)) ax = fig.add_subplot(1,1,1, projection = ccrs.PlateCarree(central_longitude = 180)) #Format map with Plate Carree projection VIIRS_Map_Settings_PC(ax) #Add and format title plt.title('NOAA-20/VIIRS\nFire Radiative Power\n8 Mar 2020', y = 1.025, ma = 'center', size = 15,
weight = 'bold') #Add FRP colorbar FRP_Colorbar() #Set up a loop to plot the data from multiple files for x in file_list: file_id = Dataset(x) #Check if data file contains FRP data; if file is empty, skip it if len(file_id.groups['Fire Pixels'].variables['FP_power'][:]) == 0: continue #Read in and process FRP data FRP_log10_Final, Lon_FRP, Lat_FRP = FRP_Data(file_id) #Create scatter plot of FRP data, with persistent fire anomalies removed Plot = ax.scatter(Lon_FRP, Lat_FRP, c = FRP_log10_Final, s = 8, linewidths = 0.5, edgecolors = 'black',
cmap = FRP_Colormap(), vmin = 0, vmax = 3, zorder = 4, transform = ccrs.PlateCarree()) #Show figure plt.show() #Save figure as a .png file fig.savefig('N20_VIIRS_FRP_20200308', bbox_inches = 'tight', dpi = 150)
N20 VIIRS FRP 20200308 Find the saved figure in your “current working directory” in this tutorial, that is your folder on Jupyter Notebook where you saved the Python Notebook .ipynb file and tar files.