AOP GEE Data Availability & Access

AOP has published a subset of AOP Level 3 (mosaicked) data products at 5 NEON sites (as of Spring 2022) on GEE. This data has been converted to Cloud Optimized GeoTIFF (COG) format. NEON L3 lidar and derived spectral indices are available in geotiff raster format, so are relatively straightforward to add to GEE, however the hyperspectral data is available in hdf5 (hierarchical data) format, and have been converted to the COG format prior to being added to GEE.

The NEON data products that have been made available on GEE can be accessed through the projects/neon folder with an appended prefix of the Data Product ID, matching the NEON data portal. The table below summarizes the prefixes to use for each data product, and is a useful reference for reading in AOP GEE datasets. You will see how to access and read in these data products in the next part of this lesson.

The table below summarizes the sites, products, and years of NEON AOP data that can currently be accessed in GEE.

Get Started with Google Earth Engine

Once you have set up your Google Earth Engine account you can navigate to the Earth Engine Code Editor. The diagram below, from the Earth-Engine Playground, shows the main components of the code editor. If you have used other programming languages such as R, Python, or Matlab, this should look fairly similar to other Integrated Development Environments (IDEs) you may have worked with. The main difference is that this has an interactive map at the bottom, similar to Google Maps and Google Earth. This editor is fairly intuitive. We encourage you to play around with the interactive map, or explore the ee documentation, linked above, to gain familiarity with the various features.

Read AOP Data Collections into GEE using ee.ImageCollection

AOP data can be accessed through GEE through the projects/neon folder. In the remainder of this lesson, we will look at the three AOP datasets, or ImageCollections in this folder.

An ImageCollection is simply a group of images. To find publicly available datasets (primarily satellite data), you can explore the Earth Engine Data Catalog. Currently, NEON AOP data cannot be discovered in the main GEE data catalog, so the following steps will walk you through how to find available AOP data.

In your code editor, copy and run the following lines of code to create 3 ImageCollection variables containing the Surface Directional Reflectance (SDR), Camera Imagery (RGB) and Digital Surface and Terrain Model (DEM) raster data sets.

//read in the AOP image collections as variables

var aopSDR = ee.ImageCollection('projects/neon/DP3-30006-001_SDR')

var aopRGB = ee.ImageCollection('projects/neon/DP3-30010-001_RGB') 

var aopDEM = ee.ImageCollection('projects/neon/DP3-30024-001_DEM')

A few tips for the Code Editor:

• In the left panel of the code editor, there is a Docs tab which includes API documentation on built in functions, showing the expected input arguments. We encourage you to refer to this documentation, as well as the GEE JavaScript Tutorial to familiarize yourself with GEE and the JavaScript programming language.
• If you have an error in your code, a red error message will show up in the Console (in the right panel), which tells you the line that failed.
• Save your code frequently! If you try to leave your code while it is unsaved, you will be prompted that there are unsaved changes in the editor.

When you Run the code above (by clicking on the Run above the code editor), you will notice that the lines of code become underlined in red, the same as you would see for a spelling error in most text editors. If you hover over each of the lines of codes, you will see a message pop up that says: <variable> can be converted to an import record. Convert Ignore.

If you click Convert, the line of code will disappear and the variable will be imported into your session directly, and will show up at the top of the code editor. Go ahead and convert the variables for all three lines of code, so you should see the following. Tip: if you type Ctrl-z, you can re-generate the line of code, and the variable will still show up in the imported variables at the top of the editor. It is a good idea to retain the original code that reads in the variable, for reproducibility. If you don't do this, and wish to share this code with someone else, or run the code outside of your own code editor, the imported variables will not be saved.

Note that each of these imported variables can now be expanded, using the arrow to the left of each. These variables now show associated information including type, id, and properties, which if you expand, shows a description. This provides more detailed information about the data product.

Information about the image collections can also be found in a slightly more user-friendly format if you click on the blue projects/neon/DP3-30006-001_SDR, as well as DP3-30010-001_RGB andDP3-30024-001_DEM, respectively. Below we'll show the window that pops-up when you click on SDR, but we encourage you to look at all three datasets.

This allows you to read the full description in a more user-friendly format. Note that the images imported into GEE may have some slight differences from the data downloaded from the data portal. For example, note that the reflectance data in GEE is scaled by 100. We highly encourage you to explore the description and associated documentation for the data products on the NEON data portal as well (eg. DP3.30006.001) for relevant information about the data products, how they are generated, and other pertinent details.

You can also click on the IMAGES tab to explore all the available NEON images for that data product. Some of the text may be cut off in the default view, but if you click in one of the table values the table will expand. This table summarizes individual sites and years that are available for the SDR Image Collection. The ImageID provides the path to read in an individual image. In the next step, we will show how to use this path to pull in a single file.

Read AOP Data into GEE using ee.Image

As a last step, we will go ahead and use the path specified in the SDR Asset Details Images table to read in a single image. Pulling in a single image uses almost identical syntax as an image collection, see below:

var TALL_2017_SDR = ee.Image('projects/neon/DP3-30006-001_SDR/DP3-30006-001_D08_TALL_SDR_2017')

Import this variable, and you can see that it pulls in to the Imports at the top, and shows (426 bands) at the right. To the right of that you will see blue eye and target icons. If you hover over the eye it displays "Show on Map". Click this eye icon to place a footprint of this data set in the Map display. If you hover over the target icon, you will see the option "Center Map on Record". Click this to center your map on this TALL SDR dataset. You should now see the footprint of the data as a layer in the Google Map.

A Quick Recap

You did it! You should now have a basic understanding of the GEE code editor and it's different components. You have also learned how to read a NEON AOP ImageCollection into a variable, import the variable into your code editor session, and navigate through the ImageCollection Asset details to find the path to an individual Image. Lastly, you learned to read in an individual SDR Image, pull the footprint of the data into a Map Layer, and center on that region.

