Since reflectance data is generated from a passive energy source (the sun), data collected in cloudy sky conditions are not directly comparable to data collected in clear-sky conditions, as overhead clouds can obscure the incoming light source. AOP aims to collect data only in optimal (<10% cloud-cover) weather conditions, but cannot always do so due to logistical constraints. The flight operators record the weather conditions during each flight, and this information is passed through to the final data product at the level of the flight line (as cloud conditions can change throughout the day). Cloud conditions are reported as green (<10% cloud cover), yellow (10-50% cloud cover), or red (>50% cloud cover). The figure below shows some examples of what the cloud conditions look like at different flights collected in the three different weather classes (green, yellow, and red).

Note that there is an important distinction between airborne and satellite reflectance data. Satellite data is collected in all weather conditions, and the clouds are below the sensor, so algorithms can be generated to filter out cloudy pixels. With aerial data, we have more control over when the data are collected, to a degree. However, clouds may be present overhead, if it were deemed necessary to collect in sub-optimal weather conditions. AOP typically will only collect in "red" sky conditions if we are running out of time in a Domain and the weather isn't forecasted to improve. Since the clouds won't appear in the actual data, maintaining this record of cloud conditions is essential for properly understanding the data, and using it for change detection or other research applications. For a more direct comparison of reflectance values, we recommend only working with the clear-weather data. This lesson outlines how to do this in GEE.

Read in the AOP SDR 2019 Dataset at SOAP

We will start at our ending point of the last tutorial. For this exercise, we will read in data from the NEON site Soaproot Saddle (SOAP) collected in 2019:

// Filter image collection by date and site to pull out a single image
var soapSDR = ee.ImageCollection("projects/neon-prod-earthengine/assets/HSI_REFL/001")
  .filterDate('2019-01-01', '2019-12-31')
  .filterMetadata('NEON_SITE', 'equals', 'SOAP')
  .first();

Display the QA Bands

From the previous lesson, recall that the SDR images include 442 bands. Bands 0-425 are the data bands, which store the spectral reflectance values for each wavelength recorded by the NEON Imaging Spectrometer (NIS). The remaining bands (426-441) contain metadata and QA information that are important for understanding and properly interpreting the hyperspectral data. The data bands all follow the naming convetion B001, B002, ..., B426, and the QA bands start with something other than the letter "B", so we can use that information to extract the QA bands.

// Pull out and display only the qa bands (these all start with something other than B)
// '[^B].*' is a regular expression to pull out bands that don't start with B
var soapSDR_qa = soapSDR.select('[^B].*') 
print('QA Bands',soapSDR_qa)

Most of these QA bands are inputs to and outputs from the Atmospheric Correction (ATCOR). We will expand upon these further, and encourage you to read more details about these data in the NEON Imaging Spectrometer Radiance to Reflectance Algorithm Theoritical Basis Document. For the purposes of this exercise, we will focus on the Weather Quality Indicator band. Note that you can explore each of the QA bands, following similar steps below, adjusting the band names and values accordingly.

Read in the Weather_Quality_Indicator Band

The weather information, called Weather_Quality_Indicator is one of the most important pieces of QA information that is collected about the NIS data, as it has a direct impact on the reflectance values.

These next lines of code pull out the Weather_Quality_Indicator band, select the "green" weather data from that band, and apply a mask to keep only the clear-weather data, which is saved to the variable soapSDR_clear.

// Extract a single band Weather Quality QA layer
var soapWeather = soapSDR.select(['Weather_Quality_Indicator']);

// Select only the clear weather data (<10% cloud cover)
var soapClearWeather = soapWeather.eq(1); // 1 = 0-10% cloud cover

// Mask out all cloudy pixels from the SDR image
var soapSDR_clear = soapSDR.updateMask(soapClearWeather);

Plot the weather quality band data

For reference, we can plot the weather band data, using AOP's stop-light (red/yellow/green) color scheme, with the code below:

// center the map at the lat / lon of the site, set zoom to 12
Map.setCenter(-119.25, 37.06, 11);

// Define a palette for the weather - to match NEON AOP's weather color conventions
var gyrPalette = [
  '00ff00', // green (<10% cloud cover)
  'ffff00', // yellow (10-50% cloud cover)
  'ff0000' // red (>50% cloud cover)
];

// Display the weather band (cloud conditions) with the green-yellow-red palette
Map.addLayer(soapWeather,
             {min: 1, max: 3, palette: gyrPalette, opacity: 0.3},
             'SOAP 2019 Cloud Cover Map');

Plot the clear-weather reflectance data

Finally, we can plot a true-color image of only the clear-weather data, from soapSDR_clear that we created earlier:

