You can then extract loss areas by selecting a value threshold on the anomaly (Fig. 48.3):

**Fig. 48.3** NDVI anomaly losses (orange) and gains (blue) based on change threshold values

```
var thresholdLoss = -0.05;
var lossfromndvi = anomaly.lte(thresholdLoss)
    .selfMask()
    .updateMask(
        mangrove2000
    ); // Only show the losses within the mangrove extent
of year 2000
Map.addLayer(lossfromndvi, {
    palette: ['orange']
}, 'Loss from Anomaly 00-20');
var thresholdGain = 0.20;
var gainfromndvi = anomaly.gte(thresholdGain)
    .selfMask()
    .updateMask(
        extentBuffer
    ); // Only show the gains within the mangrove extent
buffer of year 2000
Map.addLayer(gainfromndvi, {
    palette: ['blue']
}, 'Gain from Anomaly 00-20');
```
**Code Checkpoint A33b**. The book's repository contains a script that shows what your code should look like at this point.

**Question 4**. What are the main challenges of a map-to-image change detection approach? Think from a perspective of baseline extent, spectral index, and metric used (e.g., mean versus median).

#### **48.3 Synthesis**

With the content covered in this chapter, you will be able to reproduce a map-to-map and a map-to-image change detection to your own area of interest. Additionally, the anomaly analysis can be used to detect changes in other land cover types, including (but not limited to) forests (forest cover loss) and agricultural land (crop harvest).

**Assignment 1**. Practice your land cover classification skills by using the code in script **A33s1—Supplemental** in the book's repository. In this code, we show how to create a mangrove extent map for Guinea using manually and automatically selected samples. The code covers masking techniques and using data from Google Earth Engine Catalog to assist in the classification workflow. Using the techniques you have learned, create two 30 m Landsat-based mangrove extent maps for the year of 2000 and 2020 for Guinea, West Africa.