It doesn't look like we've done much so far, but this is a already great achievement! With just a few lines of code, you can import an entire AOP hyperspectral data set, which in most other coding environments, is not simple. One of the barriers to working with AOP data (and reflectance data in particular) is it's large data volume, which requires high-performance computing environments to carry out analysis. There are also limited open-source tools for working with the data; many of the software suites for working with hyperspectral data require licenses which can be expensive. In this lesson, we have loaded spectral data covering an entire site, and are ready for data exploration and analysis, in a free geospatial cloud-computing platform.

In the next tutorials, we will pull in spectral data, visualize RGB and false color image composites, interactively plot spectral signatures of pixels in the image, and carry out some more advanced analysis that is highly simplified by the built in GEE functions.

Get Lesson Code

Into to AOP GEE Assets

Read in and Visualize AOP SDR Data

In the first Intro to AOP data in GEE tutorial, we showed how to explore the NEON AOP GEE Image Collections. We will build off that tutorial in this lesson, to pull in and visualize some AOP hyperspectral data in GEE. Specifically, we will look at surface directional reflectance (SDR) data collected at the NEON site SRER (Santa Rita Experimental Range) for 3 years between 2018 and 2021.

Let's get started! In this tutorial we generate basic GEE (JavaScript) code to visualize hyperspectral data. We will work through the following steps:

1. Pull in an AOP hyperspectral (SDR) image
2. Set the visualization parameters
3. Mask the no data values
4. Add the AOP SDR layer to the GEE Map
5. Center on the region of interest and set zoom level

We encourage you to follow along with this code chunks in this exercise in your code editor. To run the cells, you can click the Run button at the top of the code editor. Note that until you run the last two steps (adding the data layer to the map), you will not see the AOP data show up in the Interactive Map.

Read in the SRER 2018 SDR image

Using ee.Image, you can pull in a single image if you know the path and name of the image (as opposed to filtering down to the individual images from an image collection). If you don't know this path, you can pull in the image collection and look at the Asset Details > Images tab. We will assign this image to a variable (var) called SRER_SDR2018. You can refer to the tables in the Data Access and Availability section, in the Intro to AOP data in GEE tutorial, to see how to pull in spectral data from a different site or date.

var SRER_SDR2018 = ee.Image("projects/neon/D14_SRER/L3/DP3-30006-001_D14_SRER_SDR_2018");

As we covered in the previous lesson, when you type this code, it will be underlined in red (the same as you would see with a mis-spelled word). When you hover over this line, you will see an option pop up that says "SRER_SDR2018" can be converted to an import record. Convert Ignore

If you click Convert, the line of code will disappear and the variable will be pulled into the top of the code editor, as shown below. Once imported, you can interactively explore this variable - eg. you can expand on the bands and properties to gain more information about this image, or "asset", as it's called in GEE.

Another way to learn more about this asset is to left-click on the blue projects/neon/D14_SRER/L3_DP3-30006-001-D14_SRER_SDR_2018. This will pop up a box with more detailed information about this asset, as shown below:

Click Esc to return to the code editor. Note that you can run the code either way, with the variable explicitly specified in the code editor, or imported as a variable, but we encourage you to leave the variable written out in the code, as this way is more reproducible.

Set the visualization parameters

The visualization parameters specifies the band combination that is displayed, and other display options, such as the minimum and maximum values for the histogram stretch. For more detailed information, refer to the GEE documentation on image visualization.

To set the visualization parameters, we will create a new variable (called visParams). This variable is applied to the layer and determines what is displayed. In this we are setting the RGB bands to display - for this exercise we are setting them to red, green, and blue portions of the spectrum in order to show a True Color Image. You can change these bands to show a False Color Image or any band combination of interest. You can refer to NEON's lessons on Multi-Band Rasters in R or RGB and False Color Images in Python for more background on band stacking.

var visParams = {'min':2,'max':20,'gamma':0.9,'bands':['band053','band035','band019']};

Mask the no data values

This step is optional, but recommended. AOP sets No Data Values to -9999, so if you don't mask these out you will see any missing data as black in the image (this will often result in a black boundary surrounding the site, but if any data is missing inside the site that will show up as black as well). To show only the data that was collected, we recommend masking these values using the updateMask function, keeping only values greater than or equal to zero, as shown below:

var SRER_SDR2018mask = SRER_SDR2018.updateMask(SRER_SDR2018.gte(0.0000));

Add SDR layer to the map

Now that we've defined the data, the visualization parameters, and the mask, we can add the reflectance layer to the Map! To do this, we use the Map.addLayer function with our masked data variable, SRER_SDRmask, using the visParams and assign this layer a label, which will show up in the Map.

Map.addLayer(SRER_SDR2018mask, visParams, 'SRER 2018');

Center the map on our area of interest and set zoom level

GEE by default does not know where we are interested in looking. We can center the map over our new data layer by specifiying Map.setCenter with the longitude, latitude, and zoom level (for this site, zoom of 11 works well to show the full site, but you can try other values to see how the image size changes).

Map.setCenter(-110.83549, 31.91068, 11);

Putting it All Together

The following code chunk runs all the steps we just broke down, and also adds in 2 more years of data (2019 and 2021). You can pull in this code into your code editor by clicking here, or alternately copy and paste the code below into your GEE code editor. Click Run to add the 3 SDR data layers for each year.

// This script pulls in hyperspectral data over the Santa Rita Experimental Range (SRER)
// from 2018, 2019, and 2021 and plots RGB 3-band composites

// Read in Surface Directional Reflectance (SDR) Images 
var SRER_SDR2018 = ee.Image("projects/neon/D14_SRER/L3/DP3-30006-001_D14_SRER_SDR_2018");
var SRER_SDR2019 = ee.Image("projects/neon/D14_SRER/L3/DP3-30006-001_D14_SRER_SDR_2019");
var SRER_SDR2021 = ee.Image("projects/neon/D14_SRER/L3/DP3-30006-001_D14_SRER_SDR_2021");

// Set the visualization parameters so contrast is maximized, and set display to show RGB bands 
var visParams = {'min':2,'max':20,'gamma':0.9,'bands':['band053','band035','band019']};