// Create a 3-band cloud-free image 
var soapSDR_RGB = soapSDR_clear.select(['B053', 'B035', 'B019']);

// Display the SDR image
Map.addLayer(soapSDR_RGB, {min:103, max:1160}, 'SOAP 2019 Reflectance RGB');

Plot acquisition dates

We can apply the same concepts to explore another one of the QA bands, this time let's look at the Acquisition_Date. This may be useful if you are trying to find the dates that correspond to field data you've collected, or you want to scale up to satellite data, for example. To determine the minimum and maximum dates, you can use reduceRegion with the reducer ee.Reducer.minMax() as follows. Then use these values in the visualization parameters.

// Get the minimum and maximum values of the soapDates band
var minMaxValues = soapDates.reduceRegion({reducer: ee.Reducer.minMax(),maxPixels: 1e10})
print('min and max dates', minMaxValues);
	
// Extract acquisition dates and map, don't display layer by default
var soapDates = soapSDR.select(['Acquisition_Date']);
Map.addLayer(soapDates,
            {min:20190612, max:20190616, opacity: 0.5},
            'SOAP 2019 Acquisition Dates',0);

Recap

In this lesson you learned how to read in Weather Quality Information from the Surface Directional Reflectance (SDR) QA bands in GEE. You learned to mask data to keep only data collected in the cleares sky conditions (<10% cloud cover), and plot the three weather quality classes. You also learned how to find the other QA bands, and following a similar approach could explore each of these bands similarly. Filtering by the weather quality is an important first pre-processing step to working with NEON hyperspectral data, and is essential for interpreting the data and carrying out subsequent data analysis.

Get Lesson Code

AOP GEE SDR Weather Quality

Writing a Function to Visualize AOP Reflectance Image Collections

In the earlier Reflectance pre-processing tutorial, we showed how to read in a band of data containing weather quality information and apply cloud-masking using that band. In this tutorial, we will show you a more simplified way of doing this, using functions so we can more easily apply that operation to an Image Collection. This is called "refactoring". In any coding language, if you notice you are writing very similar lines of code repeatedly, it may be an opportunity to create a function. For example, in the previous tutorial, we repeated lines of code to pull in different years of data at SRER, the only difference being the year and the variable names for each year. As you become more proficient with GEE coding, it is good practice to start writing functions to make your scripts more readable and reproducible.

GEE Function and Mapping Syntax

Let's get started! First let's take a look at the syntax for writing user-defined functions in GEE. If you are familiar with other programming languages, this should look somewhat familiar. The function requires input argument(s) args and returns an output.

var functionName = function(args) {
  // do something with input args
  return output;
};

To call the function for a full image collection, you can use a map to apply the function to all items in a collection. This is outlined in the code chunk below.

// Map the function over the collection.
var newVariable = collection.map(myFunction);

Read in AOP Directional Reflectance Image Collection

First, we'll read in the AOP Directional Reflectance Image Collection at Jones Center At Ichauway (JERC).

// specify center location of the site (JERC)
var site_center = ee.Geometry.Point([-84.468623,31.194839])

// read in the AOP Directional Reflectance (HSI_REFL/001) Image Collection
// filter to the site_center
var sdr_col = ee.ImageCollection('projects/neon-prod-earthengine/assets/HSI_REFL/001')
  .filterBounds(site_center)

print('NEON AOP Directional Reflectance Image Collection',sdr_col)

Cloud Masking Function

Building off the example from the previous tutorial, we can write a simple function to apply cloud-masking to an Image Collection:

// Function to mask out poor-weather data, keeping only the <10% cloud cover weather data
function clearSDR(image) { 
  // create a single band Weather Quality QA layer 
  var weather_qa = image.select(['Weather_Quality_Indicator']);
  // WEATHER QUALITY INDICATOR = 1 is < 10% CLOUD COVER
  var clear_qa = weather_qa.eq(1);
  // mask out all cloudy pixels from the reflectance image
  return image.updateMask(clear_qa);
}

Then we can use .map() to apply this as follows:

// Use map to apply this function on all the NISImages and return a clear-weather collection
var sdr_cloudfree = sdr_col.map(clearSDR)

Function to Add Multiple Reflectance Layers to Map

For the next example, we will write a function to add a Map Layer for each Image in an Image collection. We'll provide the full script below, including the function addNISImage, with comments explaining what each part of the function does. Note that a helpful troubleshooting technique is to add in print statements if you are unsure what the code is returning. We have included some commented-out print statements in the function, which show the outputs (which would show up in the console tab).

For a little more detail on how this function was applied, refer to this GIS Stack Exchange Post: Add/display all images of my collection in google earth engine. When writing your own GEE code, the Google earth-engine developers pages may not always have an example of what you are trying to do, so stack overflow can be a valuable resource.