**Assignment 2**. Using what you learned in this chapter, compile the areas (in km2) of mangrove change in Guinea derived from the anomaly analysis using: (a) vegetation spectral index versus a water-based spectral index; (b) mean versus median; and (c) five-year intervals (2000–2005, 2005–2010, 2010–2015, and 2015–2020.

**Assignment 3**. Practice the classification and compare the map-to-map and mapto-image change approaches for another land cover type/area of your choice, such as the following:

#### **48.4 Conclusion**

Google Earth Engine's computing infrastructure has been revolutionizing timeconsuming remote sensing processes, creating a new way for rapid land cover classification and change detection at large scales. In the case of mangrove forests—a highly dynamic ecosystem that has been under increasing anthropogenic pressure—these approaches allow for a rapid and consistent monitoring of change in extent over time. These approaches using Earth Engine may be particularly useful for developing countries where, until very recently, the high computational power and the difficulties of distributing nontrivial classification algorithms across multiple computational workstations has challenged classification and change detection at large scales.

#### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **49 Forest Degradation and Deforestation**

Carlos Souza Jr. , Karis Tenneson , John Dilger , Crystal Wespestad , and Eric Bullock

#### **Overview**

Tropical forests are being disturbed by deforestation and forest degradation at an unprecedented pace (Hansen et al. 2013; Bullock et al. 2020). Deforestation completely removes the original forest cover and replaces it with another land cover type, such as pasture or agriculture fields. Generally speaking, forest degradation is a temporary or permanent disturbance, often caused by predatory logging, fires, or forest fragmentation, where the tree loss does not entirely change the land cover type. Forest degradation leads to a more complex environment with a mixture of vegetation, soil, tree trunks and branches, and fire ash. Defining a boundary between deforestation and forest degradation is not straightforward; at the time this chapter was

C. Souza Jr. (B)

K. Tenneson · J. Dilger · C. Wespestad Spatial Informatics Group, Pleasanton, CA, USA e-mail<EMAIL_ADDRESS>
J. Dilger e-mail<EMAIL_ADDRESS>
C. Wespestad e-mail<EMAIL_ADDRESS>
E. Bullock Boston University, Boston, MA, USA e-mail<EMAIL_ADDRESS>
K. Tenneson · J. Dilger · C. Wespestad SERVIR-Amazonia, Cali, Colombia

J. Dilger Astraea, Charlottesville, VA, USA

Amazon Institute of People and the Environment, Belém, Brazil e-mail<EMAIL_ADDRESS>
written, there was no universally accepted definition for forest degradation (Aryal et al. 2021). Furthermore, the signal of forest degradation often disappears within one to two years, making degraded forests spectrally similar to undisturbed forests. Due to these factors, detecting and mapping forest degradation with remotely sensed optical data is more challenging than mapping deforestation.

The purpose of this chapter is to present a spectral unmixing algorithm and the normalized difference fraction index (NDFI) to detect and map both forest degradation and deforestation in tropical forests. This spectral unmixing model uses a set of generic endmembers (Souza et al. 2005) to process any Landsat Surface Reflectance (Tier 1) scene available in Google Earth Engine. We present two examples of change detection applications: one comparing a pair of images acquired at different times a year apart by making a temporal color composite and an empirically defined change threshold, and another using a more extensive and dense time series approach.

#### **Learning Outcomes**

#### **Helps if you know how to**

#### **49.1 Introduction to Theory**

Landsat imagery has been extensively used to monitor deforestation (Woodcock et al. 2020). However, detecting and mapping forest degradation associated with selective logging is more intricate and challenging. First efforts involved the application of spectral and textural indices to enhance the detection of canopy damage created by logging (Asner et al. 2002; Souza et al. 2005), but these turned out to be more helpful in enhancing logging infrastructure using Landsat shortwave infrared bands (i.e., roads and log landings; Matricardi et al. 2007).

Alternatively, spectral mixture analysis (SMA) has been proposed to overcome the challenge of using whole-pixel information to detect and classify forest degradation. Landsat pixels typically contain a mixture of land cover components (Adams et al. 1995). The SMA method is based on the linear spectral unmixing model, as described in Chap. 9. The identification of the nature and number of pure spectra (often referred to in this context as "endmembers") in the image scene is an important step in obtaining correct SMA models. In logged forests (and also in burned forest and forest edges), mixed pixels predominate and are expected to have a combination of green vegetation (GV), soil, non-photosynthetic vegetation (NPV), and shade-covered materials. Therefore, fraction images derived from SMA analyses are more suitable to enhance the detectability of logging infrastructure and canopy damage within degraded forests. For example, soil fractions reveal log landings and logging roads (Souza and Barreto 2000), while the NPV highlights forest canopy damage (Souza et al. 2003), and the areas decreasing in GV indicate forest canopy gaps (Asner et al. 2004).

A study has shown that it is possible to generalize the SMA model to Landsat sensors (including TM, ETM+, and OLI) (Small 2004). Souza et al. (2005) expanded the generalized SMA to handle five endmembers—GV, NPV, soil, shade, and cloud—expected within forest degradation areas and proposed a novel compositional index based on SMA fractions, the normalized difference fraction index (NDFI).

The NDFI is computed as:

$$\text{NDFI} = \frac{\text{GVshare} - (\text{NPV} + \text{Soil})}{\text{GVshare} + \text{NPV} + \text{Soil}} \tag{49.1}$$

where GVshade is the shade-normalized GV fraction given by,

$$\text{GV share} = \frac{\text{GV}}{100 - \text{Shade}} \tag{49.2}$$

NDFI values range from − 1 to 1. For intact forests, NDFI shows high values (i.e., about 1) due to the combination of high GVshade (i.e., high GV and canopy shade) and low NPV and soil values. The NPV and soil fractions increase as forests are more degraded, lowering NDFI values relative to the intact forests. Deforested areas exhibit very low GV and shade and high NPV and soil, making it possible to distinguish them from degraded forests based on NDFI magnitude.

Recent studies compared NDFI with other spectral indices. NDFI generated more accurate results in deforestation detection (Schultz et al. 2016) and forest degradation (Bullock et al. 2018) in time-series analysis. One of the key components for its success is lowering unwanted noise and accounting for illumination variability through the shade normalization applied to the GV fraction.

#### **49.2 Practicum**

#### **49.2.1 Section 1: Spectral Mixture Analysis Model**

Let us first define the Landsat endmembers based on Souza et al. (2005). These endmembers were developed and tested in the Amazon. These endmembers work well in many other environments (see example applications for calculating NDFI in non-Amazonian tropical forest in Schultz et al. 2016 [Ethiopia and Viet Nam]; Kusbach et al. 2017 [central Europe]; Hirschmugl et al. 2013 [Cameroon and Central African Republic]). If you are working in a different region, assess how well they perform for your forest types. The ratio of these endmembers that makes up the spectral signature of each pixel gives a good indication of the plant health and composition for that area. When the ratio shifts over time towards one or more of the endmembers, we can quantify how the landscape is changing.

Below, we will create a new variable endmembers by copying the values from the code block below. The six numbers in square brackets define the pure reflectance values for the blue, green, red, SWIR1, and SWIR2 bands for each endmember material.

```
// SMA Model - Section 1
// Define the Landsat endmembers (source: Souza et al.
2005)
// They can be applied to Landsat 5, 7, 8, and potentially
9.
var endmembers = [
  [0.0119,0.0475,0.0169,0.625,0.2399,0.0675], // GV
  [0.1514,0.1597,0.1421,0.3053,0.7707,0.1975], // NPV
  [0.1799,0.2479,0.3158,0.5437,0.7707,0.6646], // Soil
  [0.4031,0.8714,0.79,0.8989,0.7002,0.6607] // Cloud
];
```
We will choose a single Landsat 5 image to work with for now, and select the bands we will need for the SMA and NDFI calculation.

Create a new variable image and assign it to the Landsat 5 image from the code block below. Select the visible, near infrared, and shortwave infrared bands. Then, center the map around the Landsat 5 image with a zoom scale of 10.

```
// Select a Landsat 5 scene on which to apply the SMA
model.
var image =
ee.Image('LANDSAT/LT05/C02/T1_L2/LT05_226068_19840411')
    .multiply(0.0000275).add(-0.2);
// Center the map on the image object.
Map.centerObject(image, 10);
// Define and select the Landsat bands to apply the SMA
model.
// use ['SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B5',
'SR_B7'] for Landsat 5 and 7.
// use ['SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6',
'SR_B7'] for Landsat 8.
var bands = ['SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B5',
'SR_B7'];
image = image.select(bands);
```
Next, we will need to create a couple of functions to use the endmembers and create the NDFI image. We will need to unmix the input Landsat image.

First, we will create a new function, getSMAFractions that takes two parameters: image and endmembers. To unmix the image, we first select the visible, near-infrared, and short wave infrared bands, and call the unmix function with the endmembers used as the argument. Note: The order in which the bands are selected matters. For each endmember (GV, NPV, soil, and cloud), each value in the array represents the pure value for each band being passed in. For example, if the selected bands were ['SR\_B2', 'SR\_B3', 'SR\_B4', 'SR\_B5', 'SR\_B6', 'SR\_B7'] then the selected endmembers at position 0 would be:

```
// Unmixing image using Singular Value Decomposition.
var getSMAFractions = function(image, endmembers) {
    var unmixed = ee.Image(image)
        .select([0, 1, 2, 3, 4,
            5
        ]) // Use the visible, NIR, and SWIR bands only!
        .unmix(endmembers)
        .max(0) // Remove negative fractions, mostly Soil.
        .rename('GV', 'NPV', 'Soil', 'Cloud');
    return ee.Image(unmixed.copyProperties(image));
};
```
We will now write the function to calculate NDFI. We could write it out line by line for each time we want to perform SMA or calculate NDFI. But having the equation implemented as a function gives us many benefits, including reducing redundancy, allowing us to map the function over a collection of images, and reducing errors from typos or other little bugs that can find their way into our code.

We will use the fraction images obtained with the getSMAFractions function above to calculate shade, GVs, and NDFI using image expressions. This procedure will return a multiband image with the shade, GVs, and NDFI bands added to the input image.

First, we will create a new variable sma and pass in the image and endmembers as arguments. Then calculate the shade and GV shade-normalized (GVs) fractions from the SMA bands, and add the shade and GVs bands to the SMA image. We calculate NDFI using an expression implementing Eq. 49.1, and add the new band to the SMA image.

```
// Calculate GVS and NDFI and add them to image fractions.
// Run the SMA model passing the Landsat image and the
endmembers.
var sma = getSMAFractions(image, endmembers);
Map.addLayer(sma, {
    bands: ['NPV', 'GV', 'Soil'],
    min: 0,
    max: 0.45
}, 'sma');
```

```
// Calculate the Shade and GV shade-normalized (GVs)
fractions from the SMA bands.
var Shade = sma.reduce(ee.Reducer.sum())
    .subtract(1.0)
    .abs()
    .rename('Shade');
var GVs = sma.select('GV')
    .divide(Shade.subtract(1.0).abs())
    .rename('GVs');
// Add the new bands to the SMA image variable.
sma = sma.addBands([Shade, GVs]);
// Calculate the NDFI using image expression.
var NDFI = sma.expression(
    '(GVs - (NPV + Soil)) / (GVs + NPV + Soil)', {
        'GVs': sma.select('GVs'),
        'NPV': sma.select('NPV'),
        'Soil': sma.select('Soil')
    }).rename('NDFI');
// Add the NDFI band to the SMA image.
sma = sma.addBands(NDFI);
```
We will use a color palette that spans from white, to pink (i.e., bare land), to yellow, to green to visualize the NDFI image. Higher values of NDFI will be green, while lower values will span the colors of white, pink, and yellow. Copy the code block below into your Code Editor.

```
// Define NDFI color table.
var palettes = require(
    'projects/gee-edu/book:Part A - Applications/A3 -
Terrestrial Applications/A3.4 Forest Degradation and
Deforestation/modules/palettes'
);
var ndfiColors = palettes.ndfiColors;
```
Next, we can visualize all the hard work we have done unmixing each endmember and the NDFI bands (Fig. 49.1).

Create an image visualization object with bands 5, 4, and 3, chosen for visualization along with a min and max. Add the Landsat 5 image to the map using the image visualization object.

**Fig. 49.1** Example maps of **a** green vegetation shade normalized fraction (more vegetation is whiter); **b** shade fraction (more shade is whiter); **c** non-photosynthetic vegetation fraction (more NPV is whiter); **d** green vegetation fraction (more vegetation is whiter); **e** soil fraction (more soil is whiter)

Now add each of the SMA bands and the NDFI band to the map. Note: Rather than defining the min, max, and bands for each of these in separate variables, we can pass in the object directly to the Map.addLayer function.

```
var imageVis = {
    'bands': ['SR_B5', 'SR_B4', 'SR_B3'],
    'min': 0,
    'max': 0.4
};
// Add the Landsat color composite to the map.
Map.addLayer(image, imageVis, 'Landsat 5 RGB-543', true);
// Add the fraction images to the map.
Map.addLayer(sma.select('Soil'), {
    min: 0,
    max: 0.2
}, 'Soil');
Map.addLayer(sma.select('GV'), {
    min: 0,
    max: 0.6
}, 'GV');
Map.addLayer(sma.select('NPV'), {
    min: 0,
    max: 0.2
}, 'NPV');
Map.addLayer(sma.select('Shade'), {
    min: 0,
    max: 0.8
}, 'Shade');
Map.addLayer(sma.select('GVs'), {
    min: 0,
    max: 0.9
}, 'GVs');
Map.addLayer(sma.select('NDFI'), {
    palette: ndfiColors
}, 'NDFI');
```
The last thing we will do in this section is to create a water and cloud mask. Water and clouds will have lower NDFI values. While water may not be too much of an issue, clouds will impact how we monitor forest degradation and loss, and thus we can simply mask them out. We can mask them using a thresholding method based on the values of our fraction images.

First, create a new function variable getWaterMask that takes an SMA image as the only argument.

Next, create a water mask using threshold values for the shade, GV, and soil bands, where shade is greater than or equal to 0.65; GV is less than or equal to 0.15; and soil is less than or equal to 0.05.

Now create a cloud mask by applying a threshold of 0.1 or greater to the Cloud band.

```
var getWaterMask = function(sma) {
    var waterMask = (sma.select('Shade').gte(0.65))
        .and(sma.select('GV').lte(0.15))
        .and(sma.select('Soil').lte(0.05));
    return waterMask.rename('Water');
};
// You can use the variable below to get the cloud mask.
var cloud = sma.select('Cloud').gte(0.1);
var water = getWaterMask(sma);
```
Next, we will combine the cloud and water masks using the max reducer. Since we want to mask both water and clouds, using the max works quite nicely here as opposed to adding the images together—so we do not need to worry about overlaps.

Now add the cloud and water mask as a layer to the map.

Apply the cloud and water mask to the NDFI band using the updateMask function and invert the mask using the not function. Note: updateMask considers zeroes as invalid (i.e., to be masked) and ones as valid (i.e., to be kept). Since our original mask had values of 1 for cloud and water, we use the not function to invert the values.

```
var cloudWaterMask = cloud.max(water);
Map.addLayer(cloudWaterMask.selfMask(),
    {
        min: 1,
        max: 1,
        palette: 'blue'
    },
    'Cloud and water mask');
// Mask NDFI.
var maskedNDFI =
sma.select('NDFI').updateMask(cloudWaterMask.not());
Map.addLayer(maskedNDFI, {
    palette: ndfiColors
}, 'NDFI');
```
**Code Checkpoint A34a**. The book's repository contains a script that shows what your code should look like at this point.

#### **49.2.2 Section 2: Deforestation and Forest Degradation Change Detection**

To observe changes in the landscape over time, you need to create an NDFI image for two points in time and then calculate the difference between them. Previous studies have shown that images should not be more than one year apart because the forest degradation disturbance quickly disappears with tree foliage and understory vegetation growth. Changes in NDFI are a good indicator of forest change.

Now we will start using the functions we have built. First, we perform the SMA on a Landsat 5 image. Recall that the SMA function wants only the visible, nearinfrared, and shortwave infrared bands. When using Landsat 5, those correspond to bands 1–5 and band 7.

First, create a variable for the Landsat 5 scene specified in the code block below and select the visible, near-infrared, and shortwave infrared bands. Use the getSMAFractions function with the Landsat image and the endmembers. Rename the output SMA bands as GV, NPV, Soil, and Cloud.

```
// Select two Landsat 5 scenes on which to apply the SMA
model.
// Select Landsat bands used for forest change detection.
var imageTime0 = ee.Image(
        'LANDSAT/LT05/C02/T1_L2/LT05_226068_20000509')
    .multiply(0.0000275).add(-0.2);
var bands = ['SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B5',
'SR_B7'];
imageTime0 = imageTime0.select(bands);
// Run the SMA model.
var smaTime0 = getSMAFractions(imageTime0, endmembers);
```
Before we move on to calculate the NDFI, we will add the previous Landsat scene and the fractional images to the map to inspect them (example results in Fig. 49.2).

**Fig. 49.2** Example maps of **a** RGB composite of bands 5–4–3 (green is vegetation and brown is barren ground); **b** soil fraction (more soil is whiter); **c** green vegetation fraction (more vegetation is whiter); **d** non-photosynthetic vegetation fraction (more NPV is whiter)

```
// Center the image object.
Map.centerObject(imageTime0, 10);
// Define the visualization parameters.
var imageVis = {
    'opacity': 1,
    'bands': ['SR_B5', 'SR_B4', 'SR_B3'],
    'min': 0,
    'max': 0.4,
    'gamma': 1
};
// Scale to the expected maximum fraction values.
var fractionVis = {
    'opacity': 1,
    'min': 0.0,
    'max': 0.5
};
// Add the Landsat color composite to the map.
Map.addLayer(imageTime0, imageVis, 'Landsat 5 RGB 543',
true);
// Add the fraction images to the map.
Map.addLayer(smaTime0.select('Soil'), fractionVis, 'Soil
Fraction');
Map.addLayer(smaTime0.select('GV'), fractionVis, 'GV
Fraction');
Map.addLayer(smaTime0.select('NPV'), fractionVis, 'NPV
Fraction');
```
Next, let's set up a function to systematically reproduce the work we did computing NDFI, since we will need it a couple more times. Hereafter, we will be able to call that function instead of needing to explicitly write out each step, thus simplifying our code and making it easier to change if we need to later.

```
function getNDFI(smaImage) {
    // Calculate the Shade and GV shade-normalized (GVs)
fractions
    // from the SMA bands.
    var Shade = smaImage.reduce(ee.Reducer.sum())
        .subtract(1.0)
        .abs()
        .rename('Shade');
    var GVs = smaImage.select('GV')
        .divide(Shade.subtract(1.0).abs())
        .rename('GVs');
    // Add the new bands to the SMA image variable.
    smaImage = smaImage.addBands([Shade, GVs]);
    var ndfi = smaImage.expression(
        '(GVs - (NPV + Soil)) / (GVs + NPV + Soil)', {
            'GVs': smaImage.select('GVs'),
            'NPV': smaImage.select('NPV'),
            'Soil': smaImage.select('Soil')
        }
    ).rename('NDFI');
    return ndfi;
}
```
Then, calculate NDFI for the earlier Landsat image's SMA bands (use smaTime0 as an input) using the getNDFI function you wrote. Add this NDFI image to the map displayed in the ndfiColors palette (Fig. 49.3). This will serve as your calculated NDFI for your earlier point in time, the pre-change time (this time was defined as imageTime0).

```
// Create the initial NDFI image and add it to the map.
var ndfiTime0 = getNDFI(smaTime0);
Map.addLayer(ndfiTime0,
    {
        bands: ['NDFI'],
        min: -1,
        max: 1,
        palette: ndfiColors
    },
    'NDFI t0',
    false);
```
**Fig. 49.3** NDFI image for the previous year obtained from the Landsat 5 (path/row 226/068) scene acquired on May 9, 2000. Orange colors indicate signs of forest disturbance associated with fires and selective logging. Pink and white colors are dry vegetation and bare soil in old deforested areas. Orange colors in pasturelands mean dry vegetation

Next, you will repeat this procedure to calculate the SMA fractions and NDFI of the second Landsat 5 image (smaTime1). You will utilize the same SMA method.

First, create a new variable (imageTime1), using the Landsat 5 scene from the code block below. Note that the band numbers for Landsat 8 are different from those for Landsat 5. In general, you should always make sure to check band names when working with multiple Landsat collections. Then, you will calculate the SMA fractions. The getSMAFractions function will rename the outputs to "GV", "NPV", "Soil", and "Cloud". Then, you will calculate NDFI for the new scene. Add an RGB composite and the NDFI to the map.

```
// Select a second Landsat 5 scene on which to apply the
SMA model.
var imageTime1 = ee.Image(
        'LANDSAT/LT05/C02/T1_L2/LT05_226068_20010629')
    .multiply(0.0000275).add(-0.2)
    .select(['SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B5',
'SR_B7']);
var smaTime1 = getSMAFractions(imageTime1, endmembers);
// Create the second NDFI image and add it to the map.
var ndfiTime1 = getNDFI(smaTime1);
Map.addLayer(imageTime1, imageVis, 'Landsat 5 t1 RGB-5',
true);
Map.addLayer(ndfiTime1,
    {
        bands: ['NDFI'],
        min: -1,
        max: 1,
        palette: ndfiColors
    },
    'NDFI_t1',
    false);
```
Before looking at changes between the two NDFI images, use the opacity slider in the **Layers** panel of the map to manually view the differences between the NDFI outputs, and then the RGB outputs. Where do you expect to see the greatest changes between the two?

An easier way to enhance changes using a pair of images is to build a temporal color composite. To construct a temporal color composite, we add the first image date to the R color channel and the second to the G and B channels. An example of RGB color composite (as first described in Chap. 2) is shown in Fig. 49.4, using the NDFI\_t0 (R) and NDFI\_t1 (G and B). Deforestation appears in bright red colors since we assigned the first NDFI image to the R color channel, indicating forest in the previous year, and removal in the following (i.e., G and B colors have low NDFI values due to forest removal). In contrast, the cyan colors in the NDFI temporal color composite indicate vegetation regrowth in the second year. The gradient of dark to gray colors suggests no change in NDFI between the two dates.

**Fig. 49.4** Example of NDFI obtained for the first (**a**) and second (**b**) Landsat images, and a visualization of them using a temporal RGB color composite (**c**). Intense red colors are newly deforested areas, and light red are selectively logged forests. Cyan colors indicate vegetation regrowth in areas that were logged or burned a year ago or more. Example of NDFI difference histogram (**d**), with positive values indicating an increase in NDFI over time and negative values a decrease. Values at zero suggest no change

**Fig. 49.4** (continued)

To find the change between our images, a simple difference method can be applied by subtracting the previous image from the current image. Then, we can apply an empirically defined threshold to classify the changes based on the inspection of the histogram and the NDFI temporal color composite. For more information on two-date differencing, see Chap. 16.

Next, we will combine the two NDFI images for time t0 and t1. Create a new variable ndfiChange and subtract the current NDFI image (made using Landsat 5 imagery from 2001) from the previous NDFI image (made using Landsat 5 imagery from 2000). Note: The NDFI is calculated using the image expression method applied to the SMA fraction presented in the code above. We then combine the two NDFI bands from 2000 and 2001 in one variable to display the change over time.

```
// Combine the two NDFI images in a single variable.
var ndfi = ndfiTime0.select('NDFI')
    .addBands(ndfiTime1.select('NDFI'))
    .rename('NDFI_t0', 'NDFI_t1');
// Calculate the NDFI change.
var ndfiChange = ndfi.select('NDFI_t1')
    .subtract(ndfi.select('NDFI_t0'))
    .rename('NDFI Change');
```
Using the polygon drawing tool (as described in Chap. 6), draw a region that covers most of the ndfiChange image and rename it region in the variable import panel.

Next, make a histogram named histNDFIChange that bins the ndfiChange image using the region you just drew, and prints it to the **Console**. Optionally, click **Run** to view the histogram and identify some potential change thresholds. In your case the histogram may look slightly different.

```
var options = {
    title: 'NDFI Difference Histogram',
    fontSize: 20,
    hAxis: {
        title: 'Change'
    },
    vAxis: {
        title: 'Frequency'
    },
    series: {
        0: {
            color: 'green'
        }
    }
};
// Inspect the histogram of the NDFI change image to define
threshold
// values for classification. Make the histogram, set the
options.
var histNDFIChange = ui.Chart.image.histogram(
        ndfiChange.select('NDFI Change'), region, 30)
    .setSeriesNames(['NDFI Change'])
    .setOptions(options);
print(histNDFIChange);
```
Classify the difference image into new deforestation (red), forest degradation (orange), regrowth (cyan), and forest (green). The classification is based on slicing the NDFI difference image. Old deforested areas are detected using the NDFI first image date (t0).

Add the change classification and difference images to the map and click **Run**.

```
// Classify the NDFI difference image based on thresholds
// obtained from its histogram.
var changeClassification = ndfiChange.expression(
        '(b(0) >= -0.095 && b(0) <= 0.095) ? 1 :' +
        // No forest change
        '(b(0) >= -0.250 && b(0) <= -0.095) ? 2 :' + //
Logging
        '(b(0) <= -0.250) ? 3 :' + // Deforestation
        '(b(0) >= 0.095) ? 4 : 0') // Vegetation regrowth
    .updateMask(ndfi.select('NDFI_t0').gt(
        0.60)); // mask out no forest
// Use a simple threshold to get forest in the first image
date.
var forest = ndfi.select('NDFI_t0').gt(0.60);
```
Finally, add code to add all the new layers to the map and click **Run**.

```
// Add layers to map
Map.addLayer(ndfi, {
    'bands': ['NDFI_t0', 'NDFI_t1', 'NDFI_t1']
}, 'NDFI Change');
Map.addLayer(ndfiChange, {}, 'NDFI Difference');
Map.addLayer(forest, {}, 'Forest t0 ');
Map.addLayer(changeClassification, {
        palette: ['000000', '1eaf0c', 'ffc239', 'ff422f',
            '74fff9']
    },
    'Change Classification');
```
Figure 49.5 shows an example of forest changes between two dates displayed using the cutoffs defined by the histogram. Categorizing the whole map into a simple system of no change, old deforestation, new deforestation, partial forest disturbance, and regrowth makes the forest use patterns in the region clearer. However, you must be careful with the thresholds you chose from the histogram when creating and interpreting the NDFI changes in this way, as those values can drastically alter your map. Consider when it might be more appropriate to use an RGB temporal color composite, and in what situations the classified map using the empirically defined thresholds would be better.

**Fig. 49.5** Example of the change detection using two NDFI images. Green is remaining forest; black is old deforestation as detected in the first NDFI image; red highlights new deforestation; orange shows forest disturbance by logging or fires; and cyan is vegetation regrowth of forest since the previous date

**Code Checkpoint A34b**. The book's repository contains a script that shows what your code should look like at this point.

Save your script for your own future use, as outlined in Chap. 1. Then, refresh the page to begin with a new script for the next section.

#### **49.2.3 Section 3: Deforestation and Forest Degradation Time Series Analysis**

To assess deforestation and forest degradation with a time series, we can use a Google Earth Engine tool called CODED (Bullock et al. 2018, 2020; Bullock (2018). The algorithm is based on previous work in continuous land cover monitoring (Zhu and Woodcock 2014) and NDFI-based degradation mapping using spectral unmixing models (Souza et al. 2003, 2005). CODED has both a user interface application and an API, which can be accessed in the Earth Engine Code Editor. For this lesson, we will use the API.

**Code Checkpoint A34c**. The book's repository contains a script to use to begin this section. You will need to start with that script and script and paste code below into it. The checkpoint accesses the modules of the CODED API and support functions. Note: importing large modules can cause your browser to hang for a moment while they load.

Detecting degraded forest regions requires knowledge of the characteristics of the forest of interest when it is in its normal healthy state. Higher NDFI values, near 1, typically indicate a healthy forest, but the magnitude range of NDFI for a healthy forest is dependent on forest density, type of forest ecosystem, and the seasonality of the area. In order to determine the typical values of NDFI observed in each forest over time, CODED employs a training period. Within this period, a regression model, similar to the one initially introduced in Chap. 18, is fitted for NDFI values for each pixel. The regression model is composed of a constant for overall NDFI magnitude, a sine and cosine term encapsulating seasonal and intra-annual variability, and an RMSE to account for noise. In this way, CODED can find typical temporal patterns in the landscape, account for clouds and sensor noise, and better distinguish forests from non-forested areas.

The code below sets up the study area, accesses the Landsat collection to be used, and defines the study period. For this example, we will use the geometry of the previous Landsat scene as our study area. Then, we define a new variable studyArea and assign it to the image we used when first exploring NDFI, retrieving its geometry using the geometry method. Then, we use the utils module to call the Inputs.getLandsat function, which accesses the ImageCollection. We then filter the images to a start date and end date of interest. Note: CODED calibrates to find the typical variations in NDFI for the region by observing the NDFI patterns over time, so it is a good idea to filter the Landsat data so you have an extra six months to a year of data before the time period in which you are truly interested. For example, if you wanted to study forest change from 2000 to 2010, it would be good practice to use 1999 as the start date and 2010 as the end date. Paste the code below into your starter script you opened to begin this section.

```
// We will use the geometry of the image from the previous
section as
// the study area.
var studyArea = ee.Image(
        'LANDSAT/LT05/C02/T1_L2/LT05_226068_19840411')
    .geometry();
// Get cloud masked (Fmask) Landsat imagery.
var landsat = utils.Inputs.getLandsat()
    .filterBounds(studyArea)
    .filterDate('1984-01-01', '2021-01-01');
```
Getting our ImageCollection this way saves us quite a bit of coding. By default the returned collection will use every available Landsat mission, perform some simple cloud masking, and generate our image fractions of green vegetation (GV), soil, non-photosynthetic vegetation (NPV), shade-covered material, and NDFI, as well as additional indices.

First, make a new variable gfwImage and add the path to the Global Forest Change product from Hansen et al. (2013), which is used to create a forest mask. Define a threshold of 40 for the percent canopy cover for the mask. Apply the threshold to the treecover2000 band from the gfwImage. You could also choose to use a pre-prepared forest mask of your own instead of selecting one from the Global Forest Change product. This might be useful if you have a local highquality forest mask, as it may slightly improve your results or provide consistency with other work you have done in the area.

```
// Make a forest mask
var gfwImage =
ee.Image('UMD/hansen/global_forest_change_2019_v1_7');
// Get areas of forest cover above the threshold
var treeCover = 40;
var forestMask = gfwImage.select('treecover2000')
    .gte(treeCover)
    .rename('landcover');
```
To identify degradation and deforestation, distinguished by whether the land cover remains forest or not after the event, we will use a prepared dataset of forest and non-forested areas. This data already has the predictor data for each feature, which will speed up the computation. Alternatively, refer back to Chap. 6 on classification for details on how to create your own unique forest and nonforested area dataset and ensure each feature has a year property with the year collected as an integer.

```
var samples = ee.FeatureCollection(
    'projects/gee-book/assets/A3-
4/sample_with_pred_hansen_2010');
```
There are many parameters that can be adjusted when running CODED, as summarized below.

• Note: prepTraining tells the algorithm to add the coefficients to your samples for training the classifier. This will initiate a task to export the prepared samples for later use if you wish. In this example, we will set it to false.

```
var minObservations = 4;
var chiSquareProbability = 0.97;
var training = samples;
var forestValue = 1;
var startYear = 1990;
var endYear = 2020;
var classBands = ['NDFI', 'GV', 'Shade', 'NPV', 'Soil'];
var prepTraining = false;
```
With the parameters defined, we now have everything needed to run CODED. CODED takes a single argument, which is stored as a JavaScript dictionary. Since we defined all the parameters as variables, it may seem redundant to put them into a dictionary, but having them as variables can be helpful when you are exploring functionality and adjusting parameter values frequently. Of course, entering the values directly into the dictionary would work as well.

Create a new dictionary and assign each parameter variable to a key of the same name.

```
//---------------- CODED parameters
var codedParams = {
    minObservations: minObservations,
    chiSquareProbability: chiSquareProbability,
    training: training,
    studyArea: studyArea,
    forestValue: forestValue,
    forestMask: forestMask,
    classBands: classBands,
    collection: landsat,
    startYear: startYear,
    endYear: endYear,
    prepTraining: prepTraining
};
// -------------- Run CODED
var results = api.ChangeDetection.coded(codedParams);
print(results);
```
Run CODED by clicking **Run**. This will set up the run, and issue the call to the ChangeDetection.coded function (Fig. 49.6).

**Fig. 49.6** The result of running the CODED change detection algorithm is an object with the general and change parameters used for the run and a layers object that has all the image outputs

**Code Checkpoint A34d**. The book's repository contains a script that shows what your code should look like at this point.

Next, you will relabel some of the results so that they are easier to understand and work with. You will rename the degradation layers to something more human readable, like 'degradation\_1', 'degradation\_2', etc. Rename the deforestation layers to something more human readable, like 'deforestation\_1', 'deforestation\_2', etc.

Set a variable for the mask layer. This is the same mask that was passed into CODED, so retrieving it in this method is not strictly necessary.

Set a variable for the change output that is the concatenation of the degradation and deforestation outputs. This is mostly for organizational purposes. Since change is more rare, self-masking removes all the non-change pixels, and casting to Int32 helps keep all the bands in the same type, which you would need for exporting to a geoTIFF.

Set a variable mag to the minimum magnitude. There are multiple magnitude bands that correspond to the number of temporal segments that are retrieved when running CODED. These bands are reduced by the minimum since we want to find areas where the greatest negative change has occurred.

```
// Format the results for exporting.
var degradation = results.Layers.DatesOfDegradation
    .rename(['degradation_1', 'degradation_2',
        'degradation_3', 'degradation_4'
    ]);
var deforestation = results.Layers.DatesOfDeforestation
    .rename(['deforestation_1', 'deforestation_2',
        'deforestation_3', 'deforestation_4'
    ]);
var mask = results.Layers.mask.rename('mask');
var change = ee.Image.cat([degradation,
deforestation]).selfMask()
    .toInt32();
var mag = results.Layers.magnitude.reduce(ee.Reducer.min())
    .rename('magnitude');
```
Finally, we can combine all of this information into a stratified—or classified output layer of forest, non-forest, degradation, and deforestation. We can define a function to take in each of the outputs and apply some logic to decide these categories. A new threshold we need to apply is the magnitude, magThreshold. This threshold will define the minimum amount of change we want to qualify as a degradation or deforestation event. This is a post-processing step in which we will stratify the results into categories of stable forest, stable non-forest, degradation, and deforestation.

Next, you will create a new function named makeStrata that takes an image and a threshold as the arguments. The base of the stratified image will be the mask which is remapped from [0, 1] to [2, 1]. This keeps the forest class at a value of 1 and updates the non-forest class to a value of 2.

Then, you will create a binary mask of the minimum threshold using the magnitude threshold parameter. Then, create a binary degradation image using all the degradation bands, and then multiply it by the magnitude mask. Similarly, you'll create a binary deforestation image using all the deforestation bands and then multiply it by the magnitude mask. Update the strata image using the where functions to first assign degradation to 3 and next to assign deforestation to 4. The where functions are applied sequentially. Deforestation needs to be updated last, because in this case, areas of degradation could overlap with areas that are also deforestation.

```
var makeStrata = function(img, magThreshold) {
    var strata = img.select('mask').remap([0, 1], [2, 1]);
    var mag = img.select('magnitude').lte(magThreshold);
    var deg =
img.select(['deg.*']).gt(0).reduce(ee.Reducer.max())
        .multiply(mag);
    var def =
img.select(['def.*']).gt(0).reduce(ee.Reducer.max())
        .multiply(mag);
    strata = strata.where(deg, 3).where(def, 4);
    return strata.clip(studyArea);
};
```
Concatenate the mask, change, and mag bands into a single image. Define a magnitude threshold of − 0.6. Apply the makeStrata function using the full output image and the magnitude threshold. Then, add code to export the strata to your assets.

```
var fullOutput = ee.Image.cat([mask, change, mag]);
var magnitudeThresh = -0.6;
var strata = makeStrata(ee.Image(fullOutput),
magnitudeThresh)
    .rename('strata');
Export.image.toAsset({
    image: strata,
    description: 'strata',
    region: studyArea,
    scale: 30,
    maxPixels: 1e13,
});
```
Click **Run**. The CODED algorithm is computationally heavy, so the results need to be exported before they can be viewed on the map in the Code Editor. Once the export has finished, you can add the layer to the map. We have created the asset for you and stored it in the book repository; you can access it with the code below:

```
var exportedStrata = ee.Image('projects/gee-book/assets/A3-
4/strata');
Map.addLayer(exportedStrata,
    {
        min: 1,
        max: 4,
        palette: 'green,black,yellow,red'
    },
    'strata');
Map.setCenter(-55.0828, -11.24, 11);
```
The resulting map should look something like the example in Fig. 49.7. With these stratified results, you can see the relative amounts and geographic distribution of forest that have been degraded or deforested in your time period of interest. What patterns do you observe? Does forest degradation or deforestation seem more prevalent?

**Code Checkpoint A34e**. The book's repository contains a script that shows what your code should look like at this point.

**Fig. 49.7** Example of classified CODED results for 1990–2020. Non-forested areas are black, stable forest areas are green, deforestation is red, and degradation is yellow

#### **49.3 Synthesis**

**Assignment 1**. Study how adjusting each parameter affects your results.

**Assignment 2**. Try adjusting the CODED parameters to see how they affect the detection of degradation and deforestation in a study area where you highly suspect from the imagery that these events have occurred. Observe in the map whether you are overestimating or underestimating the disturbed forest area with each parameter combination.

#### **49.4 Conclusion**

In this chapter, you learned about the spectral unmixing algorithm (SMA) and the Normalized Difference Fraction Index (NDFI) in order to map forest degradation and deforestation. We presented two examples of change detection applications: two-image differencing and an NDFI time-series approach using the CODED algorithm. NDFI is sensitive to subtle changes in forest composition, making it ideal for detection of forest degradation. CODED's use of a regression model fitting of NDFI over a training period informs us more about what the NDFI magnitudes and variability of a forest should be in a healthy state, so changes from the norm can be more confidently labeled as forest disturbance events. However, your chosen CODED parameters have an important impact on your resulting map of forest loss and forest degradation.

#### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **50 Deforestation Viewed from Multiple Sensors**

Xiaojing Tang

#### **Overview**

Combining data from multiple sensors is the best way to increase data density and hence detect change faster. The purpose of this chapter is to demonstrate a simple method of combining Landsat, Sentinel-2, and Sentinel-1 data for monitoring tropical forest disturbance. You will learn how to import, preprocess, and fuse optical and synthetic aperture radar (SAR) remote sensing data. You will also learn how to monitor change using time-series models.

### **Learning Outcomes**

### **Helps if you know how to**

X. Tang (B)

James Madison University, Harrisonburg, VA, USA e-mail<EMAIL_ADDRESS>

#### **50.1 Introduction to Theory**

Deforestation and forest degradation are large sources of carbon emissions and negatively impact biodiversity, food security, and human well-being. The ability to quickly and accurately detect forest disturbance events is essential for preventing future forest loss and mitigating the negative effects. Combining optical and radar data has the potential to achieve faster detection of forest disturbance than using an individual system. The main challenge is the methodological approach for fusing the different datasets. The Fusion Near Real-Time (FNRT) algorithm (Tang et al. 2023) is a monitoring algorithm for tropical forest disturbance that combines data from Landsat, Sentinel-2, and Sentinel-1. For each sensor system, the FNRT algorithm fits time-series models over data from a three-year training period prior to the one-year monitoring period. FNRT then produces a change score for each new observation collected during the monitoring period based on the residuals and the root-mean-square error (RMSE) of the time-series model. The change scores from all three different sensors are then combined to form one dense time series for change detection. In this chapter, we will learn how to run a simplified version of FNRT to monitor forest disturbance in 2020 for a test area located in the Brazilian Amazon.

#### **50.2 Practicum**

This lab is designed for advanced users of Earth Engine. We assume that you already know how to import and preprocess large quantities of Landsat, Sentinel-2, and Sentinel-1 data as image collections. The code for importing and preprocessing input data will be provided in the example script, but it will not be discussed in detail in this chapter. To learn more about FNRT, refer to Tang et al. (2023).

#### **50.2.1 Section 1: Understand How FNRT Works**

For this lab, we will use a graphical user interface to help us understand conceptually how FNRT combines data from different sensors and detects forest disturbance. Figure 50.1 shows the user interface of the app.

**Code Checkpoint A35a**. The book's repository contains information about accessing the app.

**Fig. 50.1** User interface of FNRT app. The red box highlights a test site selected for demonstration of the algorithm. At the bottom left is a panel for showing plots of time-series data and model fit. The control menu is at the bottom right

There are three main tools that can be used to explore the FNRT algorithm. Note that only one plotting tool can be activated at a time, and each tool can be deactivated by clicking the button again:

Let's start by looking at some examples of time series so we have a better understanding of the model. We first need to know where to look for changes. Let's try to have only "Landsat" checked, and click **Run**. This will run FNRT with only Landsat data and quickly produce a forest disturbance map for 2020 (Fig. 50.2a). Now, click **Fit** to activate the "Fit time-series model" tool, and then click on a pixel that shows forest disturbance activity according to the map. Wait a few seconds to let the app generate the time-series plot in the time-series panel. Once it is completed, you should see a plot similar to Fig. 50.2c.

**Fig. 50.2 a** Forest disturbance map of 2020 for the test area created using Landsat data only; **b**  Landsat images before and after the forest disturbance; **c** time series of Landsat NDFI and model fit; and **d** time series of change scores and change detection result

The model-fit plot (Fig. 50.2c) shows data from Landsat for the training period (blue) and monitoring period (red). It also shows a harmonic time-series model fit to the training data (orange line). The model fit extends into the monitoring period, which represents the model prediction. Notice that the satellite observations (red points) are very close to the model fit before the forest disturbance occurs. After the disturbance, the observations deviate away from the model fit. The algorithm monitors change by keeping track of differences between the observations and the model predictions.

Now click **Monitor** to activate the "near real-time monitoring" tool, and click the same pixel in the map. The app will then make a plot similar to (d) in Fig. 50.2. The near real-time monitoring plot (Fig. 50.2d) shows the change scores (*Z*-score) for the training and the monitoring period. A large change score indicates that something looks different in the satellite image. But sometimes that can be caused by a cloud or a cloud shadow that slipped through the masking process (e.g., the large change scores in the training period). Therefore, the monitoring algorithm looks for multiple large scores within a short monitoring window as an indicator of real forest disturbance.

The monitoring algorithm works through all the observations available during the monitoring period chronologically. Using baseball as an analogy, we flag each change score above a threshold as a "strike" and each change score below the threshold as a "ball." If the algorithm finds multiple strikes in a short monitoring window, a "strikeout" would then be called on the pixel and the date of the first strike would be recorded as the date of change (Fig. 50.2d).

**Question 1**. How many strikes were detected before a forest disturbance was confirmed (strikeout)?

When performing multi-sensor data fusion, FNRT would fit a separate timeseries model and calculate change scores using data from each sensor (Fig. 50.3a– c). The change scores calculated using data from the three sensors are then combined, and change monitoring is applied to the combined scores (Fig. 50.4d).

The FNRT app also allows us to run change monitoring for the entire test area by clicking the **Run** button. A forest disturbance map showing the date of the disturbances is then added to the map. The users can use the checkbox widgets to choose which data stream is included in the process of making the map (e.g., Fig. 50.4). In general, optical data is more sensitive to forest disturbance than SAR data. Combining multiple data streams can help detect the disturbance earlier, compared to using data from each individual sensor.

**Question 2**. Can you find forest disturbances that were detected earlier through the use of data from all three sensors rather than using just Landsat data?

Now that we understand conceptually how FNRT works, let's create a simple script to implement it.

**Fig. 50.3** Time-series data and model fit using **a** Landsat, **b** Sentinel-2, and **c** Sentinel-1 data; **d**  shows the combined time series of change scores from all three sensors, and the change detection result

#### **50.2.2 Section 2: Define Study Area and Model Parameters**

#### **<IP_ADDRESS> Section 2.1: Study Area**

There are several ways to define a study area: (1) import an existing FeatureCollection from an Earth Engine asset; (2) use the drawing tool of the Code Editor; and (3) create a polygon by entering the coordinates. Here, we select a box, approximately 20 km by 20 km, located in the Brazilian Amazon as our study area.

**Fig. 50.4** Forest disturbance map made using data from **a** Landsat, **b** Sentinel-2, **c** Sentinel-1, and **d** all three sensors

```
var testArea = ee.Geometry.Polygon(
    [
        [
            [-66.73156878460787, -8.662236005089952],
            [-66.73156878460787, -8.916025640576244],
            [-66.44867083538912, -8.916025640576244],
            [-66.44867083538912, -8.662236005089952]
        ]
    ]);
Map.centerObject(testArea);
```
Note that because of the complexity of the FNRT algorithm, it requires quite a lot of computational resources and can result in a "User Memory Limit" error if applied to too large of an area. Divide the study area into multiple tiles if you are trying to monitor a large area.

#### **<IP_ADDRESS> Section 2.2: Model Parameters**

A few model parameters need to be defined and can be tuned to optimize the algorithm for specific regions or monitoring conditions.

```
// Start and end of the training and monitoring period.
var trainPeriod = ee.Dictionary({
    'start': '2017-01-01',
    'end': '2020-01-01'
});
var monitorPeriod = ee.Dictionary({
    'start': '2020-01-01',
    'end': '2021-01-01'
});
// Near-real-time monitoring parameters.
var nrtParam = {
    z: 2,
    m: 5,
    n: 4
};
// Sensor specific parameters.
var lstParam = {
    band: 'NDFI',
    minRMSE: 0.05,
    strikeOnly: false
};
var s2Param = {
    band: 'NDFI',
    minRMSE: 0.05,
    strikeOnly: false
};
var s1Param = {
    band: 'VV',
    minRMSE: 0.01,
    strikeOnly: true
};
```
For demonstration of the algorithm, we define the year 2020 as our monitoring period (monitorPeriod). We then recommend using the three years prior to the monitoring period as the training period (trainPeriod) to establish baseline time-series models.

A few parameters can be used to control the sensitivity of the algorithm: a *z*threshold (z) defines the minimum change score required to flag a possible change observation; an *m*-size (m) defines the size of the monitoring window, or how many consecutive observations we check for change; and an *n*-change (n) defines the number of possible change observations within the monitoring window required to confirm a forest disturbance.

For each sensor, we also need to select the band (band) to use for monitoring; a minimum RMSE (minRMSE) to prevent overly sensitive detection caused by unusually good model fit; and a strike-only flag (strikeOnly) to allow a data stream to contribute only to confirmation of a change.

**Code Checkpoint A35b**. The book's repository contains a script that shows what your code should look like at this point.

#### **50.2.3 Section 3: Import and Preprocess Data**

We need to import data from each sensor for the training period and the monitoring period as two separate image collections. Each image in the ImageCollection needs to be preprocessed. For optical data, we need to apply spectral unmixing and calculate the Normalized Difference Fraction Index (NDFI) (Souza et al. 2005). Since spectral unmixing would be applied to both Landsat and Sentinel-2 data, let's first implement a common "unmixing" function.

Several studies (e.g., Bullock et al. 2020; Chen et al. 2021) have shown that spectral unmixing and NDFI work very well in detecting tropical forest disturbance. Here, we use the same endmembers used in Souza et al. (2005) for the spectral unmixing and calculate NDFI based on Eqs. 4 and 5 in Souza et al. (2005).

```
var unmixing = function(col) {
    // Define endmembers and cloud fraction threshold.
    var gv = [500, 900, 400, 6100, 3000, 1000];
    var npv = [1400, 1700, 2200, 3000, 5500, 3000];
    var soil = [2000, 3000, 3400, 5800, 6000, 5800];
    var shade = [0, 0, 0, 0, 0, 0];
    var cloud = [9000, 9600, 8000, 7800, 7200, 6500];
    var cfThreshold = 0.05;
    return col.map(function(img) {
        // Select the spectral bands and perform unmixing
        var unmixed = img.select(['Blue', 'Green', 'Red',
                'NIR',
                'SWIR1', 'SWIR2'
            ])
            .unmix([gv, shade, npv, soil, cloud], true,
                true)
            .rename(['GV', 'Shade', 'NPV', 'Soil',
                'Cloud'
            ]);
        // Calculate Normalized Difference Fraction Index.+
        var NDFI = unmixed.expression(
            '10000 * ((GV / (1 - SHADE)) - (NPV + SOIL)) /
' +
            '((GV / (1 - SHADE)) + (NPV + SOIL))', {
                'GV': unmixed.select('GV'),
                'SHADE': unmixed.select('Shade'),
                'NPV': unmixed.select('NPV'),
                'SOIL': unmixed.select('Soil')
            }).rename('NDFI');
        // Mask cloudy pixel.
        var maskCloud = unmixed.select('Cloud').lt(
            cfThreshold);
        // Mask all shade pixel.
        var maskShade = unmixed.select('Shade').lt(1);
        // Mask pixel where NDFI cannot be calculated.
        var maskNDFI = unmixed.expression(
            '(GV / (1 - SHADE)) + (NPV + SOIL)', {
                'GV': unmixed.select('GV'),
                'SHADE': unmixed.select('Shade'),
                'NPV': unmixed.select('NPV'),
                'SOIL': unmixed.select('Soil')
            }).gt(0);
```

```
// Scale fractions to 0-10000 and apply masks.
        return img
            .addBands(unmixed.select(['GV', 'Shade',
                     'NPV', 'Soil'
                ])
                .multiply(10000))
            .addBands(NDFI)
            .updateMask(maskCloud)
            .updateMask(maskNDFI)
            .updateMask(maskShade);
    });
};
```
For Landsat data, we need to load the Landsat 7 and Landsat 8 Surface Reflectance product as an ImageCollection. We can filter the collection based on our study area and time period. Each image in the ImageCollection is masked based on the QA band; then spectral unmixing is applied and NDFI calculated. The loadLandsatData function is provided for importing and preprocessing Landsat data.

Preprocessing of Sentinel-2 data is similar to that of Landsat data because they are both optical remote sensing data and are almost identical in many ways. We need to load the Sentinel-2 top-of-atmosphere product (surface reflectance data do not have long enough time series for our purpose) and the Sentinel-2 cloud probability product as image collections, and filter based on our study area and time period. We then join the two collections. Each image in the collection is then masked based on both the QA band and the cloud probability.

For Sentinel-1 data, we use the Sentinel-1 SAR GRD product. Preprocessing of Sentinel-1 data includes: (1) calculation of backscattering; (2) calculation of the ratio of VH and VV; (3) calculation of three-pixel spatial mean; (4) choice of the orbital pass with the most data for the region of interest; and (5) radiometric slope correction using volume model. Please refer to Mullissa et al. (2021) for more detail on the preprocessing of SAR data.

Here, we assume you have already learned how to load and preprocess large quantities of Landsat, Sentinel-2, and Sentinel-1 data. We will not discuss these steps in detail here in this lab. Instead, three functions (loadLandsatData, loadS2Data, and loadS1Data) are provided in a separate script with which we can import using the require function. Now, we can load all the required input data. We also need to apply the spectral unmixing and calculate NDFI for the Landsat and Sentinel-2 data.

```
var input = require(
    'projects/gee-edu/book:Part A - Applications/A3 -
Terrestrial Applications/A3.5 Deforestation Viewed from
Multiple Sensors/modules/Inputs'
);
var lstTraining = unmixing(input.loadLandsatData(testArea,
    trainPeriod));
var lstMonitoring =
unmixing(input.loadLandsatData(testArea,
    monitorPeriod));
var s2Training = unmixing(input.loadS2Data(testArea,
trainPeriod));
var s2Monitoring = unmixing(input.loadS2Data(testArea,
    monitorPeriod));
var s1Training = input.loadS1Data(testArea, trainPeriod);
var s1Monitoring = input.loadS1Data(testArea,
monitorPeriod);
```
In addition to the remote sensing data, we also need to create a forest mask to limit our monitoring to forested areas. Here, we will use the Hansen Global Forest Change dataset. We first use the "tree canopy cover" layer to find areas that had larger than 50% canopy cover in 2000. We then remove all the forest loss area during 2001–2019 and add in all the forest gain from 2000 to 2012 to create a mask of forested area by the beginning of our monitoring period.

```
var hansen =
ee.Image('UMD/hansen/global_forest_change_2020_v1_8')
    .unmask();
var forestMask = hansen.select('treecover2000')
    .gt(50)
    .add(hansen.select('gain'))
    .subtract(hansen.select('loss'))
    .add(hansen.select('lossyear')
        .eq(20))
    .gt(0)
    .clip(testArea);
```
Before we move on to the next section, let's do some simple checking. Print out the image collections and check whether the correct data is loaded. You can also load the forest mask and visually examine the quality of the mask.

```
var maskVis = {
    min: 0,
    max: 1,
    palette: ['blue', 'green']
};
Map.addLayer(forestMask, maskVis, 'Forest Mask');
print('lstTraining', lstTraining);
print('lstMonitoring', lstMonitoring);
print('s2Training', s2Training);
print('s2Monitoring', s2Monitoring);
print('s1Training', s1Training);
print('s1Monitoring', s1Monitoring);
```
**Question 3**. How many images are available from each sensor for our monitoring period? Are you surprised by the numbers? (Note that one sensor may have more data because the study area happened to be in the side-lap zone where two orbit passes overlap.) What bands are included in the data of each sensor?

**Code Checkpoint A35c**. The book's repository contains a script that shows what your code should look like at this point.

#### **50.2.4 Section 4: Establish Baseline Time-Series Model**

Before we can monitor change in near real time, we need to establish baseline time-series models using the data of the training period, so that we know what to expect from new observations in the monitoring period. Due to the different physical nature of optical and radar data, it is easier to fit a separate time-series model to data from each sensor. Here, we use a harmonic time-series model similar to the one used in the Continuous Change Detection and Classification (CCDC) algorithm (Zhu and Woodcock 2014). Because the timestamps of all the input images are stored in the unit of milliseconds while the harmonic model prefers a unit of years, we first need to define a function to convert any date object into the unit of fractional year (e.g., 2015-03-01 to 2015.1612).

```
var toFracYear = function(date) {
    var year = date.get('year');
    var fYear = date.difference(
        ee.Date.fromYMD(year, 1, 1), 'year');
    return year.add(fYear);
};
```
Here, we define a function (fitHarmonicMode) to fit a harmonic model to an ImageCollection of time-series data. We first construct the dependent variables according to Eq. 1 in Zhu and Woodcock (2014) and add them to each of the input images as new bands (using addDependents function). Then, we use a reducer (ee.Reducer.robustLinearRegression) to fit the model to the data using robust linear regression. The function returns the model coefficients and the model RMSE.

```
var fitHarmonicModel = function(col, band) {
    // Function to add dependent variables to an image.
    var addDependents = function(img) {
        // Transform time variable to fract ional year.
        var t = ee.Number(toFracYear(
            ee.Date(img.get('system:time_start')), 1));
        var omega = 2.0 * Math.PI;
        // Construct dependent variables image.
        var dependents = ee.Image.constant([
                1, t,
                t.multiply(omega).cos(),
                t.multiply(omega).sin(),
                t.multiply(omega * 2).cos(),
                t.multiply(omega * 2).sin(),
                t.multiply(omega * 3).cos(),
                t.multiply(omega * 3).sin()
            ])
            .float()
            .rename(['INTP', 'SLP', 'COS', 'SIN',
                 'COS2', 'SIN2', 'COS3', 'SIN3'
            ]);
        return img.addBands(dependents);
    };
```

```
// Function to add dependent variable images to a ll
images.
    var prepareData = function(col, band) {
        return ee.ImageCollection(col.map(function(img) {
             return addDependents(img.select(band))
                 .select(['INTP', 'SLP', 'COS',
                      'SIN',
                      'COS2', 'SIN2', 'COS3',
                      'SIN3',
                     band
                 ])
                 .updateMask(img.select(band)
                     .mask());
        }));
    };
    var col2 = prepareData(col, band);
    // Fit model to data using robust linear regression. // Fit linear regression.
    var ccd = col2
        .reduce(ee.Reducer.robustLinearRegression(8, 1), 4)
        .rename([band + '_coefs', band + '_rmse']);
    // Return model coefficients and model rmse.
    return ccd.select(band + '_coefs').arrayTranspose()
        .addBands(ccd.select(band + '_rmse'));
};
```
With this function (fitHarmonicModel), we can fit the harmonic timeseries model to the training data of all three sensors. It is always a good practice to also add some metadata to the results for future reference. The model fitting process will take a few minutes and quite a lot of computer memory. Therefore, it is recommended that we run them as tasks and save the results as assets before we proceed to the next step.

```
// Fit harmonic models to training data of all sensors.
var lstModel = fitHarmonicModel(lstTraining, lstParam.band)
    .set({
        region: 'test',
        sensor: 'Landsat'
    });
var s2Model = fitHarmonicModel(s2Training, s2Param.band)
    .set({
        region: 'test',
        sensor: 'Sentinel-2'
    });
var s1Model = fitHarmonicModel(s1Training, s2Param.band)
    .set({
        region: 'test',
        sensor: 'Sentinel-1'
    });
// Define function to save the results.
var saveModel = function(model, prefix) {
    Export.image.toAsset({
        image: model,
        scale: 30,
        assetId: prefix + '_CCD',
        description: 'Save_' + prefix + '_CCD',
        region: testArea,
        maxPixels: 1e13,
        pyramidingPolicy: {
            '.default': 'sample'
        }
    });
};
// Run the saving function.
saveModel(lstModel, 'LST');
saveModel(s2Model, 'S2');
saveModel(s1Model, 'S1');
```
Run the script. The model fitting results in export tasks to create assets. Instead of exporting the results for this exercise, however, you should proceed to the next section. In that section, we will continue with a precomputed asset that is the same as what the above tasks would create.

**Code Checkpoint A35d**. The book's repository contains a script that shows what your code should look like at this point.

#### **50.2.5 Section 5: Create Predicted Values for the Monitoring Period**

We will now start exploring the pre-exported results mentioned in the previous section. Place this code below the checkpoint of the previous section:

```
var models = ee.ImageCollection('projects/gee-
book/assets/A3-5/ccd');
var lstModel = models
    .filterMetadata('sensor', 'equals', 'Landsat').first();
var s2Model = models
    .filterMetadata('sensor', 'equals', 'Sentinel-
2').first();
var s1Model = models
    .filterMetadata('sensor', 'equals', 'Sentinel-
1').first();
```
In this section, we will use the time-series models that we produced in Sect. 50.2.4 to create a predicted image (also commonly referred to as a synthetic image; see Zhu et al. 2015) for any given date. The predicted image can give us a good estimate of the reflectance (or index such as NDFI) value of a place assuming no change has occurred. We can then compare the predicted image to the image collected by the satellite on the same date. Areas that look very different would likely be areas that have changed. Note that all the coefficients of the harmonic model were saved as a single-band array. So the very first step here is to define a function dearrayModel to convert the array image into a multiband image with each coefficient in one band.

```
var dearrayModel = function(model, band) {
    band = band + '_';
    // Function to extract a non-harmonic coefficients.
    var genCoefImg = function(model, band, coef) {
        var zeros = ee.Array(0).repeat(0, 1);
        var coefImg = model.select(band + coef)
            .arrayCat(zeros, 0).float()
            .arraySlice(0, 0, 1);
        return ee.Image(coefImg
            .arrayFlatten([
                [ee.String('S1_')
                    .cat(band).cat(coef)
                ]
            ]));
    };
    // Function to extract harmonic coefficients.
    var genHarmImg = function(model, band) {
        var harms = ['INTP', 'SLP', 'COS', 'SIN',
            'COS2', 'SIN2', 'COS3', 'SIN3'
        ];
        var zeros = ee.Image(ee.Array([
                ee.List.repeat(0, harms.length)
            ]))
            .arrayRepeat(0, 1);
        var coefImg = model.select(band + 'coefs')
            .arrayCat(zeros, 0).float()
            .arraySlice(0, 0, 1);
        return ee.Image(coefImg
            .arrayFlatten([
                [ee.String(band).cat('coef')], harms
            ]));
    };
    // Extract harmonic coefficients and rmse.
    var rmse = genCoefImg(model, band, 'rmse');
    var coef = genHarmImg(model, band);
    return ee.Image.cat(rmse, coef);
};
```
Second, we need to define a function that creates a predicted image for all the dates for which a real image is available during the monitoring period.

```
var createPredImg = function(modelImg, img, band, sensor) {
    // Reformat date.
    var date =
toFracYear(ee.Date(img.get('system:time_start')));
    var dateString = ee.Date(img.get('system:time_start'))
        .format('yyyyMMdd');
    // List of coefficients .
    var coefs = ['INTP', 'SLP', 'COS', 'SIN',
        'COS2', 'SIN2', 'COS3', 'SIN3'
    ];
    // Get coefficients images from model image.
    var coef = ee.Image(coefs.map(function(coef) {
        return modelImg.select(".*".concat(coef));
    })).rename(coefs);
    var t = ee.Number(date);
    var omega = 2.0 * Math.PI;
    // Construct dependent variables.
    var pred = ee.Image.constant([
            1, t,
            t.multiply(omega).cos(),
            t.multiply(omega).sin(),
            t.multiply(omega * 2).cos(),
            t.multiply(omega * 2).sin(),
            t.multiply(omega * 3).cos(),
            t.multiply(omega * 3).sin()
        ])
        .float();
    // Matrix multiply dependent variables with
coefficients.
    return pred.multiply(coef).reduce('sum')
        // Add original image and rename bands.
        .addBands(img, [band]).rename(['predicted', band])
        // Preserve some metadata.
        .set({
            'sensor': sensor,
            'system:time_start':
img.get('system:time_start'),
            'dateString': dateString
        });
};
```
The createPredImg function calculates a predicted image based on the date of the real image and the model coefficients. It returns a new image with both the predicted image and the real image.

We can then apply the createPredImg function to each image in the image collections of the monitoring data. We can define another function to map over an ImageCollection to apply createPredImg to all images in the collection.

```
var addPredicted = function(data, modelImg, band, sensor) {
    return ee.ImageCollection(data.map(function(img) {
        return createPredImg(modelImg, img, band,
            sensor);
    }));
};
```
Let's apply addPredicted to monitoring data of all three sensors. Note this process is still sensor-specific. Let's print out the result for Landsat to examine the structure of the output.

```
// Convert models to non-array images.
var lstModelImg = dearrayModel(lstModel, lstParam.band);
var s2ModelImg = dearrayModel(s2Model, s2Param.band);
var s1ModelImg = dearrayModel(s1Model, s1Param.band);
// Add predicted image to each real image.
var lstPredicted = addPredicted(lstMonitoring, lstModelImg,
    lstParam.band, 'Landsat');
var s2Predicted = addPredicted(s2Monitoring, s2ModelImg,
    s2Param.band, 'Sentinel-2');
var s1Predicted = addPredicted(s1Monitoring, s1ModelImg,
    s1Param.band, 'Sentinel-1');
print('lstPredicted', lstPredicted);
```
**Question 4**. How many images are in the ImageCollection of the predicted images? Is it the same number as the original ImageCollection of the monitoring data?

**Code Checkpoint A35e**. The book's repository contains a script that shows what your code should look like at this point.

#### **50.2.6 Section 6: Calculate Change Scores**

Assuming a good model fit, the predicted image would be very close to the real image of the same date for areas with no change. So what we need to do now is examine each pair of real and predicted images and look for areas where they differ from each other. And if we merge the three data streams, we would be able to significantly increase the data density and potentially detect changes faster. But before we can fuse the three data streams together, we need to produce a normalized measurement of the likelihood of change.

The FNRT algorithm calculates a change score (Tang et al. 2019) very similar to a *z*-score. First, it calculates the residuals defined as the difference between the predicted value and the observed value. The residuals are then scaled with the RMSE of the time-series model to get the change score. If a change score is above the preset threshold (z in the parameters), then it is flagged as "possible change" (or a "strike" in baseball terminology). The implementation of this part is quite simple:

```
// Function to calculate residuals.
var addResiduals = function(data, band) {
    return ee.ImageCollection(data.map(function(img) {
        return img.select('predicted')
            // Restrict predicted value to be under 10000
            .where(img.select('predicted').gt(10000),
                10000)
            // Calculate the residual
            .subtract(img.select(band))
            .rename('residual')
            // Save some metadata
            .set({
                 'sensor': img.get('sensor'),
                 'system:time_start': img.get(
                     'system:time_start'),
                 'dateString': img.get(
                     'dateString')
            });
    }));
};
// Function to calculate change score and flag change.
var addChangeScores = function(data, rmse, minRMSE,
    threshold, strikeOnly) {
    // If strikeOnly then we need a mask for balls.
    var mask = ee.Image(0);
    if (strikeOnly) {
        mask = ee.Image(1);
    }
```

```
return ee.ImageCollection(data.map(function(img) {
        // Calculate change score
        var z = img.divide(rmse.max(minRMSE));
        // Check if score is above threshold
        var strike = z.multiply(z.gt(threshold));
        // Create the output image.
        var zStack = ee.Image.cat(z, strike).rename([
                 'z', 'strike'
            ])
            .set({
                 'sensor': img.get('sensor'),
                 'system:time_start': img.get(
                     'system:time_start')
            });
        // Mask balls if strikeOnly.
        return zStack.updateMask(strike.gt(0).or(
            mask));
    }));
};
```
Note that here, we restrict the predicted value to be under 10,000 because that is the theoretical maximum for NDFI.

We also need to have a minimum RMSE (minRMSE) when scaling the residuals. This is because when using indices such as NDFI, it is sometimes possible to have perfect model fits resulting in very low or even zero RMSE, which would then cause the model to be overly sensitive to subtle variations in the monitoring data. The minRMSE is calculated for each pixel as a percentage of the average values of all observations in the training period.

For this particular implementation of FNRT, the determination of a strike is directional, meaning that we are looking for changes only in one direction. This is because we already expect that a forest disturbance will cause a decrease in NDFI or VV backscattering. Modification would be required if monitoring in another direction (e.g., if using SWIR reflectance) or bidirectional monitoring is desired. Let's apply the two new functions to our three data streams:

```
// Add residuals to collection of predicted images.
var lstResiduals = addResiduals(lstPredicted,
lstParam.band);
var s2Residuals = addResiduals(s2Predicted, s2Param.band);
var s1Residuals = addResiduals(s1Predicted, s1Param.band);
// Add change score to residuals.
var lstScores = addChangeScores(
    lstResiduals, lstModelImg.select('.*rmse'),
    lstPredicted.select(lstParam.band).mean()
    .abs().multiply(lstParam.minRMSE),
    nrtParam.z, lstParam.strikeOnly);
var s2Scores = addChangeScores(
    s2Residuals, s2ModelImg.select('.*rmse'),
    s2Predicted.select(s2Param.band).mean()
    .abs().multiply(s2Param.minRMSE),
    nrtParam.z, s2Param.strikeOnly);
var s1Scores = addChangeScores(
    s1Residuals, s1ModelImg.select('.*rmse'),
    s1Predicted.select(s1Param.band).mean()
    .abs().multiply(s1Param.minRMSE),
    nrtParam.z, s1Param.strikeOnly);
print('lstScores', lstScores);
```
**Code Checkpoint A35f**. The book's repository contains a script that shows what your code should look like at this point.

**Question 5**. What bands are included in the change scores images?

#### **50.2.7 Section 7: Multisensor Data Fusion and Change Detection**

Now, we are finally ready to fuse the three data streams together. This step is quite straightforward: we just need to merge the three image collections of the change scores, and then sort them by their timestamps.

```
var fused = lstScores.merge(s2Scores).merge(s1Scores)
    .sort('system:time_start');
```
The change monitoring function uses a binary array to keep track of the number of strikes within a monitoring window. The size of the monitoring window is defined by the model parameters (nrtParam.m). The default monitoring window size is five.

We first initiate the monitoring window with a binary array of 0 with a size of five bits (00000) for each pixel. When we iterate through all the change score images chronologically, we shift the binary array left for those pixels that were not masked in the current change score image (00000→0000). We then append the current change flag (strike or not) to the shifted binary array (0000 →00001). The process continues, as we iterate through all the change scores images (e.g., 00001→00011 →00110→01101→11011). Every time a new value is appended to the binary array, we also check how many strikes (1s) are there. If the number of strikes exceeds the required number (nrtParam.n) to flag a change, then a change is flagged and the date of the change recorded.

```
var monitorChange = function(changeScores, nrtParam) {
    // Initialize an empty image.
    var zeros = ee.Image(0).addBands(ee.Image(0))
        .rename(['change', 'date']);
    // Determine shift size based on size of monitoring
window.
    var shift = Math.pow(2, nrtParam.m - 1) - 1;
    // Function to monitor.
    var monitor = function(img, result) {
        // Retrieve change image at last step.
        var change = ee.Image(result).select('change');
        // Retrieve change date image at last step.
        var date = ee.Image(result).select('date');
        // Create a shift image to shift the change binary
array
        // left for one space so that new one can be
appended.
        var shiftImg = img.select('z').mask().eq(0)
            .multiply(shift + 1).add(shift);
        change = change.bitwiseAnd(shiftImg)
            .multiply(shiftImg.eq(shift).add(1))
            .add(img.select('strike').unmask().gt(0));
```

```
// Check if there are enough strike in the current
        // monitoring window to flag a change.
        date = date.add(change.bitCount().gte(nrtParam.n)
            // Ignore pixels where change already detected.
            .multiply(date.eq(0))
            // Record change date where change is flagged.
            .multiply(ee.Number(toFracYear(
                ee.Date(img.get(
                    'system:time_start')), 1))));
        // Combine change and date layer for next
iteration.
        return (change.addBands(date));
    };
    // Iterate through the time series and look for change.
    return ee.Image(changeScores.iterate(monitor, zeros))
        // Select change date layer and selfmask.
        .select('date').rename('Alerts').selfMask();
};
```
Now, we apply the monitoring function (monitorChange) to the fused data (fused) using the default model parameters (nrtParam). We can then add the result, the forest disturbance map, to the map panel to visually examine it.

```
var alerts = monitorChange(fused,
nrtParam).updateMask(forestMask);
print('alerts', alerts);
// Define a visualization parameter.
var altVisParam = {
    min: 2020.4,
    max: 2021,
    palette: ['FF0080', 'EC1280', 'DA2480', 'C83680',
'B64880',
        'A35B80', '916D80', '7F7F80', '6D9180', '5BA380',
        '48B680', '36C880', '24DA80', '12EC80', '00FF80',
        '00EB89', '00D793', '00C49D', '00B0A7', '009CB0',
        '0089BA', '0075C4', '0062CE', '004ED7', '003AE1',
        '0027EB', '0013F5', '0000FF'
    ]
};
Map.centerObject(testArea, 10);
Map.addLayer(alerts, altVisParam, 'Forest Disturbance Map
(2020)');
Map.setOptions('SATELLITE');
```
**Code Checkpoint A35g**. The book's repository contains a script that shows what your code should look like at this point.

If we run the whole script now, we can create a map of forest disturbance that occurred in 2020 for our test area. The visualization of the map shows us the timing of all the forest disturbance detected by the FNRT algorithm (Fig. 50.5). You will notice a few large disturbance events. There may also be some tiny speckles, which are likely errors of commission. Try adjusting the monitoring parameters (nrtParam) to see if you can improve the map. Another common practice is to apply a spatial filter to filter out scattered single pixels of change (Bullock et al. 2020).

**Question 6**. Can you run the change detection without Sentinel-1 input or with just Landsat input? Was forest disturbance detected faster when we used data from all three sensors?

**Fig. 50.5** Map of forest disturbance occurring in 2020 for a selected test area located on the frontline of deforestation in the Brazilian Amazon

#### **50.3 Synthesis**

Now that you have learned how to implement the simplified version of FNRT for the test area that we selected, let's try to use it to do some real monitoring of tropical forest disturbance.

**Assignment 1**. Modify the code and apply it in a different area of your choice. Note that you will need to keep the area small in order to be able to run everything easily in the Code Editor.

**Assignment 2**. Try to monitor more recent changes by changing the monitoring period to 2021. Note you would need to change the training period to 2018–2020, and also update the forest mask.

**Assignment 3**. If you want to try a more challenging task, you can integrate your results into the FNRT app that we used in Sect. 50.2.1.

#### **50.4 Conclusion**

In this chapter, we experimented with a simple way to combine data from three different sensors (Landsat, Sentinel-2, and Sentinel-1) for the purpose of detecting tropical forest disturbance as early as possible. We started with loading and preprocessing the input data. We learned how to fit a harmonic model to timeseries data, and how to make predicted observations based on a model fit. We also learned how to look for changes with a moving monitoring window over a time series of change scores. The full FNRT algorithm is more complex than the "lite" version that we implemented in this lab. If interested, please refer to Tang et al. (2023) for more details.

#### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **51 Working with GPS and Weather Data**

Peder Engelstad , Daniel Carver , and Nicholas E. Young

#### **Overview**

The purpose of this chapter is to demonstrate how to use Google Earth Engine as a means of associating remotely sensed data (weather observations) with open-source GPS point locations. These methods will provide a quick and easy way to access and analyze large amounts of information relative to your own research and efficiently move your data outside Earth Engine.

#### **Learning Outcomes**

#### **Helps if you know how to**

P. Engelstad (B) · N. E. Young Natural Resource Ecology Laboratory, Colorado State University, Fort Collins, CO, USA e-mail<EMAIL_ADDRESS>
N. E. Young e-mail<EMAIL_ADDRESS>
D. Carver Geospatial Centroid, Colorado State University, Fort Collins, CO, USA e-mail<EMAIL_ADDRESS>
#### **51.1 Introduction to Theory**

Knowing how animals interact with their environment is critical to understanding how to manage them. The choices animals make are influenced by basic survival needs (e.g., food, shelter, water) and dynamic factors such as local weather conditions. Without direct observation, it can be difficult to understand the relationship between animal movement and weather conditions.

In this chapter, we will explore information retrieved from a GPS collar worn by a cougar and explore its relationship to daily temperature estimates from the Daymet climatological dataset available in Google Earth Engine. This will require us to bring an asset into Earth Engine, connect the weather values to the point locations, and bring those value-added data back out of Earth Engine for further analysis.

#### **51.2 Practicum**

#### **51.2.1 Section 1: GPS Location Data**

Using GPS collars, Mahoney and Young (2017) tracked the movement of 2 cougars and 16 coyotes in Central Utah. These data were used to understand some of the behavioral patterns of the individuals. These data have been freely shared with the broader research community and the public through Movebank, an online repository for animal movement datasets from across the globe. While some Movebank datasets list only the contact information of the authors, others (like those from the Mahoney study) allow you to display and download the information via an interactive web map.

We will pair the data on cougar movement with Daymet weather data. According to the Daymet website, the dataset "provides gridded estimates of daily weather parameters. Seven surface weather parameters are available at a daily time step, 1 km ×1 km spatial resolution, with a North American spatial extent allowing a rich resource of daily surface meteorology."

With data for every day at a 1 km2 spatial resolution, the Daymet data are a great resource for the temporal and spatial scales at which a cougar would interact with the landscape. There are seven measured values in total, allowing us to check multiple aspects of the weather to assess how it may be affecting behavior (Fig. 51.1).

#### **51.2.2 Section 2: Bringing Data into Earth Engine**

In this chapter, we will discuss how to import assets into Earth Engine, extract values from a dataset, and export those values out of Earth Engine. The processes by which you can bring data into Earth Engine change frequently, so it is best to go directly to the documentation to view the latest updates.

**Fig. 51.1** Metadata for Daymet imagery within Earth Engine

#### **<IP_ADDRESS> Section 2.1: Bringing in an Asset**

Using the following code, start a fresh script and import Mahoney's cougar movement data from Movebank (which has already been uploaded for you in the assets of this book). Here, we will focus on a single animal (cougar ID F53). Because the original dataset was in CSV format, it has been converted into a shapefile outside of GEE. During this process, it is important to note that data have been projected to the WGS 1984 (EPSG: 4326) coordinate reference system. Projecting to EPSG: 4326 is recommended in Earth Engine to minimize the reprojection errors during the uploading process.

```
// Import the data and add it to the map and print.
var cougarF53 = ee.FeatureCollection(
    'projects/gee-book/assets/A3-6/cougarF53');
Map.centerObject(cougarF53, 10);
Map.addLayer(cougarF53, {}, 'cougar presence data');
print(cougarF53, 'cougar data');
```
You can use the **Inspector** tool to look at the attribute data associated with the new asset. With the points visualized, make a geometry feature that encompasses our area of interest by selecting the square geometry tool and drawing a box that encompasses the points (Fig. 51.2). We will use the geometry feature to filter our climate data.

**Fig. 51.2** Draw a geometry feature around the points to spatially filter the climate data

#### **<IP_ADDRESS> Section 2.2: Defining Weather Variables**

In this chapter, we are using Earth Engine as a means of associating remotely sensed data (i.e., our rasters) with our point locations. While the process is conceptually straightforward, it does take some work to accomplish. With our points loaded, the next step is to import the Daymet weather variables.

We are using the NASA-derived dataset Daymet V4 due to its 1 km2 spatial resolution and the fact that it measures the environmental variables we think might be related to cougar movement or behavior. We will import these variables by calling the unique ID of the dataset, filtering it to our bounding box geometry and the dates during which the data were collected.

```
// Call in image collection and filter.
var Daymet = ee.ImageCollection('NASA/ORNL/DAYMET_V4')
    .filterDate('2014-02-11', '2014-11-02')
    .filterBounds(geometry)
    .map(function(image) {
        return image.clip(geometry);
    });
print(Daymet, 'Daymet');
```
**Fig. 51.3** View of the structure of the Daymet V4 data from the print statement

From the print statement (Fig. 51.3), we can see that this is an ImageCollection with 264 images (though your total number of images may be different, as the dataset changes over time). Each image has seven bands relating to specific weather measurements (see Fig. 51.1). Now that both datasets are loaded, we will associate the cougar occurrences data with the weather data.

#### **<IP_ADDRESS> Section 2.3: Extracting Values**

With our points and imagery loaded, we can call a function to extract values from the underlying raster based on the known locations of the cougar. We will do this using the ee.Image.sampleRegions function. Search for the ee.Image.sampleRegions function under the **Docs** tab to familiarize yourself with the parameters it requires.

If we tried to call this function on the Daymet ImageCollection, we would get an error, because ee.Image.sampleRegions is a function of an image. To get around this, we will convert the Daymet ImageCollection into a multiband image (Fig. 51.4). Each of the seven measurements for each day will become a specific band in our multiband image. This process will help us in the end, because each band is defined by the date it was collected and the variable it shows. We can use this information to determine which data connects to the positions of the cougar on a specific day.

**Fig. 51.4** Print statement showing the result of converting the Daymet ImageCollection into a multiband image

It is important to note that, with many images in the ImageCollection, we are going to create a single image with a large number of bands. Because Earth Engine is good at data manipulations, it can handle this type of request.

```
// Convert to a multiband image.
var DaymetImage = Daymet.toBands();
print(DaymetImage, 'DaymetImage');
```
Now that we have a single multiband image, we can use the sampleRegions function to extract information across multiple dates at each of our GPS point locations (Fig. 51.5). There are three parameters of this function that should be considered:

**Fig. 51.5** Print statement from the sampleRegions function showing that our GPS point locations now have weather measurements associated with them

```
// Call the sample regions function.
var samples = DaymetImage.sampleRegions({
    collection: cougarF53,
    properties: ['id'],
    scale: 1000
});
print(samples, 'samples');
```
**Code Checkpoint A36a**. The book's repository contains a script that shows what your code should look like at this point.

#### **<IP_ADDRESS> Section 2.4: Exporting Exporting Points**

We now have a series of daily weather data associated with the known locations of the cougar known as F53. While we could work more with these data in Earth Engine, it will be easy to bring them into R, Python, or Excel for further analysis. There are a few options to define where the exported data will be created. Generally speaking, saving the data to a Google Drive account provides an easy way to access the data with another program. Here, we will use a dictionary, denoted by the curly brackets {}, to define the parameters of the Export.table.toDrive function.

We would have preferred to export to create a shapefile, but a shapefile can contain only 255 columns. The samples variable has 1361 columns, so we will export the data as a CSV file.

```
// Export value added data to your Google Drive.
Export.table.toDrive({
    collection: samples,
    description: 'cougarDaymetToDriveExample',
    fileFormat: 'csv'
});
```
When you export something, the **Tasks** panel will light up. To run the task, click the **Run** button (Fig. 51.6).

When you click the **Run** button, the pop-up shown in Fig. 51.7 will appear. This allows you to edit the details of the export.

#### **Exporting a Raster**

While working with these spatial data, you may have realized that a raster showing the median values over the time period when data were collected on the cougar could be useful information to have. To do this, we will apply a median reducer function to the Daymet ImageCollection to generate a median value for each parameter in each cell. As with the tabular data, we will export this multiband image to Google Drive. Once we convert the ImageCollection to an image using the median function, we can clip it to the geometry feature object. This feature will be exported as a multiband raster.

**Fig. 51.6** Example of the **Tasks** bar once a script with an export function has been run

**Fig. 51.7** Example of the user-defined parameters present when exporting a feature from Earth Engine

```
// Apply a median reducer to the dataset.
var daymet1 = Daymet
    .median()
    .clip(geometry);
print(daymet1);
// Export the image to drive.
Export.image.toDrive({
    image: daymet1,
    description: 'MedianValueForStudyArea',
    scale: 1000,
    region: geometry,
    maxPixels: 1e9
});
```
In Earth Engine, there are many options for exporting images. One of the most important options when exporting data is the maxPixels parameter. Generally speaking, Earth Engine will not allow you to export a raster with more than 109 pixels. With the maxPixels parameter, you can bump this up to around 1012 pixels per image. If you are exporting data for an area larger than 1012 pixels, you will need to be creative about how to get the information out of Earth Engine. Sometimes, this involves splitting the image into smaller pieces (as described in Chap. 29) or even reevaluating the usefulness of such a large image outside of Earth Engine.

**Code Checkpoint A36b**. The book's repository contains a script that shows what your code should look like at this point.

**Question 1**. With the maxPixels count set to 1e9, what square area would you expect to be able to export for imagery with 1000, 250, 30, 10 m cell size? Determine how this area measure changes when using 1012 pixels for one of the cell sizes to determine the percent increase.

**Question 2**. We selected all available variables from the Daymet collection. What specific bands from Daymet might be most relevant to cougars? How would you edit the code to select just one or a few bands?

**Question 3**. When we extracted the raster data, we set the scale parameter to 1000 to match the known spatial resolution of the Daymet dataset. Would it be reasonable to change that scale parameter to 500 to get final resolution data? What does Earth Engine tell you when you attempt to change the scale of an export?

### **51.3 Synthesis**

**Assignment 1**. Utilize the chapter's script with a different dataset to process tracking information from a different animal from Movebank or another occurrence database of your choice (e.g., GBIF.org). Evaluate whether the temperature ranges between the two animals are different and make a hypothesis about why this might be.

#### **51.4 Conclusion**

While Google Earth Engine can be used for planetary-scale analyses, it is also an effective resource for quickly accessing and analyzing large amounts of information across time at a scale relative to your own data. The method presented in this chapter is a great way to add value to your own research. Here, we worked with weather data, but that is by no means the only option. You can connect your data to many other datasets within Earth Engine. It is up to you to explore and find what is relevant to your work.

#### **Reference**

Mahoney PJ, Young JK (2017) Uncovering behavioural states from animal activity and site fidelity patterns. Methods Ecol Evol 8:174–183. https://doi.org/10.1111/2041-210X.12658

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **52 Creating Presence and Absence Points**

Peder Engelstad , Daniel Carver , and Nicholas E. Young

### **Overview**

The purpose of this chapter is to demonstrate a method to generate your own presence and absence data and distribute those samples using specific ecological characteristics found in remotely sensed imagery. You will see that even when field data is unavailable, you can still digitally sample a landscape and gather information on current or past ecological conditions.

#### **Learning Outcomes**

P. Engelstad (B) · N. E. Young

Natural Resource Ecology Laboratory, Colorado State University, Fort Collins, CO 80523-1499, USA

e-mail<EMAIL_ADDRESS>
N. E. Young e-mail<EMAIL_ADDRESS>
D. Carver

Geospatial Centroid, Colorado State University, Fort Collins, CO 80523-1499, USA e-mail<EMAIL_ADDRESS>
#### **Helps if you know how to**

#### **52.1 Introduction to Theory**

Herbivore grazing can have a negative effect on aspen regeneration in some areas, as aspens tend to grow in large monotypic stands (Halofsky and Ripple 2008). Excluding elk, deer, and cow grazing from an area has observable effects on aspen regrowth. In a hypothetical study, we may be interested in understanding the effect of herbivory from these species on various ecological aspects of aspen stands. But how can we monitor aspen stands without setting foot in the field? In this chapter, we will use multiple datasets and high-resolution resolution imagery (1 m2) to generate sampling locations for this theoretical field survey. We will also build a presence/absence dataset that could be used to train a spatial model of aspen coverage.

#### **52.2 Practicum**

The National Land Cover Database (NLCD) is a Landsat-derived land cover database with 30 m spatial resolution, available at multiple time periods. Data from 1992 are primarily based on an unsupervised classification of Landsat data. The other years rely on a decision tree classification to identify the land cover classes.

The National Agriculture Imagery Program (NAIP) acquires aerial imagery during the agricultural growing seasons across the continental United States. NAIP projects are contracted each year based upon available funding and the Farm Service Agency imagery acquisition cycle. NAIP imagery is acquired at a one-meter ground sample distance with a horizontal accuracy that matches within six meters of photo-identifiable ground control points. NAIP imagery is often collected in four bands: blue, green, red, and near infrared. The near infrared band is helpful in distinguishing between different types of vegetation.

The USGS National Elevation Dataset (NED) 1/3 arc-second is an elevation dataset produced by the USGS. This coast-to-coast dataset is available at 0.33 arc-second spatial resolution (30 m) across the lower 48 states, parts of Alaska, Hawaii, and US Territories.

**Fig. 52.1** Example of the general region of interest to use for this module

#### **52.2.1 Section 1: Developing Your Own Sampling Locations**

We will start by developing potential field sampling locations based on the relative physical and ecological conditions.

#### **<IP_ADDRESS> Section 1.1: Region of Interest**

The geographic region for this chapter is the Grand Mesa in western Colorado. Considered to be the largest mesa in the world by area, the Grand Mesa rises from 4000 feet near Grand Junction, Colorado, to over 10,000 feet at its highest point. Historically, the Grand Mesa was glaciated, leaving its top flattened and spotted with lakes. While heavily forested at higher elevation, there are distinct transitions between shrubland, aspen, and conifer forest along the steep slopes. As one of the westernmost extensions of the montane environment in Colorado, the Grand Mesa represents important habitat for many ungulate species.

Our first step is to open a new script in Earth Engine. First, create a region of interest that encompasses the Grand Mesa (you can search for it by name in the search bar at the top). Do so using the geometry tool. Once you create the feature, rename it roi (Fig. 52.1).

#### **<IP_ADDRESS> Section 1.2: Working with NAIP**

We will rely on NAIP imagery for multiple steps in this process. NAIP imagery is not collected every year, so it makes sense to load in multiple years to determine what time frame is available. We will use a print statement to see what years are available for our area from 2015 and 2017. Copy the following code to start your script.

```
// Call in NAIP imagery as an image collection.
var naip = ee.ImageCollection('USDA/NAIP/DOQQ')
    .filterBounds(roi)
    .filterDate('2015-01-01', '2017-12-31');
Map.centerObject(naip);
print(naip);
```
When you run the script you will see a print statement similar to the one shown in Fig. 52.2.

The print statement shows us that imagery is available for both 2015 and 2017, yet it is difficult to determine the extent of coverage present without visualizing the data. For now, we will add both collections to the map to see what the image availability looks like. Add the following code to your existing script.

**Fig. 52.2** Print statement identifying the years of available imagery for this region of interest

```
// Filter the data based on date.
var naip2017 = naip
    .filterDate('2017-01-01', '2017-12-31');
var naip2015 = naip
    .filterDate('2015-01-01', '2015-12-31');
// Define viewing parameters for multi band images.
var visParamsFalse = {
    bands: ['N', 'R', 'G']
};
var visParamsTrue = {
    bands: ['R', 'G', 'B']
};
// Add both sets of NAIP imagery to the map to compare
coverage.
Map.addLayer(naip2015, visParamsTrue, '2015_true', false);
Map.addLayer(naip2017, visParamsTrue, '2017_true', false);
```
Compared to the 2015 imagery (Fig. 52.3), the 2017 imagery (Fig. 52.4) has been captured later in the season and we can start to see some yellow in the aspen. A challenge with this image is that there were some distinct time lags between neighboring flight paths during the initial image collection. Notice the vertical band of green in the 2017 image. While it is usually best practice to use images that closely match the time that you will be in the field, in this case, the consistency of the 2015 imagery outweighs those temporal concerns.

We will also add a false color image to the 2015 data because this band combination visualization is very helpful for distinguishing between deciduous and coniferous forests. We do this by adding the naip2015 object with a different set of visualization parameters to the map. Add the following code to your existing script.

```
// Add 2015 false color imagery.
Map.addLayer(naip2015, visParamsFalse, '2015_false',
false);
```
#### **<IP_ADDRESS> Section 1.3: Aspen Exclosures**

In our hypothetical study, let us imagine land managers have established some aspen exclosures on the southern extent of the Grand Mesa. Let us also say that the land managers did not have specific shapefiles of the exclosures but did have GPS locations of the four corners. We will use these data to add the exclosures

**Fig. 52.3** NAIP imagery from 2015 has captured the aspen forest across the extent of our study area

**Fig. 52.4** NAIP imagery from 2017 has captured the yellowing of the aspen forest but shows a distinct difference in vegetation cover due to variability in the time at which the given images were captured

to the map by creating a geometry feature within our script. In this case we are creating an ee.Geometry.MultiPolygon feature. Add the following code to your existing script. If your roi object does not encompass the exclosures, edit its bounds or redraw it.

```
// Creating a geometry feature.
var exclosure = ee.Geometry.MultiPolygon([
    [
        [-107.91079184, 39.012553345],
        [-107.90828129, 39.012553345],
        [-107.90828129, 39.014070552],
        [-107.91079184, 39.014070552],
        [-107.91079184, 39.012553345]
    ],
    [
        [-107.9512176, 39.00870162],
        [-107.9496834, 39.00870162],
        [-107.9496834, 39.00950196],
        [-107.95121765, 39.00950196],
        [-107.95121765, 39.00870162]
    ]
]);
print(exclosure);
Map.addLayer(exclosure, {}, 'exclosures');
```
The structure of the ee.Geometry.MultiPolygon feature is a bit complex, but it is effectively a set of nested lists. There are three tiers of lists present:

```
// Nested list example.
['tier 1' ['tier 2' ['tier 3']]]
```

**Fig. 52.5** View of the exclosure object coordinate details from the print statement

If you are trying to make a geometry feature and are having trouble, you can always create one with the draw shape tool and look at the values in a print statement as a template (Fig. 52.5).

#### **<IP_ADDRESS> Section 1.4: Determining Similar Areas for Sampling**

Now that we have our aspen exclosures loaded, we are going to bring in some additional layers to help quantify the landscape characteristics of the exclosures. We will use those values to find similar areas nearby to use as sampling sites outside of the exclosures. By keeping the environmental conditions similar, we can make stronger statements about the comparative effects of herbivory on aspen stands.

We will use three datasets to describe the conditions within the sampled areas:

#### **<IP_ADDRESS> Section 1.5: Loading in the Data**

First, we will call the NED by its unique ID. We can gather these details by searching for the feature and reading through the metadata. Add the following code to your existing script.

```
// Load in elevation dataset; clip it to general area.
var elev = ee.Image('USGS/NED')
    .clip(roi);
Map.addLayer(elev, {
    min: 1500,
    max: 3300
}, 'elevation', false);
```
We already have NAIP imagery loaded but we need to convert it to an image and calculate NDVI. We have already filtered our NAIP imagery to a single year, so there is only one image per area. We could apply a reducer to convert the ImageCollection to an image, but a reducer is not necessary for a single layer. A more logical option for this is to apply the mosaic function to convert the ImageCollection to an image. Add the following code to your existing script.

```
// Apply mosaic, clip, then calculate NDVI.
var ndvi = naip2015
    .mosaic()
    .clip(roi)
    .normalizedDifference(['N', 'R'])
    .rename('ndvi');
Map.addLayer(ndvi, {
    min: -0.8,
    max: 0.8
}, 'NDVI', false);
print(ndvi, 'ndvi');
```
The last dataset we are bringing in is the NLCD. Add the following code to your existing script.

```
// Add National Land Cover Database (NLCD).
var dataset = ee.ImageCollection('USGS/NLCD');
print(dataset, 'NLCD');
```
From the print statement (Fig. 52.6), we can see that some of the images have a band called 'landcover'*.* Rather than trying to pull it out from the feature collection, we will call the 2016 NLCD image directly by its unique ID and select the 'landcover' band.

**Fig. 52.6** Print statement showing the structure of the NLCD ImageCollection

```
// Load the selected NLCD image.
var landcover = ee.Image('USGS/NLCD/NLCD2016')
    .select('landcover')
    .clip(roi);
Map.addLayer(landcover, {}, 'Landcover', false);
```
**Question 1.** Investigate the other available years in the NLCD collection. What parts of the study area have land cover classes that change? What might be driving that change?

#### **<IP_ADDRESS> Section 1.6: Generating Random Points**

With our three datasets loaded, we can now generate a series of potential survey sites. We will do this by generating random points within a given area (Fig. 52.7). We want these sites to be accessible, near the two exclosures, and within the public land boundary. We will create another geometry feature that we will use to contain the randomly generated points. Hover over the **Geometry imports** box and click + **new layer**. Be sure to name this second geometry feature sampleArea.

With the geometry feature in place, we can create points using the randomPoints function. Add the following code to your existing script.

**Fig. 52.7** An example of what your second sample area geometry feature might look like

```
// Generate random points within the sample area.
var points = ee.FeatureCollection.randomPoints({
    region: sampleArea,
    points: 1000,
    seed: 1234
});
print(points, 'points');
Map.addLayer(points, {}, 'Points', false);
```
Here, we are using a dictionary to define the parameters of the function. The region parameter is the area in which the points are created. In our case, we are going to set this to sampleArea polygon. The points parameter specifies the total number of points to generate. The seed parameter is used to indicate specific random values. Think of this as a unique ID for a set of random values. The seed number (1234 in this example) refers to an existing random list of values. Setting the seed is very helpful because you are still using random values, but the code will place the points in the same locations every time until the seed is changed.

#### **<IP_ADDRESS> Section 1.7: Extracting Values to Points**

To associate the landscape characteristics with the point locations, we are going to call the ee.Image.sampleRegions function. This function requires a single image. Rather than calling three times on our three unique images (elevation, NDVI, and NLCD), we are going to add all those images together to create a multiband image so we can call this function a single time. Add the following code to your existing script.

```
// Add bands of elevation and NAIP.
var ndviElev = ndvi
    .addBands(elev)
    .addBands(landcover);
print(ndviElev, 'Multi band image');
```
With the multiband image set, we will call the sampleRegions function. Add the following code to your existing script.

```
// Extract values to points.
var samples = ndviElev.sampleRegions({
    collection: points,
    scale: 30,
    geometries: true
});
print(samples, 'samples');
Map.addLayer(samples, {}, 'samples', false);
```
Within this function, the collection parameter is set to the feature collection where the extracted values will be added. In this example, it is a point dataset. The scale parameter refers to the spatial resolution (pixel size) of the data. The geometries parameter indicates whether or not you want to maintain the x,y coordinate pairs for each element in the collection. The default is false, but we set it to true here because we eventually want to show these points on the map, as they will represent suitable sampling site locations.

#### **Segue on Scale**

The idea of *scale* is very important in remote sensing, but you will often encounter multiple definitions of it. *Map scale*, the relationship between a measured distance on a map and the actual distance on the landscape, is the most common in remote sensing. The image *extent* can also be referred to as the *scale of the image* or the *spatial scale*. Confusing, right?

In the case of the scale argument in the function above, it is referring to yet another definition of that word—the pixel size of the image. In our example, the multiband image has two bands with a pixel size of 30 m, and one with a pixel size of 1 m. It is best practice to use the largest pixel size when working with data having different spatial resolutions. Here, this means you are upscaling the 1 m image to 30 m. This means a potential loss of precision in your data. However, you can be confident that the number in that upscaled cell is representative of the mean value across all the cells. If you go in the opposite direction and downscale an image, you effectively make up the data to fill the gap. Because NAIP imagery has a 1 m pixel size, it contains 30× 30 more pixels than a single 30 m pixel. Take a look at the examples below to help visualize these processes.

Upscaling takes the available data and finds the mean value.

**Question 2.** Calculate the mean for the matrix above. How would this value change your understanding of the pixel? Is it still representative of the overall "story" the data is telling, or would a different reducer (i.e., median) be more appropriate?

**Fig. 52.8** Result of our sampling function has given us 1000 potential sites to choose from. We will limit this pool by comparing the measurable data we have at these sites to the mean values of those data from our exclosures

Downscaling takes a single value and places the value for all locations in the grid.

From our print statement (Fig. 52.8), we can see that each of our 1000 point locations has three properties: elevation, land cover, and NDVI. We want to use these values to filter out sites that do not match the conditions of the exclosures. The land cover data are categorical, so it is easy to filter. We know from the metadata on NLCD that the land cover class for deciduous forest is a single value (41).

We will use the filter function to select all sites that are within aspen forests. Add the following code to your existing script.

```
// Filter metadata for sites in the NLCD deciduous forest
layer.
var aspenSites =
samples.filter(ee.Filter.equals('landcover', 41));
print(aspenSites, 'Sites');
```
From the print statement, we can see this reduces the total number of potential sites (1000) by about 25%, depending on your study area shape. However, to get down to roughly 10 potential new monitoring sites, there is still a lot of trimming to do.

Filtering based on elevation and NDVI is a bit trickier because both of these variables are continuous data. You want to find sites that have similar values to those in the exclosures, but they do not need to be exactly the same. For this example, we will say that if a value is within±10% of the mean value found within the exclosure, we will group it as similar. There are a couple of things that need to be calculated before we can filter our potential sampling sites:

**Question 3.** The above use of±10% is somewhat arbitrary. How might the range change given the inclusion of additional exclosures in our study area? What other methods could be used to generate this range?

We will work with the NDVI image first and then apply this process to the elevation dataset.

• Step 1: To find the mean value, we are going to apply a mean reducer over the area that is inside of the exclosure. Add the following code to your existing script.

```
// Set the NDVI range.
var ndvi1 = ndvi
    .reduceRegion({
        reducer: ee.Reducer.mean(),
        geometry: exclosure,
        scale: 1,
        crs: 'EPSG:4326'
    });
print(ndvi1, 'Mean NDVI');
```
The reduceRegion function takes in an image and returns a dictionary where the key is the name of the band, and the value is the output of the reducer.

• Step 2: Determine the acceptable variability around the mean.

We are relying on simple mathematical functions to find the ± 10% values. Add the following code to your existing script.

```
// Generate a range of acceptable NDVI values.
var ndviNumber = ee.Number(ndvi1.get('ndvi'));
var ndviBuffer = ndviNumber.multiply(0.1);
var ndviRange = [
    ndviNumber.subtract(ndviBuffer),
    ndviNumber.add(ndviBuffer)
];
print(ndviRange, 'NDVI Range');
```
The first step in this process is all about understanding data types. The ndvi1 object is a dictionary so we call the get function to pull a value based on a known key. We then convert that value to an ee.Number type object so we can apply the math functions. Our output is a list with the minimum and maximum values (±10% of the mean NDVI values).

#### **<IP_ADDRESS> Section 1.8: Create Your Own Function**

We now have our range for NDVI, but we also need to apply the same process to elevation. In this case, we are only applying this process twice, but it seems like a useful chunk of code that might be handy down the road. Rather than just copying and pasting the code again and again, we are going to create a function with flexible parameters so we can apply this useful bit of code efficiently.

A function has the following requirements: parameters, statements, and a return value. Here is an example of the structure of a function from the official Earth Engine documentation. The following pseudocode demonstrates the syntax and structure of a function in JavaScript. *Do not add this code to your script.*

```
var myFunction = function(parameter1, parameter2,
parameter3) {
    statement;
    statement;
    var value = statement;
    return value;
};
```
If we apply this structure to our goal of reducing the elevation by area (like we did for NDVI with the code we created above), we would need to consider the following:

When using functions, it is important to use informative names within your parameters that give some indication of the data type that is required. If we want our function to be reproducible, we can provide some more information as a longer comment when we define the function. Add the following code to your existing script.

```
/*
This function is used to determine the mean value of an
image within a given area.
image: an image with a single band of ordinal or interval
level data
geom: geometry feature that overlaps with the image
pixelSize: a number that defines the cell size of the image
Returns a dictionary with the median values of the band,
the key is the band name.
*/
var reduceRegionFunction = function(image, geom, pixelSize)
{
    var dict = image.reduceRegion({
        reducer: ee.Reducer.mean(),
        geometry: geom,
        scale: pixelSize,
        crs: 'EPSG:4326'
    });
    return (dict);
};
```
Let us check our function definition by verifying that it gives the same answer for the NDVI range as when we did it step by step:

```
// Call function on the NDVI dataset to compare.
var ndvi_test = reduceRegionFunction(ndvi, exclosure, 1);
print(ndvi_test, 'ndvi_test');
```
This is a very clean method for coding when you need to apply a process multiple times. The function has been defined and now we can call it on the elevation dataset. Add the following code to your existing script.

```
// Call function on elevation dataset.
var elev1 = reduceRegionFunction(elev, exclosure, 30);
print(elev1, 'elev1');
```
We will define a second function to determine±10% around a mean value.

Add the following code to your existing script.

```
/*
Generate a range of acceptable values.
dictionary: a dictionary object
key: key to the value of interest, must be a string
proportion: a percentile to define the range of the values
around the mean
Returns a list with a min and max value for the given
range.
*/
var effectiveRange = function(dictionary, key, proportion)
{
    var number = ee.Number(dictionary.get(key));
    var buffer = number.multiply(proportion);
    var range = [
        number.subtract(buffer),
        number.add(buffer)
    ];
    return (range);
};
```
We will call the effectiveRange function on the output of the reduceRegionFunction function. Add the following code to your existing script.

```
// Call function on elevation data.
var elevRange = effectiveRange(elev1, 'elevation', 0.1);
print(elevRange);
```
Now that we have an effective range for both our NDVI and elevation values, we can apply an additional set of filters to thin the list of potential sample sites. We can do this by chaining multiple ee.Filter calls together. Add the following code to your existing script.

```
// Apply multiple filters to get at potential locations.
var combinedFilter = ee.Filter.and(
    ee.Filter.greaterThan('ndvi', ndviRange[0]),
    ee.Filter.lessThan('ndvi', ndviRange[1]),
    ee.Filter.greaterThan('elevation', elevRange[0]),
    ee.Filter.lessThan('elevation', elevRange[1])
);
var aspenSites2 = aspenSites.filter(combinedFilter);
print(aspenSites2, 'aspenSites2');
Map.addLayer(aspenSites2, {}, 'aspenSites2', false);
```
Depending on how you have drawn your study area, this filtering process should reduce the 1000 original sites by about 90%. From the ~ 100 remaining sites, we can manually select 10 either by something we already know about the landscape features of our study area, or by using a more restrictive range in our function (e.g.,±5%). This approach to site selection can be a good first step in ensuring you are sampling similar conditions in the field.

**Code Checkpoint A37a.** The book's repository contains a script that shows what your code should look like at this point.

#### **52.2.2 Section 2: Generating Your Own Training Dataset**