// Mask layers to only show values > 0 (this hides the no data values of -9999) 
var SRER_SDR2018mask = SRER_SDR2018.updateMask(SRER_SDR2018.gte(0.0000));
var SRER_SDR2019mask = SRER_SDR2019.updateMask(SRER_SDR2019.gte(0.0000));
var SRER_SDR2021mask = SRER_SDR2021.updateMask(SRER_SDR2021.gte(0.0000));

// Add the 3 years of SRER SDR data as layers to the Map:
Map.addLayer(SRER_SDR2018mask, visParams, 'SRER 2018');
Map.addLayer(SRER_SDR2019mask, visParams, 'SRER 2019');
Map.addLayer(SRER_SDR2021mask, visParams, 'SRER 2021');

// Center the map on SRER & zoom to desired level (11 = zoom level)
Map.setCenter(-110.83549, 31.91068, 11);

Once you have the three years of data added, you can look at the different years one at a time by selecting each layer in the Layers box inside the Map:

If you click anywhere inside the AOP map (where there is data), you will see the 426 spectral bands as a bar chart displayed for each of the layers in the Inspector window (top-right corner of the code editor). You can see the spectral values for different layers by clicking on the arrow to the left of the layer name under Pixels (eg. SRER 2018). Note that these values are just shown as band #s, and you can't tell from the chart what the actual wavelength values are. We will convert the band numbers to wavelengths in the next lesson, so stay tuned!

Get Lesson Code

Importing and Visualizing SRER SDR Data

Writing a Function to Visualize AOP SDR Image Collections

In the previous Introduction to AOP Hyperspectral Data in GEE tutorial, we showed how to read in SDR data for images from 3 years. In this tutorial, we will show you a different, more simplified way of doing the same thing, using functions. 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.

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 arguments args and returns an output.

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

To call the function for a full image collection, you can use a map to iterate over items in a collection. This is shown in the script below.

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

For this example, we will 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 print statements in this function below, and show the outputs (which would show up in the console tab). Note that these print statements are commented out in the code linked with this tutorial, since they are not required for the function to run.

// Specify center location of SRER
var mySite = ee.Geometry.Point([-110.83549, 31.91068])

// Read in the SDR Image Collection
var NISimages = ee.ImageCollection('projects/neon/DP3-30006-001_SDR').filterBounds(mySite)

// Create a 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()
  print("systemID: "+sysID_serverStr)
  // truncate the string to show only the fileName (NEON domain + site code + product code + year)
  var fileName = sysID_serverStr.slice(46,100); 
  print("fileName: "+fileName)
  // mask out no-data values and set visualization parameters to show RGB composite
  var imageRGB = imageId.updateMask(imageId.gte(0.0000)).select(['band053', 'band035', 'band019']);
  // add this layer to the map
  Map.addLayer(imageRGB, {min:2, max:20}, fileName)
}

// call the addNISimages function
NISimages.evaluate(function(NISimages) {
  NISimages.features.map(addNISImage);
})

// Center the map on SRER
Map.setCenter(-110.83549, 31.91068, 11);

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. You should recognize some of the same syntax from the previous tutorial in the last two lines of code in the function, defining the variable imageRGB, using updateMask, and finally using Map.addLayer to add the layer to the Map window. Note that this function is subsetting the SDR image to only pull in the red, green, and blue bands, as opposed to the previous tutorial where we read in the full hyperspectral cube, and then displayed only the RGB composite in the visParam variable.

You can see that the print statements are showing up in the console, displaying the systemID and fileName for each image in the collection. The fileName is applied to the name of the layers in the Map window.

You could alter this function to include the visualization paramters, to subset by other bands, or modify it to work for a different image collection. We encourage you to do this on your own!

Get Lesson Code

Function to display AOP SDR Image Collections in GEE

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) datasets, as well as the Ecosystem Structure (CHM) data to look at interannual differences at the site SRER. We will discuss some of the acquisition parameters and other factors that affect data quality and interperatation.

Functions to Read in SDR and CHM Image Collections

Let's get started. In this first chunk of code, we are setting the center location of SRER, reading in the SDR image collection, and then creating a function to display the NIS image in GEE. Optionally, we've included lines of code that also calculate and display NDVI from the SDR data, so you can uncomment this and run on your own, if you wish. For more details on how this function works, you can refer to the tutorial Intro to GEE Functions.

// Specify center location of SRER
var mySite = ee.Geometry.Point([-110.83549, 31.91068])

// Read in the SDR Image Collection
var NISimages = ee.ImageCollection('projects/neon/DP3-30006-001_SDR')
  .filterBounds(mySite)

// Function to display NIS Image + NDVI
function addNISImage(image) { 
  var imageId = ee.Image(image.id); // get the system:id and convert to string
  var sysID = ee.String(imageId.get("system:id")).getInfo();  // get the system:id - this is an object on the server
  var fileName = sysID.slice(46,100); // extract the fileName (NEON domain + site code + product code + year)
  var imageMasked = imageId.updateMask(imageId.gte(0.0000)) // mask out no-data values
  var imageRGB = imageMasked.select(['band053', 'band035', 'band019']);  // select only RGB bands for display
  
  Map.addLayer(imageRGB, {min:2, max:20}, fileName, 0)   // add RGB composite to the map, 0 means it is deselected initially
  
  // uncomment the lines below to calculate NDVI from the SDR and display NDVI
  // var nisNDVI =imageMasked.normalizedDifference(["band097", "band055"]).rename("ndvi") // band097 = NIR , band055 = red
  // var ndviViz= {min:0.05, max:.95,palette:['brown','white','green']}
  
  // Map.addLayer(nisNDVI, ndviViz, "NDVI "+fileName, 0) // optionally add NDVI to the map
}

// call the addNISimages function to add SDR and NDVI layers to map
NISimages.evaluate(function(NISimages) {
  NISimages.features.map(addNISImage);
})

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. Note that we can use the .filterMetadata property to select only the CHM data from the DEM image collection, since the CHM is stored in that collection, along with the DTM and DSM.

// Read in only the CHM Images (using .filterMetadata by Type)
var CHMimages =  ee.ImageCollection('projects/neon/DP3-30024-001_DEM')
  .filterBounds(mySite)
  .filterMetadata('Type','equals','CHM')