// Define visualization parameters for the reflectance data, showing a true-color image
// B053 ~ 642 nm, B035 ~ 552 nm , B019 ~ 472nm 
var sdr_vis_params = {'min':0, 'max':1200, 'gamma': 0.9, 'bands': ['B053','B035','B019']};

// Function to display each NIS Image in the NEON AOP Image Collection
function addNISImage(image) { 
  // get the system:id and convert to string
  var imageId = ee.Image(image.id);
  // get the system:id - this is an object on the server
  var sysID_serverObj = ee.String(imageId.get("system:id"));
  // getInfo() converts to string on the server
  var sysID_serverStr = sysID_serverObj.getInfo()
  // truncate the string to show only the fileName (NEON domain + site code + product code + year)
  var fileName = sysID_serverStr.slice(51,100); 
  // print("fileName: "+fileName) // optionally print the file name, can uncomment

  // add this layer to the map using the true-color (RGB) visualization parameters
  Map.addLayer(imageId, sdr_vis_params, fileName + ' Refl RGB - All Flightlines')
}

Note that the first half of this function is just pulling out relevant information about the site in order to properly label the layer on the Map display.

Function including Cloud-Masking and Weather QA Layers

Next we can build upon this function to include some small pre-processing steps, such as selecting the Weather_Quality_Indicator band, plotting it, and masking the reflectance data to include only the clear-weather (<10% cloud cover) data and add that masked dataset to the Map.

// Next we can build upon this function to add the Weather QA layer and
// call the cloud-masking function and add the clear-weather data layers to the Map

// Define a palette for the weather - to match NEON AOP's weather color conventions
// This will be used in the visualization parameters for the Weather QA layer
// green (<10% cloud cover), yellow (10-50% cloud cover), red (>50% cloud cover)
var gyr_palette = ['green','yellow','red']

// Build upon the function to add a layer of the weather band, and 
// mask out poor-weather data (>10% cloud cover), keeping only the clear weather (<10% cloud cover)
function addClearNISImages(image) { 
  // get the system:id and convert to string
  var imageId = ee.Image(image.id);
  // get the system:id - this is an object on the server
  var sysID_serverObj = ee.String(imageId.get("system:id"));
  // getInfo() converts to string on the server
  var sysID_serverStr = sysID_serverObj.getInfo()
  // truncate the string to show only the fileName (NEON domain + site code + product code + year)
  var fileName = sysID_serverStr.slice(51,100); // optionally print, can uncomment
  
  // create a single band Weather Quality QA layer for 2016
  var weather_qa = imageId.select(['Weather_Quality_Indicator']);

  // apply the function from the beginning of this script to generate a cloud-free image
  // you can apply a function directly on a single image in this way
  var sdr_cloudfree = clearSDR(imageId)

  // add the weather QA bands to the map with the green-yellow-red palette
  Map.addLayer(weather_qa, {min: 1, max: 3, palette: gyr_palette, opacity: 0.3}, fileName + ' Weather QA Band')

  // add the clear weather reflectance layer to the map - the 0 means the layer won't be turned on by default
  Map.addLayer(sdr_cloudfree, sdr_vis_params, fileName + ' Refl RGB - Clear Skies', 0)
}

//call the clearNISImages function
sdr_col.evaluate(function(sdr_col) {
  sdr_col.features.map(addClearNISImages);
})

// call the addNISimages function
// see this link for an explanation of what is occurring:
// https://gis.stackexchange.com/questions/284610/add-display-all-images-of-mycollection-in-google-earth-engine
sdr_col.evaluate(function(sdr_col) {
  sdr_col.features.map(addNISImage);
})

// Center the map on site and set zoom level (11)
Map.centerObject(site_center, 11);

In the "Layers" tab, select and de-select the Map Layers to look at the data and the weather QA information for all of the years. This figure shows the weather quality information at JERC in 2019, where the data were collected in mixed cloud conditions. Flight operators prioritize flying the area over NEON plots in the best weather conditions when possible. As explained in the previous lesson, when working with AOP reflectance data, the weather conditions during the flights are one of the most important quality considerations.

Recap

This lesson showed examples of two functions that demonstrate how to use functional programming to add layers to a Map. As explained in introduction-to-functional-programming, "Earth Engine uses a parallel processing system to carry out computation across a large number of machines...The use of for-loops is discouraged in Earth Engine. The same results can be achieved using a map() operation where you specify a function that can be independently applied to each element. This allows the system to distribute the processing to different machines." The examples presented in this lesson show some examples of how to make use of that map operation in lieu of a for loop.

Get Lesson Code

Functions to display AOP Reflectance Image Collections in GEE

Read in the AOP Directional Reflectance Image

As should be familiar by now from the previous tutorials in this series, we'll start by pulling in the AOP data. For this exercise we will only read directional reflectance data from SOAP collected in 2021:

// Filter image collection by date and site
var soapSDR = ee.ImageCollection("projects/neon-prod-earthengine/assets/HSI_REFL/001")
  .filterDate('2021-01-01', '2021-12-31')
  .filterMetadata('NEON_SITE', 'equals', 'SOAP')
  .first();

// Create a 3-band true-color image 
var soapSDR_RGB = soapSDR.select(['B053', 'B035', 'B019']);

// Display the SDR image
Map.addLayer(soapSDR_RGB, {min:103, max:1160}, 'SOAP 2021 Reflectance RGB');

// Center the map at the lat / lon of the site, set zoom to 12
Map.setCenter(-119.25, 37.06, 12);

Extract data bands

Next we will extract only the "data" bands in order to plot the spectral information. The reflectance data contains 426 data bands, and a number of QA/Metdata bands that provide additional information that can be useful in interpreting and analyzing the data (such as the Weather Quality Information). For plotting the spectra, we only need the data bands.

// Pull out only the data bands (these all start with B, eg. B001)
var soapSDR_data = soapSDR.select('B.*')
print('SOAP SDR Data',soapSDR_data)

// Read in the properties as a dictionary
var properties = soapSDR.toDictionary()

Extract wavelength information from the properties

Similar to the code above, we can use a regular expression to pull out the wavelength information from the properties. The wavelength and Full Width Half Max (FWHM) information is stored in the properties starting with WL_FWHM_B. These are stored as strings, so the nex step is to write a funciton that converts the string to a float, and only pulls out the center wavelength value (by splitting on the "," and pulling out only the first value). This is all we need for now, but if you needed the FWHM information, you could write a similar function. Lastly, we'll apply the function using GEE .map to pull out the wavelength information. We an then print some information about what we've extracted

// Select the WL_FWHM_B*** band properties (using regex)
var wl_fwhm_dict = properties.select(['WL_FWHM_B+\\d{3}']);

// Pull out the wavelength, fwhm values to a list
var wl_fwhm_list = wl_fwhm_dict.values()
print('Wavelength FWHM list:',wl_fwhm_list)

// Function to pull out the wavelength values only and convert the string to float
var get_wavelengths = function(x) {
  var str_split = ee.String(x).split(',')
  var first_elem = ee.Number.parse((str_split.get(0)))
  return first_elem
}

// apply the function to the wavelength full-width-half-max list
var wavelengths = wl_fwhm_list.map(get_wavelengths)

print('Wavelengths:',wavelengths)
print('# of data bands:',wavelengths.length())

Interactively plot the spectral signature of a pixel

Lastly, we'll create a plot in the Map panel, and use the Map.onClick function to create a spectral signature of a given pixel that you click on. Most of the code below specifies formatting, figure labels, etc.

// Create a panel to hold the spectral signature plot
var panel = ui.Panel();
panel.style().set({width: '600px',height: '300px',position: 'top-left'});
Map.add(panel);
Map.style().set('cursor', 'crosshair');

// Create a function to draw a chart when a user clicks on the map.
Map.onClick(function(coords) {
  panel.clear();
  var point = ee.Geometry.Point(coords.lon, coords.lat);
    wavelengths.evaluate(function(wvlnghts) {
      var chart = ui.Chart.image.regions({
        image: soapSDR_data, 
        regions: point, 
        scale: 1,
        seriesProperty: 'λ (nm)', 
        xLabels: wavelengths.getInfo()
    });
    chart.setOptions({
      title: 'Reflectance',
      hAxis: {title: 'Wavelength (nm)', 
      vAxis: {title: 'Reflectance'},
      gridlines: { count: 5 }}
    });
    // Create and update the location label 
    var location = 'Longitude: ' + coords.lon.toFixed(2) + ' ' +
                   'Latitude: ' + coords.lat.toFixed(2);
    panel.widgets().set(1, ui.Label(location));
    panel.add(chart);
  })
});

When you run this code, (linked at the bottom), you will see the SOAP 2021 directional reflectance layer show up in the Map panel, along with a white figure panel. When you click anywhere in the image, the empty figure will be populated with the spectral signature of the pixel you clicked on.

Recap

In this lesson you learned how to read in wavelength information from the Surface Directional Reflectance properties in GEE, created functions to convert from one data format to another, and created an interactive plot to visualize the spectral signature of a selected pixel. You can quickly see how GEE is a powerful tool for interactive data visualization and exploratory analysis.

Get Lesson Code

AOP GEE Reflectance Plot Spectra