// Function to display Single Band Images setting display range to linear 2%
function addSingleBandImage(image) { // display each image in collection
  var imageId = ee.Image(image.id); // get the system:id and convert to string
  var sysID = ee.String(imageId.get("system:id")).getInfo(); 
  var fileName = sysID.slice(46,100); // extract the fileName (NEON domain + site code + product code + year)
  // print(fileName) // optionally print the filename, this will be the name of the layer
  
  // dynamically determine the range of data to display
  // sets color scale to show all but lowest/highest 2% of data
  var pctClip = imageId.reduceRegion({
    reducer: ee.Reducer.percentile([2, 98]),
    scale: 10,
    maxPixels: 3e7});

  var keys = pctClip.keys();
  var pct02 = ee.Number(pctClip.get(keys.get(0))).round().getInfo()
  var pct98 = ee.Number(pctClip.get(keys.get(1))).round().getInfo()
    
  var imageDisplay = imageId.updateMask(imageId.gte(0.0000));
  Map.addLayer(imageDisplay, {min:pct02, max:pct98}, fileName, 0)
// print(image)
}

// call the addSingleBandImage function to add CHM layers to map 
// you could also add all DEM images (DSM/DTM/CHM) but for now let's just add CHM
CHMimages.evaluate(function(CHMimages) {
  CHMimages.features.map(addSingleBandImage);
})

// Center the map on SRER and set zoom level
Map.setCenter(-110.83549, 31.91068, 11);

Now that you've read in these two datasets over all the years, we encourage you to explore the different layers and see if you notice any patterns!

Creating CHM Difference Layers

Next let's create a new raster layer of the difference between the CHMs from 2 different years. We will difference the CHMs from 2018 and 2021 because these years were both collected with the Riegl Q780 system and so have a vertical resolution (CHM height cutoff) of 2/3 m. By contrast the Gemini system (which was used in 2017 and 2020) has a 2m cutoff, so some of the smaller shrubs are not resolved with that sensor. It is important to be aware of factors such as these that may affect the interpretation of the data! We encourage all AOP data users to read the associated metadata pdf documents that are provided with the data products (when downloading from the data portal or using the API).

For more information on the vertical resolution, read the footnotes at the end of this lesson.

var SRER_CHM2018 = ee.ImageCollection('projects/neon/DP3-30024-001_DEM')
  .filterDate('2018-01-01', '2018-12-31')
  .filterBounds(mySite).first(); 

var SRER_CHM2021 = ee.ImageCollection('projects/neon/DP3-30024-001_DEM')
  .filterDate('2021-01-01', '2021-12-31')
  .filterBounds(mySite).first();

var CHMdiff_2021_2018 = SRER_CHM2021.subtract(SRER_CHM2018);
// print(CHMdiff_2021_2018)
Map.addLayer(CHMdiff_2021_2018, {min: -1, max: 1, palette: ['#FF0000','#FFFFFF','#008000']},'CHM diff 2021-2018')

You can see some broad differences, but there also appear to be some noisy artifacts. We can smooth out some of this noise by using a spatial filter.

// Smooth out the difference raster (filter out high-frequency patterns)
// Define a boxcar or low-pass kernel.
var boxcar = ee.Kernel.square({
  radius: 1.5, units: 'pixels', normalize: true
});

// Smooth the image by convolving with the boxcar kernel.
var smooth = CHMdiff_2021_2018.convolve(boxcar);
Map.addLayer(smooth, {min: -1, max: 1, palette: ['#FF0000','#FFFFFF','#008000']}, 'CHM diff, smoothed');

CHM Difference Histograms

Next let's plot histograms of the CHM differences, between 2021-2018 as well as between 2021-2019 and 2019-2018. For this example, we'll just look at the values over a small area of the site. Looking at these 3 sets of years, we will see some of the artifacts related to the lidar sensor used (Riegl Q780 or Optech Gemini). If you didn't know about the differences between the sensors, it would look like the canopy was growing and shrinking from year to year.

Before running this chunk of code, you'll need to create a polygon of a region of interest. For this example, I selected a region in the center of the map, shown below, although you can select any region within the site. To create the polygon, select the rectangle out of the shapes in the upper left corner of the map window (hovering over it should say "Draw a rectangle"). Then drag the cursor over the area you wish to cover.

When you load in this geometry, you should see the geometry variable imported at the top of the code editor window, under Imports. It should look something like this:

Once you have this geometry polygon variable set, you can run the following code to generate histograms of the CHM differences over this area.

// read in CHMs from 2019 and 2017
var SRER_CHM2019 = ee.ImageCollection('projects/neon/DP3-30024-001_DEM')
  .filterDate('2019-01-01', '2019-12-31')
  .filterBounds(mySite).first();

var SRER_CHM2017 = ee.ImageCollection('projects/neon/DP3-30024-001_DEM')
  .filterDate('2017-01-01', '2017-12-31')
  .filterBounds(mySite).first();

// calculate the CHM difference histograms (2021-2019 & 2019-2018)
var CHMdiff_2021_2019 = SRER_CHM2021.subtract(SRER_CHM2019);
var CHMdiff_2019_2018 = SRER_CHM2019.subtract(SRER_CHM2018);

// Define the histogram charts for each CHM difference image, print to the console.
var hist1 =
    ui.Chart.image.histogram({image: CHMdiff_2021_2019, region: geometry, scale: 50})
        .setOptions({title: 'CHM Difference Histogram, 2021-2019',
                    hAxis: {title: 'CHM Difference (m)',titleTextStyle: {italic: false, bold: true},},
                    vAxis: {title: 'Count', titleTextStyle: {italic: false, bold: true}},});
print(hist1);

var hist2 =
    ui.Chart.image.histogram({image: CHMdiff_2019_2018, region: geometry, scale: 50})
        .setOptions({title: 'CHM Difference Histogram, 2019-2018',
                    hAxis: {title: 'CHM Difference (m)',titleTextStyle: {italic: false, bold: true},},
                    vAxis: {title: 'Count', titleTextStyle: {italic: false, bold: true}},});
print(hist2);