GEE is a great place to conduct exploratory analysis to better understand the datasets you are working with. In this lesson, we will show how to pull in AOP Surface Directional Reflectance (SDR) data, as well as the Ecosystem Structure (Canopy Height Model - CHM) data to look at interannual differences at the NEON site Great Smokey Mountains (GRSM), where the Chimney Tops 2 Fire broke out in late November 2016. NEON data over the GRSM site collected in June 2016 and October 2017 captures most of the burned area and presents a unique opportunity to study wildfire effects on the ecosystem and analysis of post-wildfire vegetation recovery. In this lesson, we will calculate the differenced Normalized Burn Ratio (dNBR) between 2017 and 2016, and also create a CHM difference raster to highlight vegetation structure differences in the burned area. We will also pull in Landsat satellite data and create a time-series of the NBR within the burn perimeter to look at annual differences.

Using remote sensing data to better understand wildfire impacts is an active area of research. In April 2023, Park and Sim published an Open Access paper titled "Characterizing spatial burn severity patterns of 2016 Chimney Tops 2 fire using multi-temporal Landsat and NEON LiDAR data". We encourage you to read this paper for an example of wildfire research using AOP remote sensing and satellite data. This lesson provides an introduction to conducting this sort of analysis in Google Earth Engine.

Functions to Read in SDR and CHM Image Collections

Let's get started. The code in the beginning of this lesson should look familiar from the previous tutorials in this series. In this first chunk of code, we will the center location of GRSM, and read in the fire perimeter as a FeatureCollection.

// Specify center location and flight box for GRSM (https://www.neonscience.org/field-sites/grsm)
var site_center = ee.Geometry.Point([-83.5, 35.7])

// Read in the Chimney Tops fire perimeter shapefile
var ct_fire_boundary = ee.FeatureCollection('projects/neon-sandbox-dataflow-ee/assets/chimney_tops_fire')

Next, we'll read in the SDR image collection, and then write a function to mask out the cloudy weather data, and use the map feature to apply this to our SDR collection at GRSM.

// Read in the SDR Image Collection at GRSM
var grsm_sdr_col = ee.ImageCollection('projects/neon-prod-earthengine/assets/HSI_REFL/001')
  .filterBounds(site_center)
  
// Function to mask out poor-weather data, keeping only the <10% cloud cover weather data
function clearSDR(image) { 
  // create a single band Weather Quality QA layer 
  var weather_qa = image.select(['Weather_Quality_Indicator']);
  // WEATHER QUALITY INDICATOR = 1 is < 10% CLOUD COVER
  var clear_qa = weather_qa.eq(1);
  // mask out all cloudy pixels from the SDR image
  return image.updateMask(clear_qa);
}

// Use map to apply the clearSDR function to the SDR collection and return a clear-weather subset of the data
var grsm_sdr_cloudfree = grsm_sdr_col.map(clearSDR)

Next let's write a function to display the NIS images from 2016, 2017, and 2021 in GEE. For more details on how this function works, you can refer to the tutorial Functions in Google Earth Engine (GEE).

// Function to display individual (yearly) SDR Images
function addSDRImage(image) { 
  var image_id = ee.Image(image.id); // get the system:id and convert to string
  var sys_id = ee.String(image_id.get("system:id")).getInfo();  // get the system:id - this is an object on the server
  var filename = sys_id.slice(52,100); // extract the fileName (NEON domain + site code + product code + year)
  var image_rgb = image_id.select(['B053', 'B035', 'B019']);  // select only RGB bands for display
  
  Map.addLayer(image_rgb, {min:220, max:1600}, filename, 1)   // add RGB composite to the map
}

// call the addNISimages function to add SDR layers to map
grsm_sdr_col.evaluate(function(grsm_sdr_col) {
  grsm_sdr_col.features.map(addSDRImage);
})

Next we can create a similar function for reading in the CHM dataset over all the years. The main differences between this function and the previous one are that 1) it is set to display a single band image, and 2) instead of hard-coding in the minimum and maximum values to display, we dynamically determine them from the data itself, so it will scale appropriately.

// Read in the CHM Image collection at GRSM
var grsm_chm_col =  ee.ImageCollection('projects/neon-prod-earthengine/assets/CHM/001')
  .filterBounds(site_center)