var hist3 =
    ui.Chart.image.histogram({image: CHMdiff_2021_2018, region: geometry, scale: 50})
        .setOptions({title: 'CHM Difference Histogram, 2021-2018',
                    hAxis: {title: 'CHM Difference (m)',titleTextStyle: {italic: false, bold: true},},
                    vAxis: {title: 'Count', titleTextStyle: {italic: false, bold: true}},});
print(hist3);

Let's take a minute to understand what's going on here. In each case, we subtracted the earlier year from the later year. So from 2019 to 2021, it looks like the vegetation grew on average by ~0.6m, but from 2018 to 2019 it shrunk by the same amount. This is because in 2021 there was a lower vertical cutoff, so shrubs of at least 0.67m were resolved, where before anything below 2m was obscured. These low shrubs are likely the dominant source of the change we're seeing. We can see the same pattern, but in reverse between 2018 and 2019. The difference histogram from 2021 to 2018 more accurately represents the change, which is centered around 0, and the map we displayed shows local changes in certain areas, related to actual vegetation growth and ecological drivers. Note that 2021 was a particularly wet year, and AOP's flight was in optimal peak greenness, as you can see when comparing the SDR imagery to earlier years.

NDVI Time Series

Last but not least, we can take a quick look at NDVI changes over the four years of data. A quick way to look at the interannual changes are to make a line plot, which we'll do shortly. First let's take a step back and see the weather conditions during the collections. For every mission, the AOP flight operators assess the cloud conditions and note whether the cloud clover is <10% (green), 10-50% (yellow), or >50% (red). This information gets passed through to the reflectance hdf5 data, and is also available in the summary metadata documents, delivered with all the spectrometer data products. The weather conditions have direct implications for data quality, and while we strive to collect data in "green" weather conditions, it is not always possible, so the user must take this into consideration when working with the data.

The figure below shows the weather conditions at SRER for each of the 4 collections. In 2017 and 2021, the full site was collected in <10% cloud conditions, while in 2018 and 2019 there were mixed weather conditions. However, for all four years, the center of the site was collected in optimal cloud conditions.

With this in mind, let's use the same geometry we used before, centered in the middle of the plot, to look at the mean NDVI over the four years in a small part of the site. Here's the GEE code for doing this:

// calculate NDVI for the geometry
var ndvi = NISimages.map(function(image) {
    var ndviClip = image.clip(geometry)
    return ndviClip.addBands(ndviClip.normalizedDifference(["band097", "band055"]).rename('NDVI'))
});

// Create a time series chart, with image, geometry & median reducer
var plotNDVI = ui.Chart.image.seriesByRegion(ndvi, geometry, ee.Reducer.median(), 
              'NDVI', 100, 'system:time_start') // band, scale, x-axis property
              .setChartType('LineChart').setOptions({
                title: 'Median NDVI for Selected Geometry',
                hAxis: {title: 'Date'},
                vAxis: {title: 'NDVI'},
                legend: {position: "none"},
                lineWidth: 1,
                pointSize: 3
});

// Display the chart
print(plotNDVI);

We can see how much NDVI has increased in 2021 relative to the earlier years, which makes sense when we look at the reflectance RGB composites - it is much greener in 2021! While this line plot doesn't show us a lot of information now, as the AOP data set builds up in years to come, this may become a more interesting figure.

On your own, we encourage you to dig into the code from this tutorial and modify according to your scientific interests. Think of some questions you have about this dataset, and modify these functions or try writing your own function to answer your question. For example, try out a different reducer, repeat the plots for different areas of the site, and see if there are any other datasets that you could bring in to help you with your analysis. You can also pull in satellite data and see how the NEON data compares. This is just the starting point!

Get Lesson Code

AOP GEE Internannual Change Exploratory Analysis

Footnotes

• To download the metadata documentation without downloading all the data products, you can go through the process of downloading the data product, and when you get to the files, select only the ".pdf" extension. The Algorithm Theoretical Basis Documents (ATBDs), which can be downloaded either from the data product information page or from the data portal, also discuss the uncertainty and important information pertaining to the data.
• The vertical resolution is related to the outgoing pulse width of the lidar system. Optech Gemini has a 10ns outgoing pulse, while the Riegl Q780 and Optech Galaxy Prime sensors have a 3ns outgoing pulse width. At the nominal flying altitude of 1000m AGL, 10ns translates to a range resolution of ~2m, while 3ns corresponds to 2/3m.
• From 2021 onward, all NEON lidar collections have the improved vertical resolution of .67m, as NEON started operating the Optech Galaxy Prime, which replaced one of the Optech Gemini sensors. This has a 3ns outgoing pulse width, matching the Riegl Q780 system.

Google Earth Engine has a number of built in machine learning tools that are designed to work with multi-band raster data. This simplifies more complex analyses like classification (eg. classifying land types or species). In this example, we demonstrate species classification using a random forest machine learning model, using NEON AOP reflectance and ecosystem structure (CHM) data, and TOS (Terrestrial Observation System) woody vegetation data to train the model. For this example, we'll use airshed boundary of the site CLBJ (Lyndon B. Johnson National Grassland in north-central Texas).

Random Forest Machine Learning Workflow

In this tutorial, we will take you through the following steps to classify species at the NEON site CLBJ. Note that these steps closely align with the more general supervised classification steps, described in the Earth Engine Classification Guide .

Workflow Steps:

1. Specify the input data for display and use in analysis
2. Combine the AOP spectrometer reflectance data with the lidar-derived canopy height model (CHM) to create the predictor dataset
3. Create reference (training/test) data for plant species and land cover types based on NEON vegetation structure data
4. Fit a random forest model based on spectral/CHM composite to map the distribution of plant species at CLBJ
5. Evaluate the model accuracy

Load in and Display the AOP and TOS Data

Let's get started. In this first chunk of code, we'll specify the CLBJ location, read in the pre-processed woody vegetation data, as well as the TOS and Airshed boundaries. Details about these boundaries are described on the NEON Flight Box Design webpage.