// Function to display Single Band Images setting display range to linear 2%
function addSingleBandImage(image) { // display each image in collection
  var image_id = ee.Image(image.id); // get the system:id and convert to string
  var sys_id = ee.String(image_id.get("system:id")).getInfo(); 
  var filename = sys_id.slice(52,100); // extract the fileName (NEON domain + site code + product code + year)

  // Dynamically determine the range of data to display
  // Sets color scale to show all but lowest/highest 2% of data
  var pct_clip = image_id.reduceRegion({
    reducer: ee.Reducer.percentile([2, 98]),
    scale: 10,
    maxPixels: 3e7});

  var keys = pct_clip.keys();
  var pct02 = ee.Number(pct_clip.get(keys.get(0))).round().getInfo()
  var pct98 = ee.Number(pct_clip.get(keys.get(1))).round().getInfo()
    
  Map.addLayer(image_id, {min:pct02, max:pct98, palette: chm_palette}, filename, 0)
}

// Call the addSingleBandImage function to add CHM layers to map 
grsm_chm_col.evaluate(function(grsm_chm_col) {
  grsm_chm_col.features.map(addSingleBandImage);
})

// Center the map on GRSM and set zoom level to 12
Map.setCenter(-83.5, 35.6, 12);

Now that you've read in these two datasets (SDR and CHM) over all the years of available data, we encourage you to explore the different layers and see what you notice! Toggle between the layers, play with the opacity. Visual inspection is an important first step in exploratory analysis - see if you can recognize patterns and form new questions based off what you see.

CHM Difference Layers

Next let's create a new raster layer of the difference between the CHMs from 2 different years.

// Difference the CHMs from 2017 and 2016 and 2021
var grsm_chm2021 = grsm_chm_col.filterDate('2021-01-01', '2021-12-31').first();
var grsm_chm2017 = grsm_chm_col.filterDate('2017-01-01', '2017-12-31').first();
var grsm_chm2016 = grsm_chm_col.filterDate('2016-01-01', '2016-12-31').first();

// Subtract the CHMs to create difference CHM rasters
var chm_diff_2017_2016 = grsm_chm2017.subtract(grsm_chm2016);
var chm_diff_2021_2017 = grsm_chm2021.subtract(grsm_chm2017);
var chm_diff_2021_2016 = grsm_chm2021.subtract(grsm_chm2016);

// Display the first CHM difference raster (2017-2016) and add as a layer to the Map
print('CHM Difference 2017-2016',chm_diff_2017_2016)

Map.addLayer(chm_diff_2017_2016, {min: -10, max: 10, palette: dchm_palette}, 'CHM diff 2017-2016');

CHM Difference Stats and Histograms

Next let's calculate the mean difference in Canopy Height inside the fire perimeter for the various years. We'll also plot histograms of the CHM differences.

// Calculate the mean dCHM between the various years:
print('Mean dCHM in the Chimney Tops Fire Perimeter')

print('Mean dCHM 2017-2016',chm_diff_2017_2016.reduceRegion({
      reducer: ee.Reducer.mean(),
      geometry: ct_fire_boundary,
      scale: 30}));
      
print('Mean dCHM 2021-2017',chm_diff_2021_2017.reduceRegion({
      reducer: ee.Reducer.mean(),
      geometry: ct_fire_boundary,
      scale: 30}));
      
print('Mean dCHM 2021-2016',chm_diff_2021_2016.reduceRegion({
      reducer: ee.Reducer.mean(),
      geometry: ct_fire_boundary,
      scale: 30}));

In the console, if you expand the objects, you can see that from 2016-2017, there was a net loss in canopy height of ~6.6m, and between 2017-2021 there was a net growth of ~3m, suggesting a considerable amount of re-growth in the 5 years after the fire.

We can also look at the histograms of the CHM differences to provide a little more information about the ecosystem structure dynamics immediately after the fire and in the subsequent years. To plot histograms, first write a function to create a histogram given a difference CHM raster (calculated above) and a string of the years that are being differenced, as inputs. The string is just used to include in the histogram chart title.

// Function to create histogram charts for each CHM difference layer, clipped by the chimney tops fire perimeter
function chmDiffHist(img,years_str) { 
  var hist =
      ui.Chart.image.histogram({image: img.clip(ct_fire_boundary), region: ct_fire_boundary, scale: 50})
          .setOptions({title: 'CHM Difference Histogram ' + years_str,
                      hAxis: {title: 'CHM Difference (m)',titleTextStyle: {italic: false, bold: true},},
                      vAxis: {title: 'Count', titleTextStyle: {italic: false, bold: true}},});
  return hist
}

// Apply the function to the three CHM difference rasters
var chm_diff_hist_2017_2016 = chmDiffHist(chm_diff_2017_2016,'2017-2016')
var chm_diff_hist_2021_2016 = chmDiffHist(chm_diff_2021_2016,'2021-2016')
var chm_diff_hist_2021_2017 = chmDiffHist(chm_diff_2021_2017,'2021-2017')

// Display the CHM difference histograms charts on the Console
print(chm_diff_hist_2017_2016);
print(chm_diff_hist_2021_2017);
print(chm_diff_hist_2021_2016);

On your own, try to interpret what these difference histograms are showing.

Normalized Burn Ratio (NBR)

Last but not least, we can take a quick look at the NBR and dNBR. Refer to the CU Earth Lab dNBR Lesson for a nice explanation of this metric in the context of multispectral satellite data.

// Read in clear SDR images at GRSM in 2016, 2017, and 2021 
var grsm_sdr2016_clear = grsm_sdr_cloudfree.filterDate('2016-01-01', '2016-12-31').first()
var grsm_sdr2017_clear = grsm_sdr_cloudfree.filterDate('2017-01-01', '2017-12-31').first();
var grsm_sdr2021_clear = grsm_sdr_cloudfree.filterDate('2021-01-01', '2021-12-31').first();
  
//------------------------- Normalized Difference Burn Ratio ----------------------------
// The normalized burn ratio (NBR) is a normalized difference index using the shortwave-infrared (SWIR) and near-infrared (NIR) portions of the electromagnetic spectrum. dNBR can be used as a metric to map fire extent and burn severity when calculating the difference between pre and post fire conditions.

// calculate NBR for the 3 years
// B097: B365: 
var sdr_pre_nbr_2016 = grsm_sdr2016_clear.normalizedDifference(['B097', 'B365']);
var sdr_post_nbr_2017 = grsm_sdr2017_clear.normalizedDifference(['B097', 'B365']);
var sdr_post_nbr_2021 = grsm_sdr2021_clear.normalizedDifference(['B097', 'B365']);

// calculate dNBR 2016-2017 and 2016-2021
var sdr_dNBR_2016_2017 = sdr_pre_nbr_2016
                            .subtract(sdr_post_nbr_2017)
                            .clip(ct_fire_boundary);

var sdr_dNBR_2016_2021 = sdr_pre_nbr_2016
                            .subtract(sdr_post_nbr_2021)
                            .clip(ct_fire_boundary);

// Remove comment-symbols (//) below to display pre- and post-fire NBR as layers
// Map.addLayer(sdr_pre_nbr_2016, {min: -1, max: 1, palette: red_ylw_grn}, 'Pre-fire (June 2016) Normalized Burn Ratio');
// Map.addLayer(sdr_post_nbr_2017, {min: -1, max: 1, palette: red_ylw_grn}, 'Post-fire (Oct 2017) Normalized Burn Ratio');
// Map.addLayer(sdr_post_nbr_2021, {min: -1, max: 1, palette: red_ylw_grn}, 'Post-fire (June 2021) Normalized Burn Ratio');

// add dNBR layers
Map.addLayer(sdr_dNBR_2016_2017, {min: -1, max: 1, palette: dnbr_palette}, 'dNBR 2016-2017');
Map.addLayer(sdr_dNBR_2016_2021, {min: -1, max: 1, palette: dnbr_palette}, 'dNBR 2016-2021');

The differenced Normalized Burn Ratio (dNBR) does a great job of highlighting the burned areas (in red).

On your own, we encourage you to dig into the code from this tutorial and expand upon it according to your scientific interests. Think of some questions you have about this dataset and think about how you might answer it using GEE. Modify these functions or try writing your own function to answer your question(s). For example, try out different reducers to compile other statistis to summarize the CHM and NBR differences, or see if there are any other datasets that you could bring in to expand your analysis. This is just the starting point!

Get Lesson Code

GRSM Wildfire Change Exploration

In this lesson, we'll continue to use the Great Smokey Mountains site as an example, this time creating a time series of the mean NDVI within the Chimney Tops 2 Fire perimeter between 2016-2022. We will plot this along with the NDVI time-series derived from Landsat 8 data, and use this to fill in some more detailed temporal information.

AOP strives to collect every site during peak-greenness, when the predominant vegetation is most photosynthetically active. This is so that when comparing data from year to year, the differences are due to actual changes and not just due to the time of year. This is not always possible, so it's important to consider the time of year when you are conducting your analysis.

Read in AOP and Landsat 8 Surface Reflectance Image Collections

First read in the AOP SDR data and the Landsat 8 data, filtering by the site center point / region of interest, and by date.

// Specify center location and for GRSM
var site_center = ee.Geometry.Point([-83.5, 35.7]); 

// Create region of interest (roi)
var roi = ee.FeatureCollection('projects/neon-sandbox-dataflow-ee/assets/chimney_tops_fire')

// Read in the SDR Image Collection
var sdr_col = ee.ImageCollection('projects/neon-prod-earthengine/assets/HSI_REFL/001')
  .filterBounds(site_center);