The plant species data, contained in the CLBJ_veg_2017_filtered feature collection, was derived from the NEON woody plant vegetation structure data product (DP1.10098.001). Land cover classes (grassland, water, shade) were subsequently added to the Feature Collection from visual inspection of the NEON spectrometer reflectance data product (DP3.30006.001).

// Specify CLBJ location
var geo = ee.Geometry.Point([-97.5706464, 33.4045729])
// Load species/land cover samples feature collection to variable (originally a .csv file extracted from  NEON woody plant vegetation structure data product (DP1.10098.001))
var CLBJ_veg = ee.FeatureCollection('projects/neon/AOP_NEON_Plots/CLBJ_veg_2017_filtered')
// Load terrestrial observation system (TOS) boundary to a variable
var CLBJ_TOS = ee.FeatureCollection('projects/neon/AOP_TOS_Boundaries/D11_CLBJ_TOS_dissolve')
//Load tower airshed boundary to a variable
var CLBJ_Airshed = ee.FeatureCollection('projects/neon/Airsheds/CLBJ_90percent_footprint')

Next, let's display the Digital Terrain Model (DTM) and Canopy Height Model (CHM) from the 2017 CLBJ collection, masking out the no-data values (-9999).

// Display DTM for CLBJ. First, filter the DEM image collection by year, DEM type and geographic location
var CLBJ_DTM2017 = ee.ImageCollection('projects/neon/DP3-30024-001_DEM')
  .filterDate('2017-01-01', '2017-12-31')
  .filterMetadata('Type', 'equals', 'DTM')
  .filterBounds(geo)
  .first();
// Then mask out the no-data values (-9999) in the image and add to the map using a histogram stretch based on lower and upper data values
var CLBJ_DTM2017mask = CLBJ_DTM2017.updateMask(CLBJ_DTM2017.gte(0.0000));
Map.addLayer(CLBJ_DTM2017mask, {min:285, max:1294}, 'CLBJ DTM 2017',0);

// Display CHM for CLBJ. First, filter the DEM image collection by year, DEM type and geographic location
var CLBJ_CHM2017 = ee.ImageCollection('projects/neon/DP3-30024-001_DEM')
  .filterDate('2017-01-01', '2017-12-31')
  .filterMetadata('Type', 'equals', 'CHM')
  .filterBounds(geo)
  .first();
  
// Mask out the no-data values (-9999) in the image and add to the map using a histogram stretch based on lower and upper data values
var CLBJ_CHM2017mask = CLBJ_CHM2017.updateMask(CLBJ_CHM2017.gte(0.0000));
Map.addLayer(CLBJ_CHM2017mask, {min:0, max:33}, 'CLBJ CHM 2017',0);

We also want to pull in the Surface Directional Reflectance (SDR) data. When we do this, we want to keep only the valid bands. Water vapor absorbs light between wavelengths 1340-1445 nm and 1790-1955 nm, and the atmospheric correction that converts radiance to reflectance subsequently results in spikes in reflectance in these two band windows. For more information on the water vapor bands, refer to the lesson Plot Spectral Signatures in Python. We will also remove the last 10 bands, as the bands in this region also tend to be noisy.