// Read in Landsat 8 Surface Reflectance Image Collection
var l8sr = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
                .filterBounds(roi)
                .filter(ee.Filter.calendarRange(2016, 2022, 'year'));

Cloud Mask Function for Landsat 8 Data

Next we can create a function to pre-process the Landsat 8 data, which applies scaling and cloud / saturated data masking. Don't worry too much about the details here, but the main thing to be aware of is that different satellite image collections handle the QA information differently. Landsat (and other satellites) often use something called "bitmasking" to store QA information, which is a space-efficient method for storing additional information. This cloud-masking function can be found on the earthengine-api GitHub examples here.

// cloud masking function for Landsat 8 Collection 2 
function maskL8sr(image) {

  var qaMask = image.select('QA_PIXEL').bitwiseAnd(parseInt('11111', 2)).eq(0);
  var saturationMask = image.select('QA_RADSAT').eq(0);

  // Apply the scaling factors to the appropriate bands.
  var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
  var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0);

  // Replace the original bands with the scaled ones and apply the masks.
  return image.addBands(opticalBands, null, true)
      .addBands(thermalBands, null, true)
      .updateMask(qaMask)
      .updateMask(saturationMask);
}

// Apply the cloud masking function
l8sr = l8sr.filterBounds(roi).map(maskL8sr)

Merge AOP and Landsat 8 NDVI Collections

Next we can plot the two datasets on the same chart. This code was modified from the Stack Overflow post: how-to-combine-time-series-datasets-with-different-timesteps-in-a-single-plot-on.

// compute ndvi bands to add to each collection
var addL8Bands = function(image){
  var l8_ndvi = image.normalizedDifference(['SR_B5','SR_B4']).rename('l8_ndvi')
  var aop_ndvi = ee.Image().rename('aop_ndvi') 
  return image.addBands(l8_ndvi).addBands(aop_ndvi)
}

var addAOPBands = function(image){
  var aop_ndvi = image.normalizedDifference(['B097', 'B055']).rename('aop_ndvi')
  var l8_ndvi = ee.Image().rename('l8_ndvi') 
  return image.addBands(aop_ndvi).addBands(l8_ndvi)
}

l8sr = l8sr.map(addL8Bands)
sdr_col = sdr_col.map(addAOPBands)

print('NIS Images',sdr_col)

// merge the collections
var merged = l8sr.merge(sdr_col).select(['l8_ndvi', 'aop_ndvi'])

Plot NDVI Time Series

Lastly we can create and plot (print) the time-series chart. Most of this is just setting the chart style.

// Set chart style properties.
// https://developers.google.com/earth-engine/guides/charts_style
var chartStyle = {
  title: 'NDVI at Chimney Tops Fire ROI - AOP + Landsat 8',
  hAxis: {
    title: 'Date',
    titleTextStyle: {italic: false, bold: true},
    gridlines: {color: 'FFFFFF'}
  },
  vAxis: {
    title: 'Mean NDVI',
    titleTextStyle: {italic: false, bold: true},
    gridlines: {color: 'FFFFFF'},
    format: 'short',
    baselineColor: 'FFFFFF'
  },
  series: {
    0: {lineWidth: 3, color: 'E37D05', pointSize: 7},
    1: {lineWidth: 3, color: '1D6B99'}
  },
  chartArea: {backgroundColor: 'EBEBEB'}
};

// Plot the merged (AOP + Landsat 8) Image Collection NDVI Time Series
var ndvi_timeseries = ui.Chart.image.series({
  imageCollection: merged,
  region: roi,
  reducer: ee.Reducer.mean(),
  scale: 30 // 	Scale to use with the reducer in meters
}).setOptions(chartStyle);

print(ndvi_timeseries)

This figure demonstrates some important points. First, we can see that there is good alignment between the mean NDVI calculated from AOP reflectance data and Landsat 8 at these times. While the Airborne Observation Platform seeks to collect data in peak green conditions, this is not always possible due to logistical or other constraints. For this site, one of the AOP collections in 2017 was in October, past peak-greenness, and the leaves already started senescing in some areas. Generating a time-series plot like this can help highlight the data in the context of the larger temporal trends. In this example, the NDVI time-series generated from Landsat 8 highlights strong seasonal trends. This brings us to the second point: while AOP data has very high spectral (426 bands) and spatial (1m) resolution, the temporal resolution (annual or bi-annual) may not fully suffice for your research needs. This is where scaling with satellite data - either to expand your analysis to a larger area, or to achieve a more comprehensive temporal understanding - can be very informative. This tutorial demonstrates a simple example of scaling up the temporal resolution, but demonstrates how GEE makes this sort of scalable analysis much simpler!

Get Lesson Code

GRSM NDVI Time Series - Landsat 8 & AOP