To remove bands in GEE, you can specify the bands to exclude (here we named this bandsToRemove) and use the .removeAll function to keep only the valid bands. Note we are including as much spectral information as possible for this tutorial, but you could select a smaller subset of bands and likely obtain similar results. We encourage you to test this out on your own. When running the classification on a larger area, it may be a valuable trade-off to include a smaller number of bands so the code runs faster (or doesn't run out of memory).

// Display Surface Directional Reflectance (SDR) image for CLBJ. First, filter the image collection by year, type and geographic location
var CLBJ_SDR2017 = ee.ImageCollection('projects/neon/DP3-30006-001_SDR')
  .filterDate('2017-01-01', '2017-12-31')
  .filterBounds(geo)
  .first()

// Then select all bands except water absorption bands (band195-band205 and band287-band310), as well as the last 10 bands, which also tend to be noisy
var bandNames = CLBJ_SDR2017.bandNames()
var bandsToRemove = ['band195','band196','band197','band198','band199','band200','band201','band202','band203','band204','band205','band287','band288','band289','band290','band291','band292','band293','band294','band295','band296','band297','band298','band299','band300','band301','band302','band303','band304','band305','band306','band307','band308','band309','band310','band416','band417','band418','band419','band420','band421','band422','band423','band424','band425']
var bandsToKeep = bandNames.removeAll(bandsToRemove)
var CLBJ_SDR2017subset = CLBJ_SDR2017.select(bandsToKeep)

print('CLBJ_SDR2017 Valid Band Subset')
print(CLBJ_SDR2017subset)

//Then mask out the no-data values (-9999) in the image and add to the map using a histogram stretch based on lower and upper data values
var CLBJ_SDR2017mask = CLBJ_SDR2017subset.updateMask(CLBJ_SDR2017subset.gte(0.0000)).select(['band053', 'band035', 'band019']);
Map.addLayer(CLBJ_SDR2017mask, {min:.5, max:10}, 'CLBJ SDR 2017', 0);

Next we can combine the SDR bands (381 after we removed the water vapor bands + the last 10 bands) and the CHM band to create a composite multi-band raster that will become our predictor variable (what we use to generate the random forest classification model). We can then crop this to the tower airshed boundary so we can work with a smaller area. This will speed up the process considerably. Optionally, you could classify the full airshed - this code is commented out, if you want to try that instead.

//Combine the SDR (381 bands) and CHM (1 band) for classification
var CLBJ_SDR_CHMcomposite = CLBJ_SDR2017subset.addBands(CLBJ_CHM2017mask).aside(print, 'Composite SDR and CHM image');

//Crop the reflectance image/CHM composite to the NEON tower airshed boundary and add to the map
var CLBJ_SDR2017_airshed = CLBJ_SDR_CHMcomposite.clip(CLBJ_Airshed) // comment if uncommenting next line of code
//var CLBJ_SDR2017_airshed = CLBJ_SDR_CHMcomposite // uncomment to classify entire SDR scene
Map.addLayer(CLBJ_SDR2017_airshed, {bands:['band053', 'band035', 'band019'], min:.5, max:10}, 'CLBJ-Airshed SDR/CHM 2017');

The next chunk of code just displays the TOS and Airshed polygon layers, and prints some details about the NEON vegetation structure data to the Console.

// Display the TOS boundary polygon layer
Map.addLayer(CLBJ_TOS.style({width: 3, color: "blue", fillColor: "#00000000"}),{},"CLBJ TOS", 0)

// Display the Airshed polygon layer
Map.addLayer(CLBJ_Airshed.style({width: 3, color: "white", fillColor: "#00000000"}),{},"CLBJ Airshed", 0)

// Print details about the NEON vegetation structure Feature Collection to the Console
print(CLBJ_veg, 'All woody plant samples')

If you run all the code so far, you should be able to see the following layers in the map. Expand the variables printed to the console to make sure the # of bands are correct, and the "bad" (water vapor / noisy) bands are removed.

Create Training Data Variables

Now that we've added the relevant AOP data, let's start preparing the training data, which we pulled in at the beginning of the script to the variable CLBJ_veg. This next chunk of code pulls out each species into separate variables (by their taxon ID), and adds a layer to the Map for each of these variables.

var CLBJ_QUMA3 = CLBJ_veg.filter(ee.Filter.inList('taxonID', ['QUMA3']))
Map.addLayer(CLBJ_QUMA3, {color: 'orange'}, 'Blackjack oak', 0);

var CLBJ_QUST = CLBJ_veg.filter(ee.Filter.inList('taxonID', ['QUST']))
Map.addLayer(CLBJ_QUST, {color: 'darkgreen'}, 'Post oak', 0);

var CLBJ_ULAL = CLBJ_veg.filter(ee.Filter.inList('taxonID', ['ULAL']))
Map.addLayer(CLBJ_ULAL, {color: 'cyan'}, 'Winged elm', 0);

var CLBJ_ULCR = CLBJ_veg.filter(ee.Filter.inList('taxonID', ['ULCR']))
Map.addLayer(CLBJ_ULCR, {color: 'purple'}, 'Cedar elm', 0);

var CLBJ_GRSS = CLBJ_veg.filter(ee.Filter.inList('taxonID', ['GRSS']))
Map.addLayer(CLBJ_GRSS, {color: 'lightgreen'}, 'Grassland', 0);

var CLBJ_WATR = CLBJ_veg.filter(ee.Filter.inList('taxonID', ['WATR']))
Map.addLayer(CLBJ_WATR, {color: 'blue'}, 'Water', 0);

var CLBJ_SHADE = CLBJ_veg.filter(ee.Filter.inList('taxonID', ['SHADE']));
Map.addLayer(CLBJ_SHADE, {color: 'black'}, 'Shade', 0);

Train/Test Split

Once we have the training data for each species, we can split the data for each species into training and test data, using an 80/20 split. The training data will be used later on to train the random forest model, and the test data is used to test the accuracy of the model results on an independent data set.

// Create training and test subsets for each class (i.e., species types) using stratified random sampling (80/20%)

var new_table = CLBJ_CELA.randomColumn({seed: 1});
var CELAtraining = new_table.filter(ee.Filter.lt('random', 0.80));
var CELAtest = new_table.filter(ee.Filter.gte('random', 0.80));

var new_table = CLBJ_JUVI.randomColumn({seed: 1});
var JUVItraining = new_table.filter(ee.Filter.lt('random', 0.80));
var JUVItest = new_table.filter(ee.Filter.gte('random', 0.80));

var new_table = CLBJ_PRME.randomColumn({seed: 1});
var PRMEtraining = new_table.filter(ee.Filter.lt('random', 0.80));
var PRMEtest = new_table.filter(ee.Filter.gte('random', 0.80));

var new_table = CLBJ_QUMA3.randomColumn({seed: 1});
var QUMA3training = new_table.filter(ee.Filter.lt('random', 0.80));
var QUMA3test = new_table.filter(ee.Filter.gte('random', 0.80));

var new_table = CLBJ_QUST.randomColumn({seed: 1});
var QUSTtraining = new_table.filter(ee.Filter.lt('random', 0.80));
var QUSTtest = new_table.filter(ee.Filter.gte('random', 0.80));

var new_table = CLBJ_ULAL.randomColumn({seed: 1});
var ULALtraining = new_table.filter(ee.Filter.lt('random', 0.80));
var ULALtest = new_table.filter(ee.Filter.gte('random', 0.80));

var new_table = CLBJ_ULCR.randomColumn({seed: 1});
var ULCRtraining = new_table.filter(ee.Filter.lt('random', 0.80));
var ULCRtest = new_table.filter(ee.Filter.gte('random', 0.80));

var new_table = CLBJ_GRSS.randomColumn({seed: 1});
var GRSStraining = new_table.filter(ee.Filter.lt('random', 0.80));
var GRSStest = new_table.filter(ee.Filter.gte('random', 0.80));

var new_table = CLBJ_WATR.randomColumn({seed: 1});
var WATRtraining = new_table.filter(ee.Filter.lt('random', 0.80));
var WATRtest = new_table.filter(ee.Filter.gte('random', 0.80));

var new_table = CLBJ_SHADE.randomColumn({seed: 1});
var SHADEtraining = new_table.filter(ee.Filter.lt('random', 0.80));
var SHADEtest = new_table.filter(ee.Filter.gte('random', 0.80));

Now we can merge all the training data for each species together to create the training data (variable), and similarly merge the test data to create the full test data. From those data we'll create a Features variable containing the predictor data (from the spectral and CHM composite) for the training data.

// Combine species-type reference points for training partition
var training = (CELAtraining).merge(JUVItraining).merge(PRMEtraining).merge(QUMA3training).merge(QUSTtraining).merge(ULALtraining).merge(ULCRtraining).merge(GRSStraining).merge(WATRtraining).merge(SHADEtraining).aside(print, 'Training partition');
var test = (CELAtest).merge(JUVItest).merge(PRMEtest).merge(QUMA3test).merge(QUSTtest).merge(ULALtest).merge(ULCRtest).merge(GRSStest).merge(WATRtest).merge(SHADEtest).aside(print, 'Test partition');
var points = training.merge(test).aside(print,'All points');

// Extract spectral signatures from airshed reflectance/CHM image for training sample based on species ID 
var Features = CLBJ_SDR2017_airshed.sampleRegions({
  collection: training,
  properties: ['taxonIDnum'],
  scale: 1,
  tileScale:16
});
print('Features with spectral signatures:', Features)

Generate and Apply the Random Forest Model

Now that we've assembled the training and test data, and generated predictor data for all of the training data, we can train our random forest model to creat the trainedClassifier variable, and apply it to the reflectance-CHM composite data covering the airshed using .classify(trainedClassifier). Generating the random forest model is pretty simple, once everything is set up!

// Train a 100-tree random forest classifier based on spectral signatures from the training sample for each species 
var trainedClassifier = ee.Classifier.smileRandomForest(100)
.train({
  features: Features,
  classProperty: 'taxonIDnum',
  inputProperties: CLBJ_SDR2017_airshed.bandNames()
});

// Classify the reflectance/CHM image from the trained classifier
var classified = CLBJ_SDR2017_airshed.classify(trainedClassifier);

Assess Model Performance

Next we can assess the performance of this classificiation, by looking at some different metrics. The train accuracy should be high (close to 1, or 100%), as it is testing the performance on the data used to generate the model - however, it is not an accurate representation of the actual accuracy. Instead, the test accuracy tends to be a little more reliable, as it is an independent assessment on the separate (test) data set.

// Calulate a confusion matrix and overall accuracy for the training sample
// Note that this overestimates the accuracy of the model since it does not consider the test sample
var trainAccuracy = trainedClassifier.confusionMatrix().accuracy();
print('Train Accuracy', trainAccuracy);

If you look at the console, you should see a train accuracy of 0.973, pretty good! But we expect this to be good, because we are calculating the accuracy of the samples we used to generate the model. It is not representative of the actual model accuracy.

We can also look at some other accuracy metrics. We won't go into details on each of these, but highlight some of the main takeaways for each of the metrics. For more information on accuracy assessments, you can also refer to Qiusheng Wu's Accuracy Assessment for Image Classification Tutorial

• Confusion Matrix: 2D matrix for a classifier based on its training data, where Axis 0 of the matrix corresponds to the input classes (reference data), and axis 1 corresponds to the output classes (classified data). The rows and columns start at class 0 and increase sequentially up to the maximum class value.
• Overall Accuracy: conveys what proportion were mapped correctly out of all of the reference sites (includes training and test data)
• Kappa Coefficient: evaluates how well a classification performs as compared to randomly assigning classes
• Producer's Accuracy: frequency with which a real feature on the ground is correctly shown in the classified map (corresponds to error of omission)
• User's Accuracy: frequency with which a feature in the classified map will actually be present on the ground (corresponds to error of commission)

// Test the classification accuracy (more reliable estimation of accuracy)
// Extract spectral signatures from airshed reflectance image for test sample
var test = CLBJ_SDR2017_airshed.sampleRegions({
  collection: test,
  properties: ['taxonIDnum'],
  scale: 1,
  tileScale: 16
});

// Calculate different test accuracy estimates
var Test = test.classify(trainedClassifier);
print('Confusion Matrix', Test.errorMatrix('taxonIDnum', 'classification'));
print('Overall Accuracy', Test.errorMatrix('taxonIDnum', 'classification').accuracy());
print('Kappa Coefficient', Test.errorMatrix('taxonIDnum', 'classification').kappa());
print('Producers Accuracy', Test.errorMatrix('taxonIDnum', 'classification').producersAccuracy());
print('Users Accuracy', Test.errorMatrix('taxonIDnum', 'classification').consumersAccuracy());

When you look in the Console and expand the metrics, you can assess the values. Note each taxon ID is assigned a number so you will need to refer to the order in the code to understand the Confusion Matrix as well as the Producer's and User's Accuracy. Each class is listed below for reference, along with the number of training samples for each of the species, in parentheses. When interpreting accuracy, it's important to consider how many training samples were used to generate the model; the model accuracy is impacted if there is poor representation in the training data.

0. CELA - Sugarberry (12)
1. JUVI - Eastern Redcedar (35)
2. PRME - Mexican Plum (2)
3. QUMA3 - Blackjack Oak (26)
4. QUST - Post Oak (162)
5. ULAL - Winged Elm (5)
6. ULCR - Cedar Elm (18)
7. GRSS - Grass (8)
8. WATR - Water (6)
9. SHADE - Shade (9)

Display Model Results

Lastly, we can write a function to display the image classification

// Function used to display training data and image classification
function showTrainingData(){
  var colours = ee.List(["yellow", "white", "green", "orange", "darkgreen", "cyan", "purple","lightgreen","blue","black"]);
  var lc_type = ee.List(["CELA", "JUVI", "PRME","QUMA3","QUST","ULAL","ULCR","GRSS","WATR","SHADE"]);
  var lc_label = ee.List([0, 1, 2, 3, 4, 5, 6, 7, 8, 9]);

  var lc_points = ee.FeatureCollection(
    lc_type.map(function(lc){
      var colour = colours.get(lc_type.indexOf(lc));
      return points.filterMetadata("taxonIDnum", "equals", lc_type.indexOf(lc))
                  .map(function(point){
                    return point.set('style', {color: colour, pointShape: "diamond", pointSize: 3, width: 2, fillColor: "00000000"});
                  });
        })).flatten();

  Map.addLayer(classified, {min: 0, max: 9, palette: ["yellow","white", "green", "orange", "darkgreen", "cyan", "purple", "lightgreen", "blue", "black"]}, 'Classified image', false);
  Map.addLayer(lc_points.style({styleProperty: "style"}), {}, 'All training sample points', false);
  Map.centerObject(geo, 16)
}

// Display the training data and image classification using the showTrainingData function
showTrainingData();

Here is our final classified image:

Acknowledgements

This tutorial was modified from a lesson developed by the Organization for Tropical Studies as part of the course Google Earth Engine for Ecology and Conservation. Thank you OTS!

Get Lesson Code

AOP GEE Random Forest Classification