diff --git a/.DS_Store b/.DS_Store index bde7585..f015699 100644 Binary files a/.DS_Store and b/.DS_Store differ diff --git a/F1.qmd b/F1.qmd index af5a51c..6bd7358 100644 --- a/F1.qmd +++ b/F1.qmd @@ -1,31 +1,31 @@ -# Getting Started +## Getting Started -# Programming Basics +## Programming Basics ::: {.callout-tip} -# Chapter Information +## Chapter Information -## Author{.unlisted .unnumbered} +### Author{.unlisted .unnumbered} Ujaval Gandhi -## Overview {.unlisted .unnumbered} +### Overview {.unlisted .unnumbered} This chapter introduces the Google Earth Engine application programming interface (API) and the JavaScript syntax needed to use it. You will learn about the Code Editor environment and get comfortable typing, running, and saving scripts. You will also learn the basics of JavaScript language, such as variables, data structures, and functions. -## Learning Outcomes{.unlisted .unnumbered} +### Learning Outcomes{.unlisted .unnumbered} * Familiarity with the Earth Engine Code Editor. * Familiarity with the JavaScript syntax. * Ability to use the Earth Engine API functions from the Code Editor. -## Assumes you know how to:{.unlisted .unnumbered} +### Assumes you know how to:{.unlisted .unnumbered} * Sign up for an Earth Engine account (See the Google documentation for details). * Access the Earth Engine Code Editor (See the Google documentation for details). ::: -## Introduction {.unlisted .unnumbered} +### Introduction {.unlisted .unnumbered} In order to use Earth Engine well, you will need to develop basic skills in remote sensing and programming. The language of this book is JavaScript, and you will begin by learning how to manipulate variables using it. With that base, you’ll learn about viewing individual satellite images, viewing collections of images in Earth Engine, and how common remote sensing terms are referenced and used in Earth Engine. @@ -35,7 +35,7 @@ An API is a way to communicate with Earth Engine servers. It allows you to speci The Earth Engine platform comes with a web-based Code Editor that allows you to start using the Earth Engine JavaScript API without any installation. It also provides additional functionality to display your results on a map, save your scripts, access documentation, manage tasks, and more. It has a one-click mechanism to share your code with other users—allowing for easy reproducibility and collaboration. In addition, the JavaScript API comes with a user interface library, which allows you to create charts and web-based applications with little effort. -## Getting Started in the Code Editor +### Getting Started in the Code Editor If you have not already done so, be sure to add the book’s code repository to the Code Editor by entering [](https://www.google.com/url?q=https://code.earthengine.google.com/?accept_repo%3Dprojects/gee-edu/book&sa=D&source=editors&ust=1670414092101269&usg=AOvVaw2sJyDO_fhq1tcjG77pri7V)[https://code.earthengine.google.com/?accept_repo=projects/gee-edu/book](https://www.google.com/url?q=https://code.earthengine.google.com/?accept_repo%3Dprojects/gee-edu/book&sa=D&source=editors&ust=1670414092101852&usg=AOvVaw088kfXu4o_Mp4b8DJBPYjH) into your browser. The book’s scripts will then be available in the script manager panel. If you have trouble finding the repo, you can visit [this link](https://www.google.com/url?q=https://docs.google.com/presentation/d/1Kt6wGNoesYm__Cu3k3bnlbbyPN6m9SF4hQHK-pIDHfc/edit%23slide%3Did.g18a7b4b055d_0_624&sa=D&source=editors&ust=1670414092102526&usg=AOvVaw3ZCmkCOjrZEWqxfjRZPOCn) for help. @@ -76,11 +76,11 @@ Once the script is saved, it will appear in the script manager panel (Fig. F1.0. Now you should be familiar with how to create, run, and save your scripts in the Code Editor. You are ready to start learning the basics of JavaScript. -## JavaScript Basics +### JavaScript Basics To be able to construct a script for your analysis, you will need to use JavaScript. This section covers the JavaScript syntax and basic data structures. In the sections that follow, you will see more JavaScript code, noted in a distinct font and with shaded background. As you encounter code, paste it into the Code Editor and run the script. -### Variables {.unnumbered} +#### Variables {.unnumbered} In a programming language, variables are used to store data values. In JavaScript, a variable is defined using the var keyword followed by the name of the variable. The code below assigns the text “San Francisco” to the variable named city. Note that the text string in the code should be surrounded by quotes. You are free to use either ' (single quotes) or " (double quotes), and they must match at the beginning and end of each string. In your programs, it is advisable to be consistent—use either single quotes or double quotes throughout a given script (the code in this book generally uses single quotes for code). Each statement of your script should typically end with a semicolon, although Earth Engine’s code editor does not require it.  @@ -96,7 +96,7 @@ When you assign a text value, the variable is automatically assigned the type st var population = 873965; print(population); ``` -### Lists {.unnumbered} +#### Lists {.unnumbered} It is helpful to be able to store multiple values in a single variable. JavaScript provides a data structure called a list that can hold multiple values. We can create a new list using the square brackets [] and adding multiple values separated by a comma. ```js @@ -108,7 +108,7 @@ If you look at the output in the Console, you will see “List” with an expand ![Fig. F1.0.8 A JavaScript list](F1/image10.png) -### Objects {.unnumbered} +#### Objects {.unnumbered} Lists allow you to store multiple values in a single container variable. While useful, it is not appropriate to store structured data. It is helpful to be able to refer to each item with its name rather than its position. Objects in JavaScript allow you to store key-value pairs, where each value can be referred to by its key. You can create a dictionary using the curly braces {}. The code below creates an object called cityData with some information about San Francisco. @@ -125,7 +125,7 @@ We can use multiple lines to define the object. Only when we put in the semicolo ![Fig. F1.0.9 A JavaScript object](F1/image40.png) -### Functions {.unnumbered} +#### Functions {.unnumbered} While using Earth Engine, you will need to define your own functions. Functions take user inputs, use them to carry out some computation, and send an output back. Functions allow you to group a set of operations together and repeat the same operations with different parameters without having to rewrite them every time. Functions are defined using the `function` keyword. The code below defines a function called greet that takes an input called name and returns a greeting with Hello prefixed to it. Note that we can call the function with different input and it generates different outputs with the same code. ```js @@ -139,7 +139,7 @@ print(greet('Readers')); ``` ![Fig. F1.0.10 JavaScript function output](F1/image54.png) -### Comments {.unnumbbered} +#### Comments {.unnumbbered} While writing code, it is useful to add a bit of text to explain the code or leave a note for yourself. It is a good programming practice to always add comments in the code explaining each step. In JavaScript, you can prefix any line with two forward slashes // to make it a comment. The text in the comment will be ignored by the interpreter and will not be executed. @@ -153,7 +153,7 @@ Congratulations! You have learned enough JavaScript to be able to use the Earth Code Checkpoint F10a. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Earth Engine API Basics +### Earth Engine API Basics The Earth Engine API is vast and provides objects and methods to do everything from simple math to advanced algorithms for image processing. In the Code Editor, you can switch to the Docs tab to see the API functions grouped by object types. The API functions have the prefix ee (for Earth Engine). @@ -196,25 +196,25 @@ You just accomplished a moderately complex programming task with the help of Ear Code Checkpoint F10b. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Conclusion {.unnumbered} +### Conclusion {.unnumbered} This chapter introduced the Earth Engine API. You also learned the basics of JavaScript syntax to be able to use the API in the Code Editor environment. We hope you now feel a bit more comfortable starting your journey to become an Earth Engine developer. Regardless of your programming background or familiarity with JavaScript, you have the tools at your disposal to start using the Earth Engine API to build scripts for remote sensing analysis. -# Exploring Images +## Exploring Images ::: {.callout-tip} -# Chapter Information +## Chapter Information -## Author {.unlisted .unnumbered} +### Author {.unlisted .unnumbered} Jeff Howarth -## Overview {.unlisted .unnumbered} +### Overview {.unlisted .unnumbered} Satellite images are at the heart of Google Earth Engine’s power. This chapter teaches you how to inspect and visualize data stored in image bands. We first visualize individual bands as separate map layers and then explore a method to visualize three different bands in a single composite layer. We compare different kinds of composites for satellite bands that measure electromagnetic radiation in the visible and non-visible spectrum. We then explore images that represent more abstract attributes of locations, and create a composite layer to visualize change over time.   -## Learning Outcomes {.unlisted .unnumbered} +### Learning Outcomes {.unlisted .unnumbered} * Using the Code Editor to load an image * Using code to select image bands and visualize them as map layers @@ -222,12 +222,12 @@ Satellite images are at the heart of Google Earth Engine’s power. This chapter * Constructing new multiband images. * Understanding how additive color works and how to interpret RGB composites. -## Assumes you know how to: {.unlisted .unnumbered} +### Assumes you know how to: {.unlisted .unnumbered} * Sign up for an Earth Engine account, open the Code Editor, and save your script (Chap. F1.0). ::: -## Accessing an Image +### Accessing an Image If you have not already done so, be sure to add the book’s code repository to the Code Editor by entering [](https://www.google.com/url?q=https://code.earthengine.google.com/?accept_repo%3Dprojects/gee-edu/book&sa=D&source=editors&ust=1670414092189999&usg=AOvVaw1jWHeBmeq93I_lo_9useCA)[https://code.earthengine.google.com/?accept_repo=projects/gee-edu/book](https://www.google.com/url?q=https://code.earthengine.google.com/?accept_repo%3Dprojects/gee-edu/book&sa=D&source=editors&ust=1670414092190705&usg=AOvVaw3Z7cK8r6eOSYUceNjA8oUg) into your browser. The book’s scripts will then be available in the script manager panel. If you have trouble finding the repo, you can visit [this link](https://www.google.com/url?q=https://docs.google.com/presentation/d/1Kt6wGNoesYm__Cu3k3bnlbbyPN6m9SF4hQHK-pIDHfc/edit%23slide%3Did.g18a7b4b055d_0_624&sa=D&source=editors&ust=1670414092191415&usg=AOvVaw2eETuRpR5worezkj7citx6) for help. @@ -254,7 +254,7 @@ A satellite sensor like Landsat 5 measures the magnitude of radiation in differ An image band is an example of a raster data model, a method of storing geographic data in a two-dimensional grid of pixels, or picture elements. -## Visualizing an Image +### Visualizing an Image Now let’s add one of the bands to the map as a layer so that we can see it.   ```js @@ -324,7 +324,7 @@ By manipulating these controls, you should notice that these layers are displaye ::: {.callout-note} Code Checkpoint F11a. The book’s repository contains a script that shows what your code should look like at this point. ::: -## True-Color Composites +### True-Color Composites Using the controls in the Layers manager, explore these layers and examine how the pixel values in each band differ. Does Layer 2 (displaying pixel values from the “SR_B2” band) appear generally brighter than Layer 1 (the “SR_B1” band)? Compared with Layer 2, do the ocean waters in Layer 3 (the “SR_B3” band) appear a little darker in the north, but a little lighter in the south?   @@ -343,7 +343,7 @@ The result (Fig. F1.1.4) looks like the world we see, and is referred to as a na ![Fig. F1.1.4 True-color composite](F1/image39.png) -## False-Color Composites +### False-Color Composites As you saw when you printed the band list (Fig. F1.1.1), a Landsat image contains many more bands than just the three true-color bands. We can make RGB composites to show combinations of any of the bands—even those outside what the human eye can see. For example, band 4 represents the near-infrared band, just outside the range of human vision. Because of its value in distinguishing environmental conditions, this band was included on even the earliest 1970s Landsats. It has different values in coniferous and deciduous forests, for example, and can indicate crop health. To see an example of this, add this code to your script and run it.   ```js @@ -385,7 +385,7 @@ To compare the two false-color composites, zoom into the area shown in the two p Code Checkpoint F11b. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Attributes of Locations +### Attributes of Locations So far, we have explored bands as a method for storing data about slices of the electromagnetic spectrum that can be measured by satellites. Now we will work towards applying the additive color system to bands that store non-optical and more abstract attributes of geographic locations.   @@ -410,7 +410,7 @@ With the zoom controls on the map, you can zoom out to see the bright spot of S ![Fig. F1.1.10 Stable nighttime lights in 1993](F1/image34.png) -## Abstract RGB Composites   +### Abstract RGB Composites   Now we can use the additive color system to make an RGB composite that compares stable nighttime lights at three different slices of time. Add the code below to your script and run it.   @@ -460,32 +460,32 @@ As you explore this image, remember to check your interpretations with the Inspe Code Checkpoint F11c. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Conclusion {.unnumbered} +### Conclusion {.unnumbered} In this chapter, we looked at how an image is composed of one or more bands, where each band stores data about geographic locations as pixel values. We explored different ways of visualizing these pixel values as map layers, including a grayscale display of single bands and RGB composites of three bands. We created natural and false-color composites that use additive color to display information in visible and non-visible portions of the spectrum. We examined additive color as a general system for visualizing pixel values across multiple bands. We then explored how bands and RGB composites can be used to represent more abstract phenomena, including different kinds of change over time. -# Survey of Raster Datasets +## Survey of Raster Datasets The previous chapter introduced you to images, one of the core building blocks of remotely sensed imagery in Earth Engine. In this chapter, we will expand on this concept of images by introducing image collections. Image collections in Earth Engine organize many different images into one larger data storage structure. Image collections include information about the location, date collected, and other properties of each image, allowing you to sift through the ImageCollection for the exact image characteristics needed for your analysis. ::: {.callout-tip} -# Chapter Information +## Chapter Information -## Authors {.unlisted .unnumbered} +### Authors {.unlisted .unnumbered} Andréa Puzzi Nicolau, Karen Dyson, David Saah, Nicholas Clinton -## Overview {.unlisted .unnumbered} +### Overview {.unlisted .unnumbered} The purpose of this chapter is to introduce you to the many types of collections of images available in Google Earth Engine. These include sets of individual satellite images, pre-made composites (which merge multiple individual satellite images into one composite image), classified land use and land cover (LULC) maps, weather data, and other types of datasets. If you are new to JavaScript or programming, work through Chaps. F1.0 and F1.1 first.   -## Learning Outcomes {.unlisted .unnumbered} +### Learning Outcomes {.unlisted .unnumbered} * Accessing and viewing sets of images in Earth Engine. * Extracting single scenes from collections of images. * Applying visualization parameters in Earth Engine to visualize an image. -## Assumes you know how to: {.unlisted .unnumbered} +### Assumes you know how to: {.unlisted .unnumbered} * Sign up for an Earth Engine account, open the Code Editor, and save your script. (Chap. F1.0) * Locate the Earth Engine Inspector and Console tabs and understand their purposes (Chap. F1.0). @@ -493,13 +493,13 @@ The purpose of this chapter is to introduce you to the many types of collections ::: -## Image Collections: An Organized Set of Images +### Image Collections: An Organized Set of Images There are many different types of image collections available in Earth Engine. These include collections of individual satellite images, pre-made composites that combine multiple images into one blended image, classified LULC maps, weather data, and other non-optical data sets. Each one of these is useful for different types of analyses. For example, one recent study examined the drivers of wildfires in Australia (Sulova and Jokar 2021). The research team used the European Center for Medium-Range Weather Forecast Reanalysis (ERA5) dataset produced by the European Center for Medium-Range Weather Forecasts (ECMWF) and is freely available in Earth Engine. We will look at this dataset later in the chapter. You saw some of the basic ways to interact with an individual ee.Image in the previous chapter. However, depending on how long a remote sensing platform has been in operation, there may be thousands or millions of images collected of Earth. In Earth Engine, these are organized into an ImageCollection, a specialized data type that has specific operations available in the Earth Engine API. Like individual images, they can be viewed with Map.addLayer. -## View an Image Collection +### View an Image Collection The Landsat program from NASA and the United States Geological Survey (USGS) has launched a sequence of Earth observation satellites, named Landsat 1, 2, etc. Landsats have been returning images since 1972, making that collection of images the longest continuous satellite-based observation of the Earth's surface. We will now view images and basic information about one of the image collections that is still growing: collections of scenes taken by the Operational Land Imager aboard Landsat 8, which was launched in 2013. Copy and paste the following code into the center panel and click Run. While the enormous image catalog is accessed, it could take a couple of minutes to see the result in the Map area. If it takes more than a couple of minutes to see the images, try zooming in to a specific area to speed up the process. @@ -551,11 +551,11 @@ Code Checkpoint F12a. The book’s repository contains a script that shows what Edit your code to comment out the last two code commands you have written. This will remove the call to Map.addLayer that drew every image, and will remove the print statement that demanded more than 5000 elements. This will speed up your code in subsequent sections. Placing two forward slashes (//) at the beginning of a line will make it into a comment, and any commands on that line will not be executed. -## Filtering Image Collections +### Filtering Image Collections The ImageCollection data type in Earth Engine has multiple approaches to filtering, which helps to pinpoint the exact images you want to view or analyze from the larger collection. -### Filter by Date +#### Filter by Date One of the filters is filterDate, which allows us to narrow down the date range of the ImageCollection. Copy the following code to the center panel (paste it after the previous code you had): ```js @@ -581,7 +581,7 @@ Examine the mapped landsatWinter (Fig. F1.2.4). As described in the previous c Now look at the size of the winter Landsat 8 collection. The number is significantly lower than the number of images in the entire collection. This is the result of filtering the dates to three months in the winter of 2020–2021. -### Filter by Location +#### Filter by Location A second frequently used filtering tool is filterBounds. This filter is based on a location—for example, a point, polygon, or other geometry. Copy and paste the code below to filter and add to the map the winter images from the Landsat 8 Image Collection to a point in Minneapolis, Minnesota, USA. Note below the Map.addLayer function to add the pointMN to the map with an empty dictionary {} for the visParams argument. This only means that we are not specifying visualization parameters for this element, and it is being added to the map with the default parameters. ```js @@ -609,7 +609,7 @@ If we uncheck the Winter Landsat 8 layer under Layers, we can see that only imag The first still represents the map without zoom applied. The collection is shown inside the red circle. The second still represents the map after zoom was applied to the region. The red arrow indicates the point (in black) used to filter by bounds. -### Selecting the First Image +#### Selecting the First Image The final operation we will explore is the first function. This selects the first image in an ImageCollection. This allows us to place a single image on the screen for inspection. Copy and paste the code below to select and view the first image of the Minneapolis Winter Landsat 8 Image Collection. In this case, because the images are stored in time order in the ImageCollection, it will select the earliest image in the set. @@ -636,7 +636,7 @@ Code Checkpoint F12b. The book’s repository contains a script that shows what ::: Now that we have the tools to examine different image collections, we will explore other datasets. -## Collections of Single Images +### Collections of Single Images When learning about image collections in the previous section, you worked with the Landsat 8 raw image dataset. These raw images have some important corrections already done for you. However, the raw images are only one of several image collections produced for Landsat 8. The remote sensing community has developed additional imagery corrections that help increase the accuracy and consistency of analyses. The results of each of these different imagery processing paths is stored in a distinct ImageCollection in Earth Engine. @@ -688,7 +688,7 @@ Code Checkpoint F12c. The book’s repository contains a script that shows what ::: -## MODIS Monthly Burned Areas +### MODIS Monthly Burned Areas We’ll explore two examples of composites made with data from the MODIS sensors, a pair of sensors aboard the Terra and Aqua satellites. On these complex sensors, different MODIS bands produce data at different spatial resolutions. For the visible bands, the lowest common resolution is 500 m (red and NIR are 250 m). @@ -713,7 +713,7 @@ Code Checkpoint F12d. The book’s repository contains a script that shows what Save your script and start a new one by refreshing the page. -## Methane +### Methane Satellites can also collect information about the climate, weather, and various compounds present in the atmosphere. These satellites leverage portions of the electromagnetic spectrum and how different objects and compounds reflect when hit with sunlight in various wavelengths. For example, methane (CH4) reflects the 760 nm portion of the spectrum. Let’s take a closer look at a few of these datasets. @@ -750,7 +750,7 @@ Notice the different levels of methane over the African continent (Fig. F1.2.12) ![Fig. F1.2.12. Methane levels over the African continent on November 28, 2018](F1/image56.png) -## Global Forest Change +### Global Forest Change Another useful land cover product that has been pre-classified for you and is available in Earth Engine is the Global Forest Change dataset. This analysis was conducted between 2000 and 2020. Unlike the WorldCover dataset, this dataset focuses on the percent of tree cover across the Earth’s surface in a base year of 2000, and how that has changed over time. Copy and paste the code below to visualize the tree cover in 2000. Note that in the code below we define the visualization parameters as a variable treeCoverViz instead of having its calculation done within the Map.addLayer function. @@ -796,7 +796,7 @@ Code Checkpoint F12f. The book’s repository contains a script that shows what ::: Save your script and start a new one. -## Digital Elevation Models +### Digital Elevation Models Digital elevation models (DEMs) use airborne and satellite instruments to estimate the elevation of each location. Earth Engine has both local and global DEMs available. One of the global DEMs available is the NASADEM dataset, a DEM produced from a NASA mission. Copy and paste the code below to import the dataset and visualize the elevation band. @@ -821,11 +821,11 @@ Fig. F1.2.18. NASADEM elevation Code Checkpoint F12g. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Conclusion {.unnumbered} +### Conclusion {.unnumbered} In this chapter, we introduced image collections in Earth Engine and learned how to apply multiple types of filters to image collections to identify multiple or a single image for use. We also explored a few of the many different image collections available in the Earth Engine Data Catalog. Understanding how to find, access, and filter image collections is an important step in learning how to perform spatial analyses in Earth Engine. -## References {.unnumbered} +### References {.unnumbered} Chander G, Huang C, Yang L, et al (2009) Developing consistent Landsat data sets for large area applications: The MRLC 2001 protocol. IEEE Geosci Remote Sens Lett 6:777–781. https://doi.org/10.1109/LGRS.2009.2025244 @@ -835,36 +835,36 @@ Hansen MC, Potapov PV, Moore R, et al (2013) High-resolution global maps of 21st Sulova A, Arsanjani JJ (2021) Exploratory analysis of driving force of wildfires in Australia: An application of machine learning within Google Earth Engine. Remote Sens 13:1–23. https://doi.org/10.3390/rs13010010 -# The Remote Sensing Vocabulary +## The Remote Sensing Vocabulary ::: {.callout-tip} -# Chapter Information +## Chapter Information -## Authors {.unlisted .unnumbered} +### Authors {.unlisted .unnumbered} Karen Dyson, Andréa Puzzi Nicolau, David Saah, Nicholas Clinton -## Overview {.unlisted .unnumbered} +### Overview {.unlisted .unnumbered} The purpose of this chapter is to introduce some of the principal characteristics of remotely sensed images and how they can be examined in Earth Engine. We discuss spatial resolution, temporal resolution, and spectral resolution, along with how to access important image metadata. You will be introduced to image data from several sensors aboard various satellite platforms. At the completion of the chapter, you will be able to understand the difference between remotely sensed datasets based on these characteristics, and how to choose an appropriate dataset for your analysis based on these concepts.   -## Learning Outcomes {.unlisted .unnumbered} +### Learning Outcomes {.unlisted .unnumbered} * Understanding spatial, temporal, and spectral resolution. * Navigating the Earth Engine Console to gather information about a digital image, including resolution and other data documentation. -## Assumes you know how to: {.unlisted .unnumbered} +### Assumes you know how to: {.unlisted .unnumbered} * Navigate among Earth Engine result tabs (Chap. F1.0). * Visualize images with a variety of false-color band combinations (Chap. F1.1). ::: -## Introduction {.unlisted .unnumbered} +### Introduction {.unlisted .unnumbered} Images and image collections form the basis of many remote sensing analyses in Earth Engine. There are many different types of satellite imagery available to use in these analyses, but not every dataset is appropriate for every analysis. To choose the most appropriate dataset for your analysis, you should consider multiple factors. Among these are the resolution of the dataset—including the spatial, temporal, and spectral resolutions—as well as how the dataset was created and its quality. -## Searching for and Viewing Image Collection Information +### Searching for and Viewing Image Collection Information Earth Engine’s search bar can be used to find imagery and to locate important information about datasets in Earth Engine. Let’s use the search bar, located above the Earth Engine code, to find out information about the Landsat 7 Collection 2 Raw Scenes. First, type “landsat 7 collection 2” into the search bar (Fig. F1.3.1). Without hitting Enter, matches to that search term will appear. @@ -892,13 +892,13 @@ This more complete search results inset window contains short descriptions about Now that we know how to view this information, let’s dive into some important remote sensing terminology. -## Spatial Resolution +### Spatial Resolution Spatial resolution relates to the amount of Earth’s surface area covered by a single pixel. For example, we typically say that Landsat 7 has “30 m” color imagery. This means that each pixel is 30 m to a side, covering a total area of 900 square meters of the Earth’s surface. The spatial resolution of a given data set greatly affects the appearance of images, and the information in them, when you are viewing them on Earth’s surface. Next, we will visualize data from multiple sensors that capture data at different spatial resolutions, to compare the effect of different pixel sizes on the information and detail in an image. We’ll be selecting a single image from each ImageCollection to visualize. To view the image, we will draw them each as a color-IR image, a type of false-color image (described in detail in Chap. F1.1) that uses the infrared, red, and green bands. As you move through this portion of the course, zoom in and out to see differences in the pixel size and the image size. -### Landsat Thematic Mapper +#### Landsat Thematic Mapper Thematic Mapper (TM) sensors were flown aboard Landsat 4 and 5. TM data have been processed to a spatial resolution of 30m, and were active from 1982 to 2012. Search for “Landsat 5 TM” and import the result called “USGS Landsat 5 TM Collection 2 Tier 1 Raw Scenes”. In this dataset, the three bands for a color-IR image are called “B4” (infrared), “B3” (red), and “B2” (green). Let’s now visualize TM data over San Francisco airport. Note that we can either define the visualization parameters as a variable (as in the previous code snippet) or place them in curly braces in the Map.addLayer function (as in this code snippet). @@ -917,7 +917,7 @@ Map.addLayer(tmImage, { ``` ![Fig. F1.3.10 Visualizing the TM imagery from the Landsat 5 satellite](F1/image20.png) -### Sentinel-2 MultiSpectral Instrument +#### Sentinel-2 MultiSpectral Instrument The MultiSpectral Instrument (MSI) flies aboard the Sentinel-2 satellites, which are operated by the European Space Agency. The red, green, blue, and near-infrared bands are captured at 10m resolution, while other bands are captured at 20m and 30m. The Sentinel-2A satellite was launched in 2015 and the 2B satellite was launched in 2017. @@ -941,7 +941,7 @@ Compare the Sentinel imagery with the Landsat imagery, using the opacity slider ![Fig. F1.3.11 Visualizing the MSI imagery](F1/image1.png) -### National Agriculture Imagery Program (NAIP) +#### National Agriculture Imagery Program (NAIP) The National Agriculture Imagery Program (NAIP) is a U.S. government program to acquire imagery over the continental United States using airborne sensors. Data is collected for each state approximately every three years. The imagery has a spatial resolution of 0.5–2 m, depending on the state and the date collected.   @@ -970,11 +970,11 @@ Each of the datasets we’ve examined has a different spatial resolution. By com Code Checkpoint F13a. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Temporal Resolution +### Temporal Resolution Temporal resolution refers to the revisit time, or temporal cadence of a particular sensor’s image stream. Revisit time is the number of days between sequential visits of the satellite to the same location on the Earth’s surface. Think of this as the frequency of pixels in a time series at a given location. -### Landsat +#### Landsat The Landsat satellites 5 and later are able to image a given location every 16 days. Let’s use our existing tm dataset from Landsat 5. To see the time series of images at a location, you can filter an ImageCollection to an area and date range of interest and then print it. For example, to see the Landsat 5 images for three months in 1987, run the following code: @@ -1030,7 +1030,7 @@ When you print the chart, it will have a point each time an image was collected ![Fig. F1.3.15 A chart showing the temporal cadence, or temporal resolution of the Landsat 5 TM instrument at the San Francisco airport](F1/image47.png) -### Sentinel-2 +#### Sentinel-2 The Sentinel-2 program’s two satellites are in coordinated orbits, so that each spot on Earth gets visited about every 5 days. Within Earth Engine, images from these two sensors are pooled in the same dataset. Let’s create a chart using the MSI instrument dataset we have already imported. ```js @@ -1053,13 +1053,13 @@ Compare this Sentinel-2 graph (Fig. F1.3.16) with the Landsat graph you just pro ::: {.callout-note} Code Checkpoint F13b. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Spectral Resolution +### Spectral Resolution Spectral resolution refers to the number and width of spectral bands in which the sensor takes measurements. You can think of the width of spectral bands as the wavelength intervals for each band. A sensor that measures radiance in multiple bands is called a multispectral sensor (generally 3–10 bands), while a sensor with many bands (possibly hundreds) is called a hyperspectral sensor. Let’s compare the multispectral MODIS instrument with the hyperspectral Hyperion sensor aboard the EO-1 satellite, which is also available in Earth Engine. -### MODIS +#### MODIS There is an easy way to check the number of bands in an image: ```js @@ -1113,7 +1113,7 @@ The resulting chart is shown in Fig. F1.3.17. Use the expand button in the uppe ![Fig. F1.3.17 Plot of TOA reflectance for MODIS](F1/image50.png) -### EO-1 +#### EO-1 Now let’s compare MODIS with the EO-1 satellite’s hyperspectral sensor. Search for “eo-1” and import the “EO-1 Hyperion Hyperspectral Imager” dataset. Name it eo1. We can look at the number of bands from the EO-1 sensor. ```js @@ -1166,7 +1166,7 @@ Compare this hyperspectral instrument chart with the multispectral chart we plot ::: {.callout-note} Code Checkpoint F13c. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Per-Pixel Quality +### Per-Pixel Quality As you saw above, an image consists of many bands. Some of these bands contain spectral responses of Earth’s surface, including the NIR, red, and green bands we examined in the Spectral Resolution section. What about the other bands? Some of these other bands contain valuable information, like pixel-by-pixel quality-control data. @@ -1200,7 +1200,7 @@ Use the Inspector tool to examine some of the values. You may see values of 0 ( ::: {.callout-note} Code Checkpoint F13d. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Metadata +### Metadata In addition to band imagery and per-pixel quality flags, Earth Engine allows you to access substantial amounts of metadata associated with an image. This can all be easily printed to the Console for a single image. @@ -1228,10 +1228,10 @@ print('MSI CLOUDY_PIXEL_PERCENTAGE:', msiCloudiness); Code Checkpoint F13e. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Conclusion {.unnumbered} +### Conclusion {.unnumbered} A good understanding of the characteristics of your images is critical to your work in Earth Engine and the chapters going forward. You now know how to observe and query a variety of remote sensing datasets, and can choose among them for your work. For example, if you are interested in change detection, you might require a dataset with spectral resolution including near-infrared imagery and a fine temporal resolution. For analyses at a continental scale, you may prefer data with a coarse spatial scale, while analyses for specific forest stands may benefit from a very fine spatial scale. -## References {.unnumbered} +### References {.unnumbered} Fisher JRB, Acosta EA, Dennedy-Frank PJ, et al (2018) Impact of satellite imagery spatial resolution on land use classification accuracy and modeled water quality. Remote Sens Ecol Conserv 4:137–149. https://doi.org/10.1002/rse2.61 \ No newline at end of file diff --git a/F2.qmd b/F2.qmd index fc38c7d..13cc6be 100644 --- a/F2.qmd +++ b/F2.qmd @@ -1,38 +1,37 @@ -# Interpreting Images +## Interpreting Images Now that you know how images are viewed and what kinds of images exist in Earth Engine, how do we manipulate them? To gain the skills of interpreting images, you’ll work with bands, combining values to form indices and masking unwanted pixels. Then, you’ll learn some of the techniques available in Earth Engine for classifying images and interpreting the results. -# Image Manipulation: Bands, Arithmetic, Thresholds, and Masks +## Image Manipulation: Bands, Arithmetic, Thresholds, and Masks :::{.callout-tip} -# Chapter Information +## Chapter Information -## Author {.unlisted .unnumbered} +#### Author {.unlisted .unnumbered} Karen Dyson, Andréa Puzzi Nicolau, David Saah, and Nicholas Clinton -## Overview {.unlisted .unnumbered} +#### Overview {.unlisted .unnumbered} Once images have been identified in Earth Engine, they can be viewed in a wide array of band combinations for targeted purposes. For users who are already versed in remote sensing concepts, this chapter shows how to do familiar tasks on this platform; for those who are entirely new to such concepts, it introduces the idea of band combinations. -## Learning Outcomes {.unlisted .unnumbered} +#### Learning Outcomes {.unlisted .unnumbered} * Understanding what spectral indices are and why they are useful. * Being introduced to a range of example spectral indices used for a variety of purposes. -## Assumes you know how to:{.unlisted .unnumbered} +#### Assumes you know how to:{.unlisted .unnumbered} * Import images and image collections, filter, and visualize (Part F1). - - ::: -## Introduction {.unlisted .unnumbered} + +### Introduction {.unlisted .unnumbered} Spectral indices are based on the fact that different objects and land covers on the Earth’s surface reflect different amounts of light from the Sun at different wavelengths. In the visible part of the spectrum, for example, a healthy green plant reflects a large amount of green light while absorbing blue and red light—which is why it appears green to our eyes. Light also arrives from the Sun at wavelengths outside what the human eye can see, and there are large differences in reflectances between living and nonliving land covers, and between different types of vegetation, both in the visible and outside the visible wavelengths. We visualized this earlier, in Chaps. F1.1 and F1.3 when we mapped color-infrared images (Fig. F2.0.1). @@ -49,13 +48,13 @@ Spectral indices use math to express how objects reflect light across multiple p Indices derived from satellite imagery are used as the basis of many remote-sensing analyses. Indices have been used in thousands of applications, from detecting anthropogenic deforestation to examining crop health. For example, the growth of economically important crops such as wheat and cotton can be monitored throughout the growing season: Bare soil reflects more red wavelengths, whereas growing crops reflect more of the near-infrared (NIR) wavelengths. Thus, calculating a ratio of these two bands can help monitor how well crops are growing (Jackson and Huete 1991). -## Band Arithmetic in Earth Engine +### Band Arithmetic in Earth Engine If you have not already done so, be sure to add the book’s code repository to the Code Editor by entering [](https://www.google.com/url?q=https://code.earthengine.google.com/?accept_repo%3Dprojects/gee-edu/book&sa=D&source=editors&ust=1671458829783542&usg=AOvVaw2f8xfEZP6c0zP_Ke8jL26U)[https://code.earthengine.google.com/?accept_repo=projects/gee-edu/book](https://www.google.com/url?q=https://code.earthengine.google.com/?accept_repo%3Dprojects/gee-edu/book&sa=D&source=editors&ust=1671458829783919&usg=AOvVaw2i09J44MzpMZkjV_JLEnNR) into your browser. The book’s scripts will then be available in the script manager panel. If you have trouble finding the repo, you can visit [this link](https://www.google.com/url?q=https://docs.google.com/presentation/d/1Kt6wGNoesYm__Cu3k3bnlbbyPN6m9SF4hQHK-pIDHfc/edit%23slide%3Did.g18a7b4b055d_0_624&sa=D&source=editors&ust=1671458829784270&usg=AOvVaw1Kr82KG60ZeFLYC8cOZ67A) for help. Many indices can be calculated using band arithmetic in Earth Engine. Band arithmetic is the process of adding, subtracting, multiplying, or dividing two or more bands from an image. Here we’ll first do this manually, and then show you some more efficient ways to perform band arithmetic in Earth Engine. -### Arithmetic Calculation of NDVI +#### Arithmetic Calculation of NDVI The red and near-infrared bands provide a lot of information about vegetation due to vegetation’s high reflectance in these wavelengths. Take a look at Fig. F2.0.2 and note, in particular, that vegetation curves (graphed in green) have relatively high reflectance in the NIR range (approximately 750–900 nm). Also note that vegetation has low reflectance in the red range (approximately 630–690 nm), where sunlight is absorbed by chlorophyll. This suggests that if the red and near-infrared bands could be combined, they would provide substantial information about vegetation. @@ -124,7 +123,7 @@ Examine the resulting index, using the Inspector to pick out pixel values in are Using these simple arithmetic tools, you can build almost any index, or develop and visualize your own. Earth Engine allows you to quickly and easily calculate and display the index across a large area. -### Single-Operation Computation of Normalized Difference for NDVI +#### Single-Operation Computation of Normalized Difference for NDVI Normalized differences like NDVI are so common in remote sensing that Earth Engine provides the ability to do that particular sequence of subtraction, addition, and division in a single step, using the normalizedDifference method. This method takes an input image, along with bands you specify, and creates a normalized difference of those two bands. The NDVI computation previously created with band arithmetic can be replaced with one line of code: @@ -140,7 +139,7 @@ Map.addLayer(ndviND, { ``` Note that the order in which you provide the two bands to normalizedDifference is important. We use B8, the near-infrared band, as the first parameter, and the red band B4 as the second. If your two computations of NDVI do not look identical when drawn to the screen, check to make sure that the order you have for the NIR and red bands is correct. -### Using Normalized Difference for NDWI +#### Using Normalized Difference for NDWI As mentioned, the normalized difference approach is used for many different indices. Let’s apply the same normalizedDifference method to another index. @@ -171,11 +170,11 @@ Examine the areas of the map that NDVI identified as having a lot of vegetation. :::{.callout-note} Code Checkpoint F20a. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Thresholding, Masking, and Remapping Images +### Thresholding, Masking, and Remapping Images The previous section in this chapter discussed how to use band arithmetic to manipulate images. Those methods created new continuous values by combining bands within an image. This section uses logical operators to categorize band or index values to create a categorized image. -### Implementing a Threshold +#### Implementing a Threshold Implementing a threshold uses a number (the threshold value) and logical operators to help us partition the variability of images into categories. For example, recall our map of NDVI. High amounts of vegetation have NDVI values near 1 and non-vegetated areas are near 0. If we want to see what areas of the map have vegetation, we can use a threshold to generalize the NDVI value in each pixel as being either “no vegetation” or “vegetation”. That is a substantial simplification, to be sure, but can help us to better comprehend the rich variation on the Earth’s surface. This type of categorization may be useful if, for example, we want to look at the proportion of a city that is vegetated. Let’s create a Sentinel-2 map of NDVI near Seattle, Washington, USA. Enter the code below in a new script. @@ -229,7 +228,7 @@ Use the Inspector tool to explore this new layer. If you click on a green locati Other operators in this Boolean family include less than (lt), less than or equal to (lte), equal to (eq), not equal to (neq), and greater than or equal to (gte) and more. -### Building Complex Categorizations with .where +#### Building Complex Categorizations with .where A binary map classifying NDVI is very useful. However, there are situations where you may want to split your image into more than two bins. Earth Engine provides a tool, the where method, that conditionally evaluates to true or false within each pixel depending on the outcome of a test. This is analogous to an if statement seen commonly in other languages. However, to perform this logic when programming for Earth Engine, we avoid using the JavaScript if statement. Importantly, JavaScript if commands are not calculated on Google’s servers, and can create serious problems when running your code—in effect, the servers try to ship all of the information to be executed to your own computer’s browser, which is very underequipped for such enormous tasks. Instead, we use the where clause for conditional logic. @@ -260,7 +259,7 @@ There are a few interesting things to note about this code that you may not have ![Fig. F2.0.8 Thresholded water, forest, and non-forest image based on NDVI for Seattle, Washington, USA.](F2/image37.png) -### Masking Specific Values in an Image +#### Masking Specific Values in an Image Masking an image is a technique that removes specific areas of an image—those covered by the mask—from being displayed or analyzed. Earth Engine allows you to both view the current mask and update the mask. @@ -313,7 +312,7 @@ Map.addLayer(maskedVeg.mask(), {}, 'maskedVeg Mask'); ![Fig. F2.0.11 The updated mask. Areas of non-forest are now masked out as well (black areas of the image).](F2/image33.png) -### Remapping Values in an Image +#### Remapping Values in an Image Remapping takes specific values in an image and assigns them a different value. This is particularly useful for categorical datasets, including those you read about in Chap. F1.2 and those we have created earlier in this chapter. @@ -340,43 +339,11 @@ Use the inspector to compare values between our original seaWhere (displayed as :::{.callout-note} Code Checkpoint F20b. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Synthesis {.unnumbered} - -Assignment 1. In addition to vegetation indices and other land cover indices, you can use properties of different soil types to create geological indices. The Clay Minerals Ratio (CMR) is one of these. This index highlights soils containing clay and alunite, which absorb radiation in the SWIR portion (2.0–2.3 μm) of the spectrum. - -![](F2/image3.png) - -SWIR 1 should be in the 1.55–1.75 µm range, and SWIR 2 should be in the 2.08–2.35 µm range. Calculate and display CMR at the following point: ee.Geometry.Point(-100.543, 33.456). Don’t forget to use Map.centerObject. - -We’ve selected an area of Texas known for its clay soils. Compare this with an area without clay soils (for example, try an area around Seattle or Tacoma, Washington, USA). Note that this index will also pick up roads and other paved areas. - -Assignment 2. Calculate the Iron Oxide Ratio, which can be used to detect hydrothermally altered rocks (e.g., from volcanoes) that contain iron-bearing sulfides which have been oxidized (Segal, 1982). - -Here’s the formula: - -![](F2/image4.png) - -Red should be the 0.63–0.69 µm spectral range and Blue the 0.45–0.52 µm. Using Landsat 8, you can also find an interesting area to map by considering where these types of rocks might occur. - -Assignment 3. Calculate the Normalized Difference Built-Up Index (NDBI) for the sfoImage used in this chapter. - -The NDBI was developed by Zha et al. (2003) to aid in differentiating urban areas (e.g., densely clustered buildings and roads) from other land cover types. The index exploits the fact that urban areas, which generally have a great deal of impervious surface cover, reflect SWIR very strongly. If you like, refer back to Fig. F2.0.2. - -The formula is: - -![](F2/image5.png) - -Using what we know about Sentinel-2 bands, compute NDBI and display it. - -Bonus: Note that NDBI is the negative of NDWI computed earlier. We can prove this by using the JavaScript reverse method to reverse the palette used for NDWI in Earth Engine. This method reverses the order of items in the JavaScript list. Create a new palette for NDBI using the reverse method and display the map. As a hint, here is code to use the reverse method. - -var barePalette = waterPalette.reverse(); - -## Conclusion {.unnumbered} +### Conclusion {.unnumbered} In this chapter, you learned how to select multiple bands from an image and calculate indices. You also learned about thresholding values in an image, slicing them into multiple categories using thresholds. It is also possible to work with one set of class numbers and remap them quickly to another set. Using these techniques, you have some of the basic tools of image manipulation. In subsequent chapters you will encounter more complex and specialized image manipulation techniques, including pixel-based image transformations (Chap. F3.1), neighborhood-based image transformations (Chap. F3.2), and object-based image analysis (Chap. F3.3). -## References {.unnumbered} +### References {.unnumbered} Baig MHA, Zhang L, Shuai T, Tong Q (2014) Derivation of a tasselled cap transformation based on Landsat 8 at-satellite reflectance. Remote Sens Lett 5:423–431. https://doi.org/10.1080/2150704X.2014.915434 @@ -406,16 +373,16 @@ Souza Jr CM, Siqueira JV, Sales MH, et al (2013) Ten-year Landsat classification -# Interpreting an Image: Classification +## Interpreting an Image: Classification :::{.callout-tip} -# Chapter Information +## Chapter Information -## Author {.unlisted .unnumbered} +#### Author {.unlisted .unnumbered} @@ -423,12 +390,12 @@ Andréa Puzzi Nicolau, Karen Dyson, David Saah, Nicholas Clinton -## Overview {.unlisted .unnumbered} +#### Overview {.unlisted .unnumbered} Image classification is a fundamental goal of remote sensing. It takes the user from viewing an image to labeling its contents. This chapter introduces readers to the concept of classification and walks users through the many options for image classification in Earth Engine. You will explore the processes of training data collection, classifier selection, classifier training, and image classification. -## Learning Outcomes {.unlisted .unnumbered} +#### Learning Outcomes {.unlisted .unnumbered} * Running a classification in Earth Engine. @@ -437,14 +404,14 @@ Image classification is a fundamental goal of remote sensing. It takes the user * Learning how to collect sample data in Earth Engine. * Learning the basics of the hexadecimal numbering system. -## Assumes you know how to:{.unlisted .unnumbered} +#### Assumes you know how to:{.unlisted .unnumbered} * Import images and image collections, filter, and visualize (Part F1). * Understand bands and how to select them (Chap. F1.2, Chap. F2.0). ::: -## Introduction {.unlisted .unnumbered} +### Introduction {.unlisted .unnumbered} Classification is addressed in a broad range of fields, including mathematics, statistics, data mining, machine learning, and more. For a deeper treatment of classification, interested readers may see some of the following suggestions: Witten et al. (2011), Hastie et al. (2009), Goodfellow et al. (2016), Gareth et al. (2013), Géron (2019), Müller et al. (2016), or Witten et al. (2005). Unlike regression, which predicts continuous variables, classification predicts categorical, or discrete, variables—variables with a finite number of categories (e.g., age range). @@ -459,7 +426,7 @@ Image classification techniques for generating land cover and land use informati It is important to define land use and land cover. Land cover relates to the physical characteristics of the surface: simply put, it documents whether an area of the Earth’s surface is covered by forests, water, impervious surfaces, etc. Land use refers to how this land is being used by people. For example, herbaceous vegetation is considered a land cover but can indicate different land uses: the grass in a pasture is an agricultural land use, whereas the grass in an urban area can be classified as a park. -## Supervised Classification +### Supervised Classification If you have not already done so, be sure to add the book’s code repository to the Code Editor by entering [](https://www.google.com/url?q=https://code.earthengine.google.com/?accept_repo%3Dprojects/gee-edu/book&sa=D&source=editors&ust=1671458829866098&usg=AOvVaw16x5swm9HlorS5Mbw7E42X)[https://code.earthengine.google.com/?accept_repo=projects/gee-edu/book](https://www.google.com/url?q=https://code.earthengine.google.com/?accept_repo%3Dprojects/gee-edu/book&sa=D&source=editors&ust=1671458829866485&usg=AOvVaw0-N-JCWWgnM493BKa7Ichm) into your browser. The book’s scripts will then be available in the script manager panel. If you have trouble finding the repo, you can visit [this link](https://www.google.com/url?q=https://docs.google.com/presentation/d/1Kt6wGNoesYm__Cu3k3bnlbbyPN6m9SF4hQHK-pIDHfc/edit%23slide%3Did.g18a7b4b055d_0_624&sa=D&source=editors&ust=1671458829866823&usg=AOvVaw0ytMyRvutssBcVr2GdcBHA) for help. @@ -679,7 +646,7 @@ Inspect the result (Fig. F2.1.13). How does this classified image differ from th :::{.callout-note} Code Checkpoint F21b. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Unsupervised Classification +### Unsupervised Classification In an unsupervised classification, we have the opposite process of supervised classification. Spectral classes are grouped first and then categorized into clusters. Therefore, in Earth Engine, these classifiers are ee.Clusterer objects. They are “self-taught” algorithms that do not use a set of labeled training data (i.e., they are “unsupervised”). You can think of it as performing a task that you have not experienced before, starting by gathering as much information as possible. For example, imagine learning a new language without knowing the basic grammar, learning only by watching a TV series in that language, listening to examples, and finding patterns. @@ -737,27 +704,11 @@ Another key point of classification is the accuracy assessment of the results. T :::{.callout-note} Code Checkpoint F21c. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Synthesis {.unnumbered} - -Test if you can improve the classifications by completing the following assignments. - -Assignment 1. For the supervised classification, try collecting more points for each class. The more points you have, the more spectrally represented the classes are. It is good practice to collect points across the entire composite and not just focus on one location. Also look for pixels of the same class that show variability. For example, for the water class, collect pixels in parts of rivers that vary in color. For the developed class, collect pixels from different rooftops. - -Assignment 2. Add more predictors. Usually, the more spectral information you feed the classifier, the easier it is to separate classes. Try calculating and incorporating a band of NDVI or the Normalized Difference Water Index (Chap. F2.0) as a predictor band. Does this help the classification? Check for developed areas that were being classified as herbaceous or vice versa. - -Assignment 3. Use more trees in the Random Forest classifier. Do you see any improvements compared to 50 trees? Note that the more trees you have, the longer it will take to compute the results, and that more trees might not always mean better results. - -Assignment 4. Increase the number of samples that are extracted from the composite in the unsupervised classification. Does that improve the result? - -Assignment 5. Increase the number k of clusters for the k-means algorithm. What would happen if you tried 10 classes? Does the classified map result in meaningful classes? - -Assignment 6. Test other clustering algorithms. We only used k-means; try other options under the ee.Clusterer object. - -## Conclusion {.unnumbered} +### Conclusion {.unnumbered} Classification algorithms are key for many different applications because they allow you to predict categorical variables. You should now understand the difference between supervised and unsupervised classification and have the basic knowledge on how to handle misclassifications. By being able to map the landscape for land use and land cover, we will also be able to monitor how it changes (Part F4). -## References {.unnumbered} +### References {.unnumbered} Breiman L (2001) Random forests. Mach Learn 45:5–32. https://doi.org/10.1023/A:1010933404324 @@ -779,16 +730,16 @@ Witten IH, Frank E, Hall MA, et al (2005) Practical machine learning tools and t -# Accuracy Assessment: Quantifying Classification Quality +## Accuracy Assessment: Quantifying Classification Quality :::{.callout-tip} -# Chapter Information +## Chapter Information -## Author {.unlisted .unnumbered} +#### Author {.unlisted .unnumbered} @@ -796,12 +747,12 @@ Andréa Puzzi Nicolau, Karen Dyson, David Saah, Nicholas Clinton -## Overview {.unlisted .unnumbered} +#### Overview {.unlisted .unnumbered} This chapter will enable you to assess the accuracy of an image classification. You will learn about different metrics and ways to quantify classification quality in Earth Engine. Upon completion, you should be able to evaluate whether your classification needs improvement and know how to proceed when it does. -## Learning Outcomes {.unlisted .unnumbered} +#### Learning Outcomes {.unlisted .unnumbered} * Learning how to perform accuracy assessment in Earth Engine. @@ -809,14 +760,14 @@ This chapter will enable you to assess the accuracy of an image classification. * Understanding overall accuracy and the kappa coefficient. * Understanding the difference between user’s and producer’s accuracy, and the difference between omission and commission errors. -## Assumes you know how to:{.unlisted .unnumbered} +#### Assumes you know how to:{.unlisted .unnumbered} * ​​Create a graph using ui.Chart (Chap. F1.3). * Perform a supervised Random Forest image classification (Chap. F2.1). ::: -## Introduction {.unlisted .unnumbered} +### Introduction {.unlisted .unnumbered} Any map or remotely sensed product is a generalization or model that will have inherent errors. Products derived from remotely sensed data used for scientific purposes and policymaking require a quantitative measure of accuracy to strengthen the confidence in the information generated (Foody 2002, Strahler et al. 2006, Olofsson et al. 2014). Accuracy assessment is a crucial part of any classification project, as it measures the degree to which the classification agrees with another data source that is considered to be accurate, ground-truth data (i.e., “reality”). @@ -828,7 +779,7 @@ In Chap. F2.1, we asked whether the classification results were satisfactory. In In a thorough accuracy assessment, we think carefully about the sampling design, the response design, and the analysis (Olofsson et al. 2014). Fundamental protocols are taken into account to produce scientifically rigorous and transparent estimates of accuracy and area, which requires robust planning and time. In a standard setting, we would calculate the number of samples needed for measuring accuracy (sampling design). Here, we will focus mainly on the last step, analysis, by examining the confusion matrix and learning how to calculate the accuracy metrics. This will be done by partitioning the existing data into training and testing sets. -## Quantifying Classification Accuracy Through a Confusion Matrix +### Quantifying Classification Accuracy Through a Confusion Matrix If you have not already done so, be sure to add the book’s code repository to the Code Editor by entering [](https://www.google.com/url?q=https://code.earthengine.google.com/?accept_repo%3Dprojects/gee-edu/book&sa=D&source=editors&ust=1671458829937499&usg=AOvVaw3qqOwSX_A-Pllh6X3X31q4)[https://code.earthengine.google.com/?accept_repo=projects/gee-edu/book](https://www.google.com/url?q=https://code.earthengine.google.com/?accept_repo%3Dprojects/gee-edu/book&sa=D&source=editors&ust=1671458829937976&usg=AOvVaw0WioXIhzue8-WoaX4UtabH) into your browser. The book’s scripts will then be available in the script manager panel. If you have trouble finding the repo, you can visit [this link](https://www.google.com/url?q=https://docs.google.com/presentation/d/1Kt6wGNoesYm__Cu3k3bnlbbyPN6m9SF4hQHK-pIDHfc/edit%23slide%3Did.g18a7b4b055d_0_624&sa=D&source=editors&ust=1671458829938470&usg=AOvVaw2CH8V3-_qV99EcgMxUAaSO) for help. @@ -1018,7 +969,7 @@ How is the classification accuracy? Which classes have higher accuracy compared :::{.callout-note} Code Checkpoint F22a. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Hyperparameter tuning +### Hyperparameter tuning We can also assess how the number of trees in the Random Forest classifier affects the classification accuracy. Copy and paste the code below to create a function that charts the overall accuracy versus the number of trees used. The code tests from 5 to 100 trees at increments of 5, producing Fig. F2.2.2. (Do not worry too much about fully understanding each item at this stage of your learning. If you want to find out how these operations work, you can see more in Chaps. F4.0 and F4.1.) @@ -1063,13 +1014,7 @@ We might also want to ensure that the samples from the training set are uncorrel :::{.callout-note} Code Checkpoint F22c. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Synthesis {.unnumbered} - -Assignment 1. Based on Sect. 1, test other classifiers (e.g., a Classification and Regression Tree or Support Vector Machine classifier) and compare the accuracy results with the Random Forest results. Which model performs better? - -Assignment 2. Try setting a different seed in the randomColumn method and see how that affects the accuracy results. You can also change the split between the training and testing sets (e.g., 70/30 or 60/40). - -## Conclusion {.unnumbered} +### Conclusion {.unnumbered} You should now understand how to calculate how well your classifier is performing on the data used to build the model. This is a useful way to understand how a classifier is performing, because it can help indicate which classes are performing better than others. A poorly modeled class can sometimes be improved by, for example, collecting more training points for that class. diff --git a/F4.qmd b/F4.qmd index cf10e7c..0ce5635 100644 --- a/F4.qmd +++ b/F4.qmd @@ -1,25 +1,25 @@ -# Interpreting Image Series +## Interpreting Image Series One of the paradigm-changing features of Earth Engine is the ability to access decades of imagery without the previous limitation of needing to download all the data to a local disk for processing. Because remote-sensing data files can be enormous, this used to limit many projects to viewing two or three images from different periods. With Earth Engine, users can access tens or hundreds of thousands of images to understand the status of places across decades. -# Filter, Map, Reduce +## Filter, Map, Reduce :::{.callout-tip} -# Chapter Information +## Chapter Information -## Author {.unlisted .unnumbered} +#### Author {.unlisted .unnumbered} Jeffrey A. Cardille -## Overview {.unlisted .unnumbered} +#### Overview {.unlisted .unnumbered} The purpose of this chapter is to teach you important programming concepts as they are applied in Earth Engine. We first illustrate how the order and type of these operations can matter with a real-world, non-programming example. We then demonstrate these concepts with an ImageCollection, a key data type that distinguishes Earth Engine from desktop image-processing implementations. -## Learning Outcomes {.unlisted .unnumbered} +#### Learning Outcomes {.unlisted .unnumbered} * Visualizing the concepts of filtering, mapping, and reducing with a hypothetical, non-programming example. @@ -27,14 +27,14 @@ The purpose of this chapter is to teach you important programming concepts as th * Learning how to efficiently map a user-written function over the images of a filtered ImageCollection. * Learning how to summarize a set of assembled values using Earth Engine reducers. -## Assumes you know how to:{.unlisted .unnumbered} +#### Assumes you know how to:{.unlisted .unnumbered} * Import images and image collections, filter, and visualize (Part F1). * Perform basic image analysis: select bands, compute indices, create masks (Part F2). ::: -## Introduction {.unlisted .unnumbered} +### Introduction {.unlisted .unnumbered} Prior chapters focused on exploring individual images—for example, viewing the characteristics of single satellite images by displaying different combinations of bands (Chap. F1.1), viewing single images from different datasets (Chap. F1.2, Chap. F1.3), and exploring image processing principles (Parts F2, F3) as they are implemented for cloud-based remote sensing in Earth Engine. Each image encountered in those chapters was pulled from a larger assemblage of images taken from the same sensor. The chapters used a few ways to narrow down the number of images in order to view just one for inspection (Part F1) or manipulation (Part F2, Part F3). @@ -68,7 +68,7 @@ In addition to the greater order of magnitude of the list size, you can see that This counterexample would also have other complications—such as the need to track store locations on the list of milk prices—that could present serious problems if you did those operations in that unwise order. For now, the point is to filter, then map, then reduce. Below, we’ll apply these concepts to image collections. -## Filtering Image Collections in Earth Engine +### Filtering Image Collections in Earth Engine The first part of the filter, map, reduce paradigm is “filtering” to get a smaller ImageCollection from a larger one. As in the milk example, filters take a large set of items, limit it by some criterion, and return a smaller set for consideration. Here, filters take an ImageCollection, limit it by some criterion of date, location, or image characteristics, and return a smaller ImageCollection (Fig. F4.0.1). @@ -145,7 +145,7 @@ Code Checkpoint F40a. The book’s repository contains a script that shows what ::: Now, with an efficiently filtered collection that satisfies our chosen criteria, we will next explore the second stage: executing a function for all of the images in the set. -## Mapping over Image Collections in Earth Engine +### Mapping over Image Collections in Earth Engine In Chap. F3.1, we calculated the Enhanced Vegetation Index (EVI) in very small steps to illustrate band arithmetic on satellite images. In that chapter, code was called once, on a single image. What if we wanted to compute the EVI in the same way for every image of an entire ImageCollection? Here, we use the key tool for the second part of the workflow in Earth Engine, a .map command (Fig. F4.0.1). This is roughly analogous to the step of making phone calls in the milk example that began this chapter, in which you took a list of store names and transformed it through effort into a list of milk prices. @@ -175,7 +175,7 @@ After entering and executing this code, you will see a grayscale image. If you l :::{.callout-note} Code Checkpoint F40b. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Reducing an Image Collection +### Reducing an Image Collection The third part of the filter, map, reduce paradigm is “reducing” values in an ImageCollection to extract meaningful values (Fig. F4.0.1). In the milk example, we reduced a large list of milk prices to find the minimum value. The Earth Engine API provides a large set of reducers for reducing a set of values to a summary statistic. @@ -213,36 +213,22 @@ The reducers described here treat each pixel independently. In subsequent chapte :::{.callout-note} Code Checkpoint F40c. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Synthesis {.unnumbered} - -Assignment 1. Compare the mean and median images produced in Sect. 3 (Fig. 4.0.2). In what ways do they look different, and in what ways do they look alike? To understand how they work, pick a pixel and inspect the EVI values computed. In your opinion, which is a better representative of the data set? - -Assignment 2. Adjust the filters to filter a different proportion of clouds, or a different date range. What effects do these changes have on the number of images and the look of the reductions made from them? - -Assignment 3. Explore the ee.Filter options in the API documentation, and select a different filter that might be of interest. Filter images using it, and comment on the number of images and the reductions made from them. - -Assignment 4. Change the EVI function so that it returns the original image with the EVI band appended by replacing the return statement with this: return (oneL5Image.addBands(landsat5EVI)) - -What does the median reducer return in that case? Some EVI values are 0. What are the conditions in which this occurs? - -Assignment 5. Choose a date and location that is important to you (e.g., your birthday and your place of birth). Filter Landsat imagery to get all the low-cloud imagery at your location within 6 months of the date. Then, reduce the ImageCollection to find the median EVI. Describe the image and how representative of the full range of values it is, in your opinion. - -## Conclusion {.unnumbered} +### Conclusion {.unnumbered} In this chapter, you learned about the paradigm of filter, map, reduce. You learned how to use these tools to sift through, operate on, and summarize a large set of images to suit your purposes. Using the Filter functionality, you learned how to take a large ImageCollection and filter away images that do not meet your criteria, retaining only those images that match a given set of characteristics. Using the Map functionality, you learned how to apply a function to each image in an ImageCollection, treating each image one at a time and executing a requested set of operations on each. Using the Reduce functionality, you learned how to summarize the elements of an ImageCollection, extracting summary values of interest. In the subsequent chapters of Part 4, you will encounter these concepts repeatedly, manipulating image collections according to your project needs using the building blocks seen here. By building on what you have done in this chapter, you will grow in your ability to do sophisticated projects in Earth Engine. -# Exploring Image Collections +## Exploring Image Collections :::{.callout-tip} -# Chapter Information +## Chapter Information -## Author {.unlisted .unnumbered} +#### Author {.unlisted .unnumbered} @@ -250,12 +236,12 @@ Gennadii Donchyts -## Overview {.unlisted .unnumbered} +#### Overview {.unlisted .unnumbered} This chapter teaches how to explore image collections, including their spatiotemporal extent, resolution, and values stored in images and image properties. You will learn how to map and inspect image collections using maps, charts, and interactive tools, and how to compute different statistics of values stored in image collections using reducers. -## Learning Outcomes {.unlisted .unnumbered} +#### Learning Outcomes {.unlisted .unnumbered} * Inspecting the spatiotemporal extent and resolution of image collections by mapping image geometry and plotting image time properties. @@ -263,7 +249,7 @@ This chapter teaches how to explore image collections, including their spatiotem * Filtering image collections by using stored or computed image properties. * Exploring the distribution of values stored in image pixels of an ImageCollection through percentile reducers. -## Assumes you know how to:{.unlisted .unnumbered} +#### Assumes you know how to:{.unlisted .unnumbered} * Import images and image collections, filter, and visualize (Part F1). @@ -274,7 +260,7 @@ This chapter teaches how to explore image collections, including their spatiotem In the previous chapter (Chap. F4.0), the filter, map, reduce paradigm was introduced. The main goal of this chapter is to demonstrate some of the ways that those concepts can be used within Earth Engine to better understand the variability of values stored in image collections. Sect. 1 demonstrates how time-dependent values stored in the images of an ImageCollection can be inspected using the Code Editor user interface after filtering them to a limited spatiotemporal range (i.e., geometry and time ranges). Sect. 2 shows how the extent of images, as well as basic statistics, such as the number of observations, can be visualized to better understand the spatiotemporal extent of image collections. Then, Sects. 3 and 4 demonstrate how simple reducers such as mean and median, and more advanced reducers such as percentiles, can be used to better understand how the values of a filtered ImageCollection are distributed. -## Filtering and Inspecting an Image Collection +### Filtering and Inspecting an Image Collection We will focus on the area in and surrounding Lisbon, Portugal. Below, we will define a point, lisbonPoint, located in the city; access the very large Landsat ImageCollection and limit it to the year 2020 and to the images that contain Lisbon; and select bands 6, 5, and 4 from each of the images in the resulting filtered ImageCollection. @@ -320,7 +306,7 @@ print(chart); :::{.callout-note} Code Checkpoint F41a. The book’s repository contains a script that shows what your code should look like at this point. ::: -## How Many Images Are There, Everywhere on Earth? +### How Many Images Are There, Everywhere on Earth? Suppose we are interested to find out how many valid observations we have at every map pixel on Earth for a given ImageCollection. This enormously computationally demanding task is surprisingly easy to do in Earth Engine. The API provides a set of reducer functions to summarize values to a single number in each pixel, as described in Chap. F4.0. We can apply this reducer, count, to our filtered ImageCollection with the code below. We’ll return to the same data set and filter for 2020, but without the geographic limitation. This will assemble images from all over the world, and then count the number of images in each pixel. The following code does that count, and adds the resulting image to the map with a predefined red/yellow/green color palette stretched between values 0 and 50. Continue pasting the code below into the same script. @@ -360,7 +346,7 @@ You might have noticed that we summarized a single band from the original ImageC :::{.callout-note} Code Checkpoint F41b. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Reducing Image Collections to Understand Band Values +### Reducing Image Collections to Understand Band Values As we have seen, you could click at any point on Earth’s surface and see both the number of Landsat images recorded there in 2020 and the values of any image in any band through time. This is impressive and perhaps mind-bending, given the enormous amount of data in play. In this section and the next, we will explore two ways to summarize the numerical values of the bands—one straightforward way and one more complex but highly powerful way to see what information is contained in image collections. @@ -388,7 +374,7 @@ There is a wide range of reducers available in Earth Engine. If you are curious :::{.callout-note} Code Checkpoint F41c. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Compute Multiple Percentile Images for an Image Collection +### Compute Multiple Percentile Images for an Image Collection One particularly useful reducer that can help you better understand the variability of values in image collections is ee.Reducer.percentile. The nth percentile gives the value that is the nth largest in a set. In this context, you can imagine accessing all of the values for a given band in a given ImageCollection for a given pixel and sorting them. The 30th percentile, for example, is the value 30% of the way along the list from smallest to largest. This provides an easy way to explore the variability of the values in image collections by computing a cumulative density function of values on a per-pixel basis. The following code shows how we can calculate a single 30th percentile on a per-pixel and per-band basis for our Landsat 8 ImageCollection. Continue pasting the code below into the same script. @@ -428,36 +414,26 @@ Earth Engine provides a very rich API, allowing users to explore image collectio :::{.callout-note} Code Checkpoint F41d. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Synthesis {.unnumbered} - -In the example above, the 30th percentile composite image would be useful for typical studies that need cloud-free data for analysis. The “best” composite to use, however, will depend on the goal of a study, the characteristics of the given data set, and the location being viewed. You can imagine choosing different percentile composite values if exploring image collections over the Sahara Desert or over Congo, where cloud frequency would vary substantially (Wilson et al. 2016). - -Assignment 1. Noting that your own interpretation of what constitutes a good composite is subjective, create a series of composites of a different location, or perhaps a pair of locations, for a given set of dates. - -Assignment 2. Filter to create a relevant data set—for example, for Landsat 8 or Sentinel-2 over an agricultural growing season. Create percentile composites for a given location. Which image composite is the most satisfying, and what type of project do you have in mind when giving that response? - -Assignment 3. Do you think it is possible to generalize about the relationship between the time window of an ImageCollection and the percentile value that will be the most useful for a given project, or will every region need to be inspected separately? - -## Conclusion {.unnumbered} +### Conclusion {.unnumbered} In this chapter, you have learned different ways to explore image collections using Earth Engine in addition to looking at individual images. You have learned that image collections in Earth Engine may have global footprints as well as images with a smaller, local footprint, and how to visualize the number of images in a given filtered ImageCollection. You have learned how to explore the temporal and spatial extent of images stored in image collections, and how to quickly examine the variability of values in these image collections by computing simple statistics like mean or median, as well as how to use a percentile reducer to better understand this variability. -## References {.unnumbered} +### References {.unnumbered} Wilson AM, Jetz W (2016) Remotely sensed high-resolution global cloud dynamics for predicting ecosystem and biodiversity distributions. PLoS Biol 14:e1002415. https://doi.org/10.1371/journal.pbio.1002415 -# Aggregating Images for Time Series +## Aggregating Images for Time Series :::{.callout-tip} -# Chapter Information +## Chapter Information -## Author {.unlisted .unnumbered} +#### Author {.unlisted .unnumbered} @@ -465,7 +441,7 @@ Ujaval Gandhi -## Overview {.unlisted .unnumbered} +#### Overview {.unlisted .unnumbered} Many remote sensing datasets consist of repeated observations over time. The interval between observations can vary widely. The Global Precipitation Measurement dataset, for example, produces observations of rain and snow worldwide every three hours. The Climate Hazards Group InfraRed Precipitation with Station (CHIRPS) project produces a gridded global dataset at the daily level and also for each five-day period. The Landsat 8 mission produces a new scene of each location on Earth every 16 days. With its constellation of two satellites, the Sentinel-2 mission images every location every five days. @@ -476,14 +452,14 @@ While individual scenes are informative, many days are cloudy, and it is useful This chapter will cover the techniques for aggregating individual images from a time series at a chosen interval. We will take the CHIRPS time series of rainfall for one year and aggregate it to create a monthly rainfall time series. -## Learning Outcomes {.unlisted .unnumbered} +#### Learning Outcomes {.unlisted .unnumbered} * Using the Earth Engine API to work with dates. * Aggregating values from an ImageCollection to calculate monthly, seasonal, or yearly images. * Plotting the aggregated time series at a given location. -## Assumes you know how to:{.unlisted .unnumbered} +#### Assumes you know how to:{.unlisted .unnumbered} * Import images and image collections, filter, and visualize (Part F1). @@ -493,7 +469,7 @@ This chapter will cover the techniques for aggregating individual images from a * Inspect an Image and an ImageCollection, as well as their properties (Chap. F4.1). ::: -## Introduction {.unlisted .unnumbered} +### Introduction {.unlisted .unnumbered} CHIRPS is a high-resolution global gridded rainfall dataset that combines satellite-measured precipitation with ground station data in a consistent, long time-series dataset. The data are provided by the University of California, Santa Barbara, and are available from 1981 to the present. This dataset is extremely useful in drought monitoring and assessing global environmental change over land. The satellite data are calibrated with ground station observations to create the final product. @@ -501,7 +477,7 @@ CHIRPS is a high-resolution global gridded rainfall dataset that combines satell In this exercise, we will work with the CHIRPS dataset using the pentad. A pentad represents the grouping of five days. There are six pentads in a calendar month, with five pentads of exactly five days each and one pentad with the remaining three to six days of the month. Pentads reset at the beginning of each month, and the first day of every month is the start of a new pentad. Values at a given pixel in the CHIRPS dataset represent the total precipitation in millimeters over the pentad. -## Filtering an Image Collection +### Filtering an Image Collection We will start by accessing the CHIRPS Pentad collection and filtering it to create a time series for a single year. @@ -522,7 +498,7 @@ Each image’s pixel values store the total precipitation during the pentad. Wit :::{.callout-note} Code Checkpoint F42a. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Working with Dates +### Working with Dates To aggregate the time series, we need to learn how to create and manipulate dates programmatically. This section covers some functions from the ee.Date module that will be useful. @@ -555,7 +531,7 @@ We will use the millis function in the next section when we need to set the syst :::{.callout-note} Code Checkpoint F42b. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Aggregating Images +### Aggregating Images Now we can start aggregating the pentads into monthly sums. The process of aggregation has two fundamental steps. The first is to determine the beginning and ending dates of one time interval (in this case, one month), and the second is to sum up all of the values (in this case, the pentads) that fall within each interval. To begin, we can envision that the resulting series will contain 12 images. To prepare to create an image for each month, we create an ee.List of values from 1 to 12. We can use the ee.List.sequence function, as first presented in Chap. F1.0, to create the list of items of type ee.Number. Continuing with the script of the previous section, paste the following code: @@ -599,7 +575,7 @@ We have now successfully computed an aggregated collection from the source Image :::{.callout-note} Code Checkpoint F42c. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Plotting Time Series +### Plotting Time Series One useful application of gridded precipitation datasets is to analyze rainfall patterns. We can plot a time-series chart for a location using the newly computed time series. We can plot the pixel value at any given point or polygon. Here we create a point geometry for a given coordinate. Continuing with the script of the previous section, paste the following code: @@ -646,33 +622,11 @@ The customized chart (Fig. F4.2.3) shows the typical rainfall pattern in the cit :::{.callout-note} Code Checkpoint F42d. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Synthesis {.unnumbered} - -Assignment 1. The CHIRPS collection contains data for 40 years. Aggregate the same collection to yearly images and create a chart for annual precipitation from 1981 to 2021 at your chosen location. - -Instead of creating a list of months and writing a function to create monthly images, we will create a list of years and write a function to create yearly images. The code snippet below will help you get started. - -var chirps = ee.ImageCollection('UCSB-CHG/CHIRPS/PENTAD'); - -```js -// Create a list of years -var years = ee.List.sequence(1981, 2021); - -// Write a function that takes a year number -// and returns a yearly image -var createYearlyImage = function(beginningYear) { // Add your code -}; - -var yearlyImages = years.map(createYearlyImage); -var yearlyCollection = ee.ImageCollection.fromImages(yearlyImages); -print(yearlyCollection); - -``` -## Conclusion {.unnumbered} +### Conclusion {.unnumbered} In this chapter, you learned how to aggregate a collection to months and plot the resulting time series for a location. This chapter also introduced useful functions for working with the dates that will be used across many different applications. You also learned how to iterate over a list using the map function. The technique of mapping a function over a list or collection is essential for processing data. Mastering this technique will allow you to scale your analysis using the parallel computing capabilities of Earth Engine. -## References {.unnumbered} +### References {.unnumbered} Banerjee A, Chen R, Meadows ME, et al (2020) An analysis of long-term rainfall trends and variability in the Uttarakhand Himalaya using Google Earth Engine. Remote Sens 12:709. https://doi.org/10.3390/rs12040709 @@ -682,16 +636,16 @@ Okamoto K, Ushio T, Iguchi T, et al (2005) The global satellite mapping of preci -# Clouds and Image Compositing +## Clouds and Image Compositing :::{.callout-tip} -# Chapter Information +## Chapter Information -## Author {.unlisted .unnumbered} +#### Author {.unlisted .unnumbered} : @@ -699,19 +653,19 @@ Txomin Hermosilla, Saverio Francini, Andréa P. Nicolau, Michael A. Wulder, Joan -## Overview {.unlisted .unnumbered} +#### Overview {.unlisted .unnumbered} The purpose of this chapter is to provide necessary context and demonstrate different approaches for image composite generation when using data quality flags, using an initial example of removing cloud cover. We will examine different filtering options, demonstrate an approach for cloud masking, and provide additional opportunities for image composite development. Pixel selection for composite development can exclude unwanted pixels—such as those impacted by cloud, shadow, and smoke or haze—and can also preferentially select pixels based upon proximity to a target date or a preferred sensor type. -## Learning Outcomes {.unlisted .unnumbered} +#### Learning Outcomes {.unlisted .unnumbered} * Understanding and applying satellite-specific cloud mask functions. * Incorporating images from different sensors. * Using focal functions to fill in data gaps. -## Assumes you know how to:{.unlisted .unnumbered} +#### Assumes you know how to:{.unlisted .unnumbered} * Import images and image collections, filter, and visualize (Part F1). @@ -723,7 +677,7 @@ The purpose of this chapter is to provide necessary context and demonstrate diff * Summarize an ImageCollection with reducers (Chap. F4.0, Chap. F4.1). ::: -## Introduction {.unlisted .unnumbered} +### Introduction {.unlisted .unnumbered} In many respects, satellite remote sensing is an ideal source of data for monitoring large or remote regions. However, cloud cover is one of the most common limitations of optical sensors in providing continuous time series of data for surface mapping and monitoring. This is particularly relevant in tropical, polar, mountainous, and high-latitude areas, where clouds are often present. Many studies have addressed the extent to which cloudiness can restrict the monitoring of various regions (Zhu and Woodcock 2012, 2014; Eberhardt et al. 2016; Martins et al. 2018). @@ -747,7 +701,7 @@ The general workflow to generate a cloud-free composite involves: Additional steps may be necessary to improve the composite generated. These steps will be explained in the following sections. -## Cloud Filter and Cloud Mask +### Cloud Filter and Cloud Mask The first step is to define your AOI and center the map. The goal is to create a nationwide composite for the country of Colombia. We will use the Large Scale International Boundary (2017) simplified dataset from the US Department of State (USDOS), which contains polygons for all countries of the world. @@ -865,7 +819,7 @@ The resulting image has substantially fewer data gaps (you can zoom in to better :::{.callout-note} Code Checkpoint F43a. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Incorporating Data from Other Satellites +### Incorporating Data from Other Satellites Another option to reduce the presence of data gaps in cloudy situations is to bring in imagery from other sensors acquired during the time period of interest. The Landsat collection spans multiple missions, which have continuously acquired uninterrupted data since 1972 at different acquisition dates. Next, we will try incorporating Landsat 7 Level 2, Collection 2, Tier 1 images from 2019 to fill the gaps in the 2019 Landsat 8 composite. @@ -959,7 +913,7 @@ Comparing the composite generated considering both Landsat 7 and 8 to the Landsa :::{.callout-note} Code Checkpoint F43b. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Best-Available-Pixel Compositing Earth Engine Application +### Best-Available-Pixel Compositing Earth Engine Application This section presents an Earth Engine application that enables the generation of annual best-available-pixel (BAP) image composites by globally combining multiple Landsat sensors and images: GEE-BAP. Annual BAP image composites are generated by choosing optimal observations for each pixel from all available Landsat 5 TM, Landsat 7 ETM+, and Landsat 8 OLI imagery within a given year and within a given day range from a specified acquisition day of year, in addition to other constraints defined by the user. The data accessible via Earth Engine are from the USGS free and open archive of Landsat data. The Landsat images used are atmospherically corrected to surface reflectance values. Following White et al. (2014), a series of scoring functions ranks each pixel observation for (1) acquisition day of year, (2) cloud cover in the scene, (3) distance to clouds and cloud shadows, (4) presence of haze, and (5) acquisition sensor. Further information on the BAP image compositing approach can be found in Griffiths et al. (2013), and detailed information on tuning parameters can be found in White et al. (2014). @@ -1031,19 +985,11 @@ library.AddSLider(startYear, endYear); :::{.callout-note} Code Checkpoint F43d. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Synthesis {.unnumbered} - -Assignment 1. Create composites for other cloudy regions or less cloudy regions. For example, change the country variable to 'Cambodia' or 'Mozambique'. Are more gaps present in the resulting composite? Can you change the compositing rules to improve this (using Acquisition day of year and Day range)? Different regions of the Earth have different cloud seasonal patterns, so the most appropriate date windows to acquire cloud-free composites will change depending on location. Also be aware that the larger the country, the longer it will take to generate the composite. - -Assignment 2. Similarly, try creating composites for the wet and dry seasons of a region separately. Compare the two composites. Are some features brighter or darker? Is there evidence of drying of vegetation, such as related to leaf loss or reduction in herbaceous ground vegetation? - -Assignment 3. Test different cloud threshold values and see if you can find an optimal threshold that balances data gaps against area coverage for your particular target date. - -## Conclusion {.unnumbered} +### Conclusion {.unnumbered} We cannot monitor what we cannot see. Image compositing algorithms provide robust and transparent tools to address issues with clouds, cloud shadows, haze, and smoke in remotely sensed images derived from optical satellite data, and expand data availability for remote sensing applications. The tools and approaches described here should provide you with some useful strategies to aid in mitigating the presence of cloud cover in your data. Note that the quality of image outcomes is a function of the quality of cloud masking routines applied to the source data to generate the various flags that are used in the scoring functions described herein. Different compositing parameters can be used to represent a given location as a function of conditions that are present at a given point in time and the information needs of the end user. Tuning or optimization of compositing parameters is possible (and recommended) to ensure best capture of the physical conditions of interest. -## References {.unnumbered} +### References {.unnumbered} Braaten JD, Cohen WB, Yang Z (2015) Automated cloud and cloud shadow identification in Landsat MSS imagery for temperate ecosystems. Remote Sens Environ 169:128–138. https://doi.org/10.1016/j.rse.2015.08.006 @@ -1091,14 +1037,14 @@ Zhu Z, Woodcock CE (2012) Object-based cloud and cloud shadow detection in Lands -# Change Detection +## Change Detection :::{.callout-tip} -# Chapter Information +## Chapter Information -## Author {.unlisted .unnumbered} +#### Author {.unlisted .unnumbered} @@ -1106,12 +1052,12 @@ Karis Tenneson, John Dilger, Crystal Wespestad, Brian Zutta, Andréa P Nicolau, -## Overview {.unlisted .unnumbered} +#### Overview {.unlisted .unnumbered} This chapter introduces change detection mapping. It will teach you how to make a two-date land cover change map using image differencing and threshold-based classification. You will use what you have learned so far in this book to produce a map highlighting changes in the land cover between two time steps. You will first explore differences between the two images extracted from these time steps by creating a difference layer. You will then learn how to directly classify change based on the information in both of your images. -## Learning Outcomes {.unlisted .unnumbered} +#### Learning Outcomes {.unlisted .unnumbered} * Creating and exploring how to read a false-color cloud-free Landsat composite image @@ -1119,14 +1065,14 @@ This chapter introduces change detection mapping. It will teach you how to make * Creating a two-image difference to help locate areas of change. * Producing a change map and classifying changes using thresholding. -## Assumes you know how to:{.unlisted .unnumbered} +#### Assumes you know how to:{.unlisted .unnumbered} * Import images and image collections, filter, and visualize (Part F1). * Perform basic image analysis: select bands, compute indices, create masks (Part F2). ::: -## Introduction {.unlisted .unnumbered} +### Introduction {.unlisted .unnumbered} Change detection is the process of assessing how landscape conditions are changing by looking at differences in images acquired at different times. This can be used to quantify changes in forest cover—such as those following a volcanic eruption, logging activity, or wildfire—or when crops are harvested (Fig. F4.4.1). For example, using time-series change detection methods, Hansen et al. (2013) quantified annual changes in forest loss and regrowth. Change detection mapping is important for observing, monitoring, and quantifying changes in landscapes over time. Key questions that can be answered using these techniques include identifying whether a change has occurred, measuring the area or the spatial extent of the region undergoing change, characterizing the nature of the change, and measuring the pattern (configuration or composition) of the change (MacLeod and Congalton 1998). @@ -1160,7 +1106,7 @@ For the practicum, you will select pre-event and post-event image scenes and inv ![Fig. F4.4.2 Change detection workflow for this practicum](F4/image68.png) -## Preparing Imagery +### Preparing Imagery Before beginning a change detection workflow, image preprocessing is essential. The goal is to ensure that each pixel records the same type of measurement at the same location over time. These steps include multitemporal image registration and radiometric and atmospheric corrections, which are especially important. A lot of this work has been automated and already applied to the images that are available in Earth Engine. Image selection is also important. Selection considerations include finding images with low cloud cover and representing the same phenology (e.g., leaf-on or leaf-off). @@ -1191,7 +1137,7 @@ var preImage = landsat8 .sort('CLOUD_COVER', true) .first(); -## Creating False-Color Composites +### Creating False-Color Composites Before running any sort of change detection analysis, it is useful to first visualize your input images to get a sense of the landscape, visually inspect where changes might occur, and identify any problems in the inputs before moving further. As described in Chap. F1.1, false-color composites draw bands from multispectral sensors in the red, green, and blue channels in ways that are designed to illustrate contrast in imagery. Below, you will produce a false-color composite using SWIR-2 in the red channel, NIR in the green channel, and Red in the blue channel (Fig. F4.4.3). @@ -1205,7 +1151,7 @@ Map.addLayer(postImage, visParam, 'post'); ![Fig. F4.4.3 False-color composite using SWIR2, NIR, and red. Vegetation shows up vividly in the green channel due to vegetation being highly reflective in the NIR band. Shades of green can be indicative of vegetation density; water typically shows up as black to dark blue; and burned or barren areas show up as brown.](F4/image31.png) -## Calculating NBR +### Calculating NBR The next step is data transformation, such as calculating NBR. The advantage of using these techniques is that the data, along with the noise inherent in the data, have been reduced in order to simplify a comparison between two images. Image differencing is done by subtracting the spectral value of the first-date image from that of the second-date image, pixel by pixel (Fig. F4.4.2). Two-date image differencing can be used with a single band or with spectral indices, depending on the application. Identifying the correct band or index to identify change and finding the correct thresholds to classify it are critical to producing meaningful results. Working with indices known to highlight the land cover conditions before and after a change event of interest is a good starting point. For example, the Normalized Difference Water Index would be good for mapping water level changes during flooding events; the NBR is good at detecting soil brightness; and the NDVI can be used for tracking changes in vegetation (although this index does saturate quickly). In some cases, using derived band combinations that have been customized to represent the phenomenon of interest is suggested, such as using the Normalized Difference Fraction Index to monitor forest degradation (see Chap. A3.4). @@ -1224,7 +1170,7 @@ var nbrPost = postImage.normalizedDifference(['nir', 'swir2']) :::{.callout-note} Code Checkpoint F44a. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Single Date Transformation +### Single Date Transformation Next, we will examine the changes that have occurred, as seen when comparing two specific dates in time. @@ -1253,7 +1199,7 @@ a) b) ![Fig. F4.4.4 (a) Two-date NBR difference; (b) pre-event image (June 2013) false-color composite; (c) post-event image (June 2020) false-color composite. In the change map (a), areas on the lower range of values (blue) depict areas where vegetation has been negatively affected, and areas on the higher range of values (pink) depict areas where there has been vegetation gain; the green/orange areas have experienced little change. In the pre-event and post-event images (b and c), the green areas indicate vegetation, while the brown regions are barren ground.](F4/image4.png) -## Classifying Change +### Classifying Change Once the images have been transformed and differenced to highlight areas undergoing change, the next step is image classification into a thematic map consisting of stable and change classes. This can be done rather simply by thresholding the change layer, or by using classification techniques such as machine learning algorithms. One challenge of working with simple thresholding of the difference layer is knowing how to select a suitable threshold to partition changed areas from stable classes. On the other hand, classification techniques using machine learning algorithms partition the landscape using examples of reference data that you provide to train the classifier. This may or may not yield better results, but does require additional work to collect reference data and train the classifier. In the end, resources, timing, and the patterns of the phenomenon you are trying to map will determine which approach is suitable—or perhaps the activity you are trying to track requires something more advanced, such as a time-series approach that uses more than two dates of imagery. @@ -1295,23 +1241,11 @@ Chapters F4.5 through F4.9 present more-advanced change detection algorithms tha :::{.callout-note} Code Checkpoint F44b. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Synthesis {.unnumbered} - -Evaluating any maps you create, including change detection maps, is essential to determining whether the method you have selected is appropriate for informing land management and decision-making (Stehman and Czaplewski 1998), or whether you need to iterate on the mapping process to improve the final results. Maps generally, and change maps specifically, will always have errors. This is due to a suite of factors, such as the geometric registration between images, the calibration between images, the data resolution (e.g., temporal, spectral, radiometric) compared to the characteristics of the activity of interest, the complexity of the landscape of the study region (topography, atmospheric conditions, etc.), and the classification techniques employed (Lu et al. 2004). This means that similar studies can present different, sometimes controversial, conclusions about landscape dynamics (e.g., Cohen et al. 2017). In order to be useful for decision-making, a change detection mapping effort should provide the user with an understanding of the strengths and weaknesses of the product, such as by presenting omission and commission error rates. The quantification of classification quality is presented in Chap. F2.2. - -Assignment 1. Try using a different index, such as NDVI or a Tasseled Cap Transformation, to run the change detection steps, and compare the results with those obtained from using NBR. - -Assignment 2. Experiment with adjusting the thresholdLoss and thresholdGain values. - -Assignment 3. Use what you have learned in the classification chapter (Chap. F2.1) to run a supervised classification on the difference layer (or layers, if you have created additional ones). Hint: To complete a supervised classification, you would need reference examples of both the stable and change classes of interest to train the classifier. - -Assignment 4. Think about how things like clouds and cloud shadows could affect the results of change detection. What do you think the two-date differencing method would pick up for images in the same year in different seasons? - -## Conclusion {.unnumbered} +### Conclusion {.unnumbered} In this chapter, you learned how to make a change detection map using two-image differencing. The importance of visualizing changes in this way instead of using a post-classification comparison, where two classified maps are compared instead of two satellite images, is that it avoids multiplicative errors from the classifications and is better at observing more subtle changes in the landscape. You also learned that how you visualize your images and change maps—such as what band combinations and color ramps you select, and what threshold values you use for a classification map—has an impact on how easily and what types of changes can be seen. -## References {.unnumbered} +### References {.unnumbered} Cohen WB, Healey SP, Yang Z, et al (2017) How similar are forest disturbance maps derived from different Landsat time series algorithms? Forests 8:98. https://doi.org/10.3390/f8040098 @@ -1337,16 +1271,16 @@ Woodcock CE, Loveland TR, Herold M, Bauer ME (2020) Transitioning from change de -# Interpreting Annual Time Series with LandTrendr +## Interpreting Annual Time Series with LandTrendr :::{.callout-tip} -# Chapter Information +## Chapter Information -## Author {.unlisted .unnumbered} +#### Author {.unlisted .unnumbered} @@ -1354,12 +1288,12 @@ Robert Kennedy, Justin Braaten, Peter Clary -## Overview {.unlisted .unnumbered} +#### Overview {.unlisted .unnumbered} Time-series analysis of change can be achieved by fitting the entire spectral trajectory using simple statistical models. These allow us to both simplify the time series and to extract useful information about the changes occurring. In this chapter, you will get an introduction to the use of LandTrendr, one of these time-series approaches used to characterize time series of spectral values. -## Learning Outcomes {.unlisted .unnumbered} +#### Learning Outcomes {.unlisted .unnumbered} * Evaluating yearly time-series spectral values to distinguish between true change and artifacts. @@ -1367,14 +1301,14 @@ Time-series analysis of change can be achieved by fitting the entire spectral tr * Interpreting change segments and translating them to maps. * Applying parameters in a graphical user interface to create disturbance maps in forests. -## Assumes you know how to:{.unlisted .unnumbered} +#### Assumes you know how to:{.unlisted .unnumbered} * Calculate and interpret vegetation indices (Chap. F2.0) * Interpret bands and indices in terms of land surface characteristics (Chap. F2.0). ::: -## Introduction {.unlisted .unnumbered} +### Introduction {.unlisted .unnumbered} Land surface change happens all the time, and satellite sensors witness it. If a spectral index is chosen to match the type of change being sought, surface change can be inferred from changes in spectral index values. Over time, the progression of spectral values witnessed in each pixel tells a story of the processes of change, such as growth and disturbance. Time-series algorithms are designed to leverage many observations of spectral values over time to isolate and describe changes of interest, while ignoring uninteresting change or noise. @@ -1387,7 +1321,7 @@ For this lab, we will use a graphical user interface (GUI) to teach the concepts :::{.callout-note} Code Checkpoint F45a. The book’s repository contains information about accessing the LandTrendr interface. ::: -## Pixel Time Series +### Pixel Time Series When working with LandTrendr for the first time in your area, there are two questions you must address. @@ -1457,7 +1391,7 @@ Note: because of temporal autocorrelation, these cannot be interpreted as true f From the models that pass the p-value threshold, one is chosen as the final fit. It may be the one with the lowest p-value. However, an adjustment is made to allow more complicated models (those with more segments) to be picked even if their p-value is within a defined proportion of the best-scoring model. That proportion is set by the best model proportion parameter. As an example, a best model proportion value of 0.75 would allow a more complicated model to be chosen if its score were greater than 75% that of the best model. -## Translating Pixels to Maps +### Translating Pixels to Maps Although the full time series is the best description of each pixel’s “life history,” we typically are interested in the behavior of all of the pixels in our study area. It would be both inefficient to manually visualize all of them and ineffective to try to summarize areas and locations. Thus, we seek to make maps. @@ -1503,45 +1437,13 @@ There are multiple layers of disturbance added to the map. Use the map layers ch ![Fig. 4.5.8 Magnitude of change for the same area](F4/image14.png) -## Synthesis {.unnumbered} - -In this chapter, you have learned how to work with annual time series to interpret regions of interest. Looking at annual snapshots of the landscape provides three key benefits: (1) the ability to view your area of interest without the clouds and noise that typically obscure single-day views; (2) gauge the amount by which the noise-dampened signal still varies from year to year in response to large-scale forcing mechanisms; and (3) the ability to view the response of landscapes as they recover, sometimes slowly, from disturbance. - -To learn more about LandTrendr, see the assignments below. - -Assignment 1. Find your own change processes of interest. First, navigate the map (zooming and dragging) to an area of the world where you are interested in a change process, and the spectral index that would capture it. Make sure the UI control panel is open to the Pixel Time-Series Options section. Next, click on the map in areas where you know change has occurred, and observe the spectral trajectories in the charts. Then, describe whether the change of interest is detectable in the spectral reflectance record, and what are its characteristics in different parts of the study area. . - -Assignment 2: Find a pixel in your area of interest that shows a distinctive disturbance process, as you define it for your topic of interest. Adjust date ranges, parameters, etc. using the steps outlined in Section 1 above, and then answer these questions: - -* Question 1. Which index and date range did you use? -* Question 2. Did you need to change fitting parameters to make the algorithm find the disturbance? If so, which ones, and why? -* Question 3. How do you know this is a disturbance? - -Assignment 3. Switch the control panel in the GUI to Change Filter Options, and use the guidance in Section 2 to set parameters and make disturbance maps. - -* Question 4. Do the disturbance year and magnitude as mapped in the image match with what you would expect from the trajectory itself? -* Question 5. Can you change some of the filters to create a map where your disturbance process is not mapped? If so, what did you change? -* Question 6. Can you change filters to create a map that includes a different disturbance process, perhaps subtler, longer duration, etc.? Find a pixel and use the “Pixel time series” plotter to look at the time series of those processes. - -Assignment 4: Return to the Pixel Time-Series Options section of the control panel, and navigate to a pixel in your area of interest that you believe would show a distinctive recovery or growth process, as you define it for your topic of interest. You may want to modify the index, parameters, etc. as covered in Section 1 to adequately capture the growth process with the fitted trajectories. - -* Question 7. Did you use the same spectral index? If not, why? -* Question 8. Were the fitting parameters the same as those for disturbance? If not, what did you change, and why? -* Question 9. What evidence do you have that this is a vegetative growth signal? - -Assignment 5. For vegetation gain mapping, switch the control panel back to Change Filter Options and use the guidance in Section 2 to set parameters, etc. to make maps of growth. - -* Question 10. For the pixel or pixels you found for Assignment 3, does the year and magnitude as mapped in the “gain” image match with what you would expect from the trajectory itself? -* Question 11. Compare what the map looks like when you run it with and without the MMU filter. What differences do you see? -* Question 12. Try changing the recovery duration filter to a very high number (perhaps the full length of your archive) and to a very low number (say, one or two years). What differences do you see? - -## Conclusion {.unnumbered} +### Conclusion {.unnumbered} This exercise provides a baseline sense of how the LandTrendr algorithm works. The key points are learning how to interpret change in spectral values in terms of the processes occurring on the ground, and then translating those into maps. You can export the images you’ve made here using Download Options. Links to materials are available in the chapter checkpoints and LandTrendr documentation about both the GUI and the script-based versions of the algorithm. In particular, there are scripts that handle different components of the fitting and mapping process, and that allow you to keep track of the fitting and image selection criteria. -## References {.unnumbered} +### References {.unnumbered} Kennedy RE, Yang Z, Cohen WB (2010) Detecting trends in forest disturbance and recovery using yearly Landsat time series: 1. LandTrendr - Temporal segmentation algorithms. Remote Sens Environ 114:2897–2910. https://doi.org/10.1016/j.rse.2010.07.008 @@ -1549,16 +1451,16 @@ Kennedy RE, Yang Z, Gorelick N, et al (2018) Implementation of the LandTrendr al -# Fitting Functions to Time Series +## Fitting Functions to Time Series :::{.callout-tip} -# Chapter Information +## Chapter Information -## Author {.unlisted .unnumbered} +#### Author {.unlisted .unnumbered} @@ -1566,19 +1468,19 @@ Andréa Puzzi Nicolau, Karen Dyson, Biplov Bhandari, David Saah, Nicholas Clinto -## Overview {.unlisted .unnumbered} +#### Overview {.unlisted .unnumbered} The purpose of this chapter is to establish a foundation for time-series analysis of remotely sensed data, which is typically arranged as an ordered stack of images. You will be introduced to the concepts of graphing time series, using linear modeling to detrend time series, and fitting harmonic models to time-series data. At the completion of this chapter, you will be able to perform analysis of multi-temporal data for determining trend and seasonality on a per-pixel basis. -## Learning Outcomes {.unlisted .unnumbered} +#### Learning Outcomes {.unlisted .unnumbered} * Graphing satellite imagery values across a time series. * Quantifying and potentially removing linear trends in time series. * Fitting linear and harmonic models to individual pixels in time-series data. -## Assumes you know how to:{.unlisted .unnumbered} +#### Assumes you know how to:{.unlisted .unnumbered} * Import images and image collections, filter, and visualize (Part F1). @@ -1589,7 +1491,7 @@ The purpose of this chapter is to establish a foundation for time-series analysi * Mask cloud, cloud shadow, snow/ice, and other undesired pixels (Chap. F4.3). ::: -## Introduction {.unlisted .unnumbered} +### Introduction {.unlisted .unnumbered} Many natural and man-made phenomena exhibit important annual, interannual, or longer-term trends that recur—that is, they occur at roughly regular intervals. Examples include seasonality in leaf patterns in deciduous forests and seasonal crop growth patterns. Over time, indices such as the Normalized Difference Vegetation Index (NDVI) will show regular increases (e.g., leaf-on, crop growth) and decreases (e.g., leaf-off, crop senescence), and typically have a long-term, if noisy, trend such as a gradual increase in NDVI value as an area recovers from a disturbance. @@ -1599,7 +1501,7 @@ Earth Engine supports the ability to do complex linear and non-linear regression If you have not already done so, be sure to add the book’s code repository to the Code Editor by entering [https://code.earthengine.google.com/?accept_repo=projects/gee-edu/book](https://www.google.com/url?q=https://code.earthengine.google.com/?accept_repo%3Dprojects/gee-edu/book&sa=D&source=editors&ust=1671458868327546&usg=AOvVaw0m0oT1feMKQKRq3rtuOxfY) into your browser. The book’s scripts will then be available in the script manager panel. -## Multi-Temporal Data in Earth Engine +### Multi-Temporal Data in Earth Engine As explained in Chaps. F4.0 and F4.1, a time series in Earth Engine is typically represented as an ImageCollection. Because of image overlaps, cloud treatments, and filtering choices, an ImageCollection can have any of the following complex characteristics: @@ -1618,7 +1520,7 @@ First, we will give some very basic notation (Fig. F4.6.1). A scalar pixel at ti ![Fig. F4.6.1 Time series representation of pixel p ](F4/image38.png) -## Data Preparation and Preprocessing +### Data Preparation and Preprocessing The first step in analysis of time-series data is to import data of interest and plot it at an interesting location. We will work with the USGS Landsat 8 Level 2, Collection 2, Tier 1 ImageCollection and a cloud-masking function (Chap. F4.3), scale the image values, and add variables of interest to the collection as bands. Copy and paste the code below to filter the Landsat 8 collection to a point of interest over California (variable roi) and specific dates, and to apply the defined function. The variables of interest added by the function are: (1) NDVI (Chap. F2.0), (2) a time variable that is the difference between the image’s current year and the year 1970 (a start point), and (3) a constant variable with value 1. @@ -1693,7 +1595,7 @@ Now that we have plotted and visualized the data, lots of interesting analyses c :::{.callout-note} Code Checkpoint F46a. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Estimating Linear Trend Over Time +### Estimating Linear Trend Over Time Time series datasets may contain not only trends but also seasonality, both of which may need to be removed prior to modeling. Trends and seasonality can result in a varying mean and a varying variance over time, both of which define a time series as non-stationary. Stationary datasets, on the other hand, have a stable mean and variance, and are therefore much easier to model. @@ -1763,7 +1665,7 @@ var detrendedChart = ui.Chart.image.series(detrended, roi, null, 30) :::{.callout-note} Code Checkpoint F46b. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Estimating Seasonality with a Harmonic Model +### Estimating Seasonality with a Harmonic Model A linear trend is one of several possible types of trends in time series. Time series can also present harmonic trends, in which a value goes up and down in a predictable wave pattern. These are of particular interest and usefulness in the natural world, where harmonic changes in greenness of deciduous vegetation can occur across the spring, summer, and autumn. Now we will return to the initial time series (landsat8sr) of Fig. F4.6.2 and fit a harmonic pattern through the data. Consider the following harmonic model, where A is amplitude, ω is frequency, φ is phase, and et is a random error. @@ -1858,9 +1760,9 @@ The code uses the HSV to RGB transformation hsvToRgb for visualization purposes :::{.callout-note} Code Checkpoint F46c. The book’s repository contains a script that shows what your code should look like at this point. ::: -## An Application of Curve Fitting +### An Application of Curve Fitting -### The rich data about the curve fits can be viewed in a multitude of different ways. Add the code below to your script to produce the view in Fig. 4.6.9. The image will be a close-up of the area around Modesto, California. +#### The rich data about the curve fits can be viewed in a multitude of different ways. Add the code below to your script to produce the view in Fig. 4.6.9. The image will be a close-up of the area around Modesto, California. ```js ///////////////////// Section 5 ///////////////////////////// @@ -1894,7 +1796,7 @@ The upper image in Fig. F4.6.9 is a closer view of Fig. F4.6.8, showing an image :::{.callout-note} Code Checkpoint F46d. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Higher-Order Harmonic Models +### Higher-Order Harmonic Models Harmonic models are not limited to fitting a single wave through a set of points. In some situations, there may be more than one cycle within a given year—for example, when an agricultural field is double-cropped. Modeling multiple waves within a given year can be done by adding more harmonic terms to Eq. F4.6.2. The code at the following checkpoint allows the fitting of any number of cycles through a given point. @@ -1910,29 +1812,11 @@ Beginning with the repository script, changing the value of the harmonics variab ![Fig. F4.6.10 Fit with harmonic curves of increasing complexity, fitted for data at a given point ](F4/image64.png) -## Synthesis {.unnumbered} - -Assignment 1. Fit two NDVI harmonic models for a point close to Manaus, Brazil: one prior to a disturbance event and one after the disturbance event (Fig. F4.6.11). You can start with the code checkpoint below, which gives you the point coordinates and defines the initial functions needed. The disturbance event happened in mid-December 2014, so set filter dates for the first ImageCollection to '2013-01-01','2014-12-12', and set the filter dates for the second ImageCollection to '2014-12-13','2019-01-01'. Merge both fitted collections and plot both NDVI and fitted values. The result should look like Fig. F4.6.12. - -:::{.callout-note} -Code Checkpoint F46s1. The book’s repository contains a script that shows what your code should look like at this point. -::: -![Fig. F4.6.11 Landsat 8 images showing the land cover change at a point in Manaus, Brazil; (left) July, 6, 2014, (right) August 8, 2015](F4/image26.png) - - -![Fig. F4.6.12 Fitted harmonic models before and after disturbance events to a given point in the Brazilian Amazon](F4/image36.png) - - -What do you notice? Think about how the harmonic model would look if you tried to fit the entire period. In this example, you were given the date of the breakpoint between the two conditions of the land surface within the time series. State-of-the-art land cover change algorithms work by assessing the difference between the modeled and observed pixel values. These algorithms look for breakpoints in the model, typically flagging changes after a predefined number of consecutive observations. - -:::{.callout-note} -Code Checkpoint F46s2. The book’s repository contains a script that shows what your code should look like at this point. -::: -## Conclusion {.unnumbered} +### Conclusion {.unnumbered} In this chapter, we learned how to graph and fit both linear and harmonic functions to time series of remotely sensed data. These skills underpin important tools such as Continuous Change Detection and Classification (CCDC, Chap. F4.7) and Continuous Degradation Detection (CODED, Chap. A3.4). These approaches are used by many organizations to detect forest degradation and deforestation (e.g., Tang et al. 2019, Bullock et al. 2020). These approaches can also be used to identify crops (Chap. A1.1) with high degrees of accuracy (Ghazaryan et al. 2018). -## References {.unnumbered} +### References {.unnumbered} Bradley BA, Jacob RW, Hermance JF, Mustard JF (2007) A curve fitting procedure to derive inter-annual phenologies from time series of noisy satellite NDVI data. Remote Sens Environ 106:137–145. https://doi.org/10.1016/j.rse.2006.08.002 @@ -1944,16 +1828,16 @@ Shumway RH, Stoffer DS (2019) Time Series: A Data Analysis Approach Using R. Cha Tang X, Bullock EL, Olofsson P, et al (2019) Near real-time monitoring of tropical forest disturbance: New algorithms and assessment framework. Remote Sens Environ 224:202–218. https://doi.org/10.1016/j.rse.2019.02.003 -# Interpreting Time Series with CCDC +## Interpreting Time Series with CCDC :::{.callout-tip} -# Chapter Information +## Chapter Information -## Author {.unlisted .unnumbered} +#### Author {.unlisted .unnumbered} @@ -1961,12 +1845,12 @@ Paulo Arévalo, Pontus Olofsson -## Overview {.unlisted .unnumbered} +#### Overview {.unlisted .unnumbered} Continuous Change Detection and Classification (CCDC) is a land change monitoring algorithm designed to operate on time series of satellite data, particularly Landsat data. This chapter focuses on the portion that is the change detection component (CCD); you will learn how to run the algorithm, interpret its outputs, and visualize coefficients and change information. -## Learning Outcomes {.unlisted .unnumbered} +#### Learning Outcomes {.unlisted .unnumbered} * Exploring pixel-level time series of Landsat observations, as well as the temporal segments that CCDC fits to the observations. @@ -1976,7 +1860,7 @@ Continuous Change Detection and Classification (CCDC) is a land change monitorin * Using array image functions. * Attaching user-defined metadata to an image when exporting. -## Assumes you know how to:{.unlisted .unnumbered} +#### Assumes you know how to:{.unlisted .unnumbered} * Import images and image collections, filter, and visualize (Part F1). @@ -1986,12 +1870,12 @@ Continuous Change Detection and Classification (CCDC) is a land change monitorin * Work with array images (Chap. F3.1, Chap. F4.6). * Interpret fitted harmonic models (Chap. F4.6). -## Introduction to Theory +### Introduction to Theory “A time series is a sequence of observations taken sequentially in time. … An intrinsic feature of a time series is that, typically, adjacent observations are dependent. Time-series analysis is concerned with techniques for the analysis of this dependency.” This is the formal definition of time-series analysis by Box et al. (1994). In a remote sensing context, the observations of interest are measurements of radiation reflected from the surface of the Earth from the Sun or an instrument emitting energy toward Earth. Consecutive measurements made over a given area result in a time series of surface reflectance. By analyzing such time series, we can achieve a comprehensive characterization of ecosystem and land surface processes (Kennedy et al. 2014). The result is a shift away from traditional, retrospective change-detection approaches based on data acquired over the same area at two or a few points in time to continuous monitoring of the landscape (Woodcock et al. 2020). Previous obstacles related to data storage, preprocessing, and computing power have been largely overcome with the emergence of powerful cloud-computing platforms that provide direct access to the data (Gorelick et al. 2017). In this chapter, we will illustrate how to study landscape dynamics in the Amazon river basin by analyzing dense time series of Landsat data using the CCDC algorithm. Unlike LandTrendr (Chap. F4.5), which uses anniversary images to fit straight line segments that describe the spectral trajectory over time, CCDC uses all available clear observations. This has multiple advantages, including the ability to detect changes within a year and capture seasonal patterns, although at the expense of much higher computational demands and more complexity to manipulate the outputs, compared to LandTrendr. -## Understanding Temporal Segmentation with CCDC +### Understanding Temporal Segmentation with CCDC Spectral change is detected at the pixel level by testing for structural breaks in a time series of reflectance. In Earth Engine, this process is referred to as “temporal segmentation,” as pixel-level time series are segmented according to periods of unique reflectance. It does so by fitting harmonic regression models to all spectral bands in the time series. The model-fitting starts at the beginning of the time series and moves forward in time in an “online” approach to change detection. The coefficients are used to predict future observations, and if the residuals of future observations exceed a statistical threshold for numerous consecutive observations, then the algorithm flags that a change has occurred. After the change, a new regression model is fit and the process continues until the end of the time series. The details of the original algorithm are described in Zhu and Woodcock (2014). We have created an interface-based tool (Arévalo et al. 2020) that facilitates the exploration of time series of Landsat observations and the CCDC results. @@ -2007,7 +1891,7 @@ Pay special attention to the characteristics of each segment. For example, look Question 1. While still using the SWIR1 band, click on a pixel that is forested. What do the time series and time segments look like? -## Running CCDC +### Running CCDC The tool shown above is useful for understanding the temporal dynamics for a specific point. However, we can do a similar analysis for larger areas by first running the CCDC algorithm over a group of pixels. The CCDC function in Earth Engine can take any ImageCollection, ideally one with little or no noise, such as a Landsat ImageCollection where clouds and cloud shadows have been masked. CCDC contains an internal cloud masking algorithm and is rather robust against missed clouds, but the cleaner the data the better. To simplify the process, we have developed a function library that contains functions for generating input data and processing CCDC results. Paste this line of code in a new script: @@ -2082,7 +1966,7 @@ As a general rule, try to use pre-existing CCDC results if possible, and if you :::{.callout-note} Code Checkpoint F47b. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Extracting Break Information +### Extracting Break Information We will now start exploring the pre-exported CCDC results mentioned in the previous section. We will make use of the third-party module palettes, described in detail in Chap. F6.0, that simplifies the use of palettes for visualization. Paste the following code in a new script: @@ -2206,7 +2090,7 @@ Question 2. Compare the “first change” and “last change” layers with the Question 3. Looking at the “max magnitude of change” layer, find places showing the largest and the smallest values. What type of changes do you think are happening in each of those places? -## Extracting Coefficients Manually +### Extracting Coefficients Manually In addition to the change information generated by the CCDC algorithm, we can use the coefficients of the time segments for multiple purposes, like land cover classification. Each time segment can be described as a harmonic function with an intercept, slope, and three pairs of sine and cosine terms that allow the time segments to represent seasonality occurring at different temporal scales. These coefficients, as well as the root-mean-square error (RMSE) obtained by comparing each predicted and actual Landsat value, are produced when the CCDC algorithm is run. The following example will show you how to retrieve the intercept coefficient for a segment intersecting a specific date. In a new script, paste the code below: @@ -2380,19 +2264,11 @@ Finally, the RGB image we created is shown in Fig. 4.7.6. The intercepts are cal :::{.callout-note} Code Checkpoint F47e. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Synthesis {.unnumbered} - -Assignment 1. Use the time series from the first section of this chapter to explore the time series and time segments produced by CCDC in many locations around the world. Compare places with different land cover types, and places with more stable dynamics (e.g., lakes, primary forests) vs. highly dynamic places (e.g., agricultural lands, construction sites). Pay attention to the variability in data density across continents and latitudes, and the effect that data density has on the appearance of the time segments. Use different spectral bands and indices and notice how they capture the temporal dynamics you are observing. - -Assignment 2. Pick three periods within the temporal study period of the CCDC results we used earlier: one near to the start, another in the middle, and the third close to the end. For each period, visualize the maximum change magnitude. Compare the spatial patterns between periods, and reflect on the types of disturbances that might be happening at each stage. - -Assignment 3. Select the intercept coefficients of the middle date of each of the periods you chose in the previous assignment. For each of those dates, load an RGB image with the band combination of your choosing (or simply use the Red, Green and Blue intercepts to obtain true-color images). Using the Inspector tab, compare the values across images in places with subtle and large differences between them, as well as in areas that do not change. What do the values tell you in terms of the benefits of using CCDC to study changes in a landscape? - -## Conclusion {.unnumbered} +### Conclusion {.unnumbered} This chapter provided a guide for the interpretation of the results from the CCDC algorithm for studying deforestation in the Amazon. Consider the advantages of such an analysis compared to traditional approaches to change detection, which are typically based on the comparison of two or a few images collected over the same area. For example, with time-series analysis, we can study trends and subtle processes such as vegetation recovery or degradation, determine the timing of land-surface events, and move away from retrospective analyses to monitoring in near-real time. Through the use of all available clear observations, CCDC can detect intra-annual breaks and capture seasonal patterns, although at the expense of increased computational requirements and complexity, unlike faster and easier to interpret methods based on annual composites, such as LandTrendr (Chap. F4.5). We expect to see more applications that make use of multiple change detection approaches (also known as “Ensemble” approaches), and multisensor analyses in which data from different satellites are fused (radar and optical, for example) for higher data density. -## References {.unnumbered} +### References {.unnumbered} Arévalo P, Bullock EL, Woodcock CE, Olofsson P (2020) A suite of tools for continuous land change monitoring in Google Earth Engine. Front Clim 2. https://doi.org/10.3389/fclim.2020.576740 @@ -2408,16 +2284,16 @@ Zhu Z, Woodcock CE (2014) Continuous change detection and classification of land -# Data Fusion: Merging Classification Streams +## Data Fusion: Merging Classification Streams :::{.callout-tip} -# Chapter Information +## Chapter Information -## Author {.unlisted .unnumbered} +#### Author {.unlisted .unnumbered} @@ -2425,19 +2301,19 @@ Jeffrey A. Cardille, Rylan Boothman, Mary Villamor, Elijah Perez, Eidan Willis, -## Overview {.unlisted .unnumbered} +#### Overview {.unlisted .unnumbered} As the ability to rapidly produce classifications of satellite images grows, it will be increasingly important to have algorithms that can sift through them to separate the signal from inevitable classification noise. The purpose of this chapter is to explore how to update classification time series by blending information from multiple classifications made from a wide variety of data sources. In this lab, we will explore how to update the classification time series of the Roosevelt River found in Fortin et al. (2020). That time series began with the 1972 launch of Landsat 1, blending evidence from 10 sensors and more than 140 images to show the evolution of the area until 2016. How has it changed since 2016? What new tools and data streams might we tap to understand the land surface through time? -## Learning Outcomes {.unlisted .unnumbered} +#### Learning Outcomes {.unlisted .unnumbered} * Distinguishing between merging sensor data and merging classifications made from sensors. * Working with the Bayesian Updating of Land Cover (BULC) algorithm, in its basic form, to blend classifications made across multiple years and sensors. * Working with the BULC-D algorithm to highlight locations that changed. -## Assumes you know how to:{.unlisted .unnumbered} +#### Assumes you know how to:{.unlisted .unnumbered} * Import images and image collections, filter, and visualize (Part F1). @@ -2446,7 +2322,7 @@ As the ability to rapidly produce classifications of satellite images grows, it * Obtain accuracy metrics from classifications (Chap. F2.2). ::: -## Introduction {.unlisted .unnumbered} +### Introduction {.unlisted .unnumbered} When working with multiple sensors, we are often presented with a challenge: What to do with classification noise? It’s almost impossible to remove all noise from a classification. Given the information contained in a stream of classifications, however, you should be able to use the temporal context to distinguish noise from true changes in the landscape. @@ -2460,7 +2336,7 @@ If the Events agree strongly, they are evidence of the true condition of the lan BULC’s code is written in JavaScript, with modules that weigh evidence for and against change in several ways, while recording parts of the data-weighing process for you to inspect. In this lab, we will explore BULC through its graphical user interface (GUI), which allows rapid interaction with the algorithm’s main functionality. -## Imagery and Classifications of the Roosevelt River +### Imagery and Classifications of the Roosevelt River How has the Roosevelt River area changed in recent decades? One way to view the area’s recent history is to use Google Earth Timelapse, which shows selected annual clear images of every part of Earth’s terrestrial surface since the 1980s. (You can find the site quickly with a web search.) Enter “Roosevelt River, Brazil” in the search field. For centuries, this area was very remote from agricultural development. It was so little known to Westerners that when former US President Theodore Roosevelt traversed it in the early 1900s there was widespread doubt about whether his near-death experience there was exaggerated or even entirely fictional (Millard 2006). After World War II, the region saw increased agricultural development. Fortin et al. (2020) traced four decades of the history of this region with satellite imagery. Timelapse, meanwhile, indicates that land cover conversion continued after 2016. Can we track it using Earth Engine? @@ -2533,7 +2409,7 @@ In the thumbnail sequence, the color palette shows Forest (class 1) as green, Wa :::{.callout-note} Code Checkpoint F48a. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Basics of the BULC Interface +### Basics of the BULC Interface To see if BULC can successfully sift through these Events, we will use BULC’s GUI (Fig. F4.8.1), which makes interacting with the functionality straightforward. :::{.callout-note} Code Checkpoint F48b in the book’s repository contains information about accessing that interface. @@ -2586,7 +2462,7 @@ In the rightmost panel below the Console, the interface offers you multiple opti Question 2. Select the BULC option, then select the Movie tool to view the result, and choose a drawing speed and resolution. When viewing the full area, would you assess the additional LULC changes since 2016 as being minor, moderate, or major compared to the changes that occurred before 2016? Explain the reasoning for your assessment. -## Detailed LULC Inspection with BULC +### Detailed LULC Inspection with BULC BULC results can be viewed interactively, allowing you to view more detailed estimations of the LULC around the study area. We will zoom into a specific area where change did occur after 2016. To do that, turn on the Satellite view and zoom in. Watching the scale bar in the lower right of the Map panel, continue zooming until the scale bar says 5 km. Then, enter "-60.742, -9.844" in the Earth Engine search tool, located above the code. The text will be interpreted as a longitude/latitude value and will offer you a nearby coordinate, indicated with a value for the degrees West and the degrees South. Click that entry and Earth Engine will move to that location, while keeping at the specified zoom level. Let’s compare the BULC result in this sector against the image from Earth Engine’s satellite view that is underneath it (Fig. 4.8.4). @@ -2606,7 +2482,7 @@ With this finer-scale access to the results of BULC, you can select individual p Question 3. Run the code again with the same data, but adjust the three levelers, then view the results presented in the Map window and the Results panel. How do each of the three parameters affect the behavior of BULC in its results? Use the thumbnail to assess your subjective satisfaction with the results, and use the Inspector to view the BULC behavior in individual pixels. Can you produce an optimal outcome for this given set of input classifications? -## Change Detection with BULC-D +### Change Detection with BULC-D What if we wanted to identify areas of likely change or stability without trying to identify the initial and final LULC class? BULC-D is an algorithm that estimates, at each time step, the probability of noteworthy change. The example below uses the Normalized Burn Ratio (NBR) as a gauge: BULC-D assesses whether the ratio has meaningfully increased, decreased, or remained the same. It is then the choice of the analyst to decide how to treat these assessed probabilities of stability and change. @@ -2660,13 +2536,13 @@ Figure F4.8.8 also shows that, for that pixel, the fit of values for the years u Fig. F4.8.9 shows the NBR history for the green balloon in the southern part of the study area in Fig. F4.8.4. For that pixel, the values in the expectation collection formed a sine wave, and the values in the target collection deviated only slightly from the expectation during the target year. -## Change Detection with BULC and Dynamic World +### Change Detection with BULC and Dynamic World Recent advances in neural networks have made it easier to develop consistent models of LULC characteristics using satellite data. The Dynamic World project (Brown et al. 2022) applies a neural network, trained on a very large number of images, to each new Sentinel-2 image soon after it arrives. The result is a near-real-time classification interpreting the LULC of Earth’s surface, kept continually up to date with new imagery. What to do with the inevitable inconsistencies in a pixel’s stated LULC class through time? For a given pixel on a given image, its assigned class label is chosen by the Dynamic World algorithm as the maximum class probability given the band values on that day. Individual class probabilities are given as part of the dataset and could be used to better interpret a pixel’s condition and perhaps its history. Future work with BULC will involve incorporating these probabilities into BULC’s probability-based structure. For this tutorial, we will explore the consistency of the assigned labels in this same Roosevelt River area as a way to illustrate BULC’s potential for minimizing noise in this vast and growing dataset. -### 5.1. Using BULC To Explore and Refine Dynamic World Classifications +#### 5.1. Using BULC To Explore and Refine Dynamic World Classifications Code Checkpoint A48d. The book’s repository contains a script to use to begin this section. You will need to load the linked script and run it to begin. @@ -2698,7 +2574,7 @@ As BULC uses the classification inputs to estimate the state of the LULC at each ![Fig. F4.8.12 Still frame from the animation of changing confidence through time, near Muiraquitã.](F4/image40.png) -### 5.2. Using BULC To Visualize Uncertainty of Dynamic World in Simplified Categories +#### 5.2. Using BULC To Visualize Uncertainty of Dynamic World in Simplified Categories In the previous section, you may have noticed that there are two main types of uncertainty in BULC’s assessment of long-term classification confidence. One type is due to spatial uncertainty at the edge of two relatively distinct phenomena, like the River/Forest boundary visible in Fig. F4.8.12. These are shown in dark tones in the confidence images, and emphasized in the Probability Hillshade. The other type of uncertainty is due to some cause of labeling uncertainty, due either (1) to the similarity of the classes, or (2) to persistent difficulty in distinguishing two distinct classes that are meaningfully different but spectrally similar. An example of uncertainty due to similar labels is distinguishing flooded and non-flooded wetlands in classifications that contain both those categories. An example of difficulty distinguishing distinct but spectrally similar classes might be distinguishing a parking lot from a body of water. @@ -2723,44 +2599,13 @@ The confidence image shown in the main Map panel is instructive (Fig. 4.8.13). U Given the tools and approaches presented in this lab, you should now be able to import your own classifications for BULC (Sects. 1–3), detect changes in sets of raw imagery (Sect. 4), or use Dynamic World’s pre-created classifications (Sect. 5). The following exercises explore this potential. -## Synthesis {.unnumbered} - -Assignment 1. For a given set of classifications as inputs, BULC uses three parameters that specify how strongly to trust the initial classification, how heavily to weigh the evidence of each classification, and how to adjust the confidence at the end of each time step. For this exercise, adjust the values of these three parameters to explore the strength of the effect they can have on the BULC results. - -Assignment 2. The BULC-D framework produces a continuous three-value vector of the probability of change at each pixel. This variability accounts for the mottled look of the figures when those probabilities are viewed across space. Use the Inspector tool or the interface to explore the final estimated probabilities, both numerically and as represented by different colors of pixels in the given example. Compare and contrast the mean NBR values from the earlier and later years, which are drawn in the Layer list. Then answer the following questions: - -1. In general, how well does BULC-D appear to be identifying locations of likely change? -2. Does one type of change (decrease, increase, no change) appear to be mapped better than the others? If so, why do you think this is? - -Assignment 3. The BULC-D example used here was for 2021. Run it for 2022 or later at this location. How well do results from adjacent years complement each other? - -Assignment 4. Run BULC-D in a different area for a year of interest of your choosing. How do you like the results? - -Assignment 5. Describe how you might use BULC-D as a filter for distinguishing meaningful change from noise. In your answer, you can consider using BULC-D before or after BULC or some other time-series algorithm, like CCDC or LandTrendr. - -Assignment 6. Analyze stability and change with Dynamic World for other parts of the world and for other years. For example, you might consider: - -1. Quebec, Canada, days 150–300 for 2019 and 2020: - -[[[-71.578, 49.755], [-71.578, 49.445], [-70.483, 49.445], [-70.483, 49.755]]] - -Location of a summer 2020 fire - -2. Addis Ababa, Ethiopia: [[[38.79, 9.00], [38.79, 8.99], [38.81, 8.99], [38.81, 9.00]]] - -3. Calacalí, Ecuador: [[[-78.537, 0.017], [-78.537, -0.047], [-78.463, -0.047], [-78.463, 0.017]]] - -4. Irpin, Ukraine: [[[30.22, 50.58], [30.22, 50.525], [30.346, 50.525], [30.346, 50.58]]] - -5. A different location of your own choosing. To do this, use the Earth Engine drawing tools to draw a rectangle somewhere on Earth. Then, at the top of the Import section, you will see an icon that looks like a sheet of paper. Click that icon and look for the polygon specification for the rectangle you drew. Paste that into the location field for the Dynamic World interface. - -## Conclusion {.unnumbered} +### Conclusion {.unnumbered} In this lab, you have viewed several related but distinct ways to use Bayesian statistics to identify locations of LULC change in complex landscapes. While they are standalone algorithms, they are each intended to provide a perspective either on the likelihood of change (BULC-D) or of extracting signal from noisy classifications (BULC). You can consider using them especially when you have pixels that, despite your best efforts, periodically flip back and forth between similar but different classes. BULC can help ignore noise, and BULC-D can help reveal whether this year’s signal has precedent in past years. To learn more about the BULC algorithm, you can view this interactive probability illustration tool by a link found in script F48s1 - Supplemental in the book's repository. In the future, after you have learned how to use the logic of BULC, you might prefer to work with the JavaScript code version. To do that, you can find a tutorial at the website of the authors. -## References {.unnumbered} +### References {.unnumbered} Brown CF, Brumby SP, Guzder-Williams B, et al (2022) Dynamic World, Near real-time global 10 m land use land cover mapping. Sci Data 9:1–17. https://doi.org/10.1038/s41597-022-01307-4 @@ -2776,16 +2621,16 @@ Millard C (2006) The River of Doubt: Theodore Roosevelt’s Darkest Journey. Anc -# Exploring Lagged Effects in Time Series +## Exploring Lagged Effects in Time Series :::{.callout-tip} -# Chapter Information +## Chapter Information -## Author {.unlisted .unnumbered} +#### Author {.unlisted .unnumbered} @@ -2793,12 +2638,12 @@ Andréa Puzzi Nicolau, Karen Dyson, David Saah, Nicholas Clinton -## Overview {.unlisted .unnumbered} +#### Overview {.unlisted .unnumbered} In this chapter, we will introduce lagged effects to build on previous work in modeling time-series data. Time-lagged effects occur when an event at one point in time impacts dependent variables at a later point in time. You will be introduced to concepts of autocovariance and autocorrelation, cross-covariance and cross-correlation, and auto-regressive models. At the end of this chapter, you will be able to examine how variables relate to one another across time, and to fit time series models that take into account lagged events. -## Learning Outcomes {.unlisted .unnumbered} +#### Learning Outcomes {.unlisted .unnumbered} * Using the ee.Join function to create time-lagged collections. @@ -2806,7 +2651,7 @@ In this chapter, we will introduce lagged effects to build on previous work in m * Calculating cross-covariance and cross-correlation. * Fitting auto-regressive models. -## Assumes you know how to:{.unlisted .unnumbered} +#### Assumes you know how to:{.unlisted .unnumbered} * Import images and image collections, filter, and visualize (Part F1). @@ -2817,7 +2662,7 @@ In this chapter, we will introduce lagged effects to build on previous work in m * Fit linear and nonlinear functions with regression in an ImageCollection time series (Chap. F4.6). ::: -## Introduction {.unlisted .unnumbered} +### Introduction {.unlisted .unnumbered} While fitting functions to time series allows you to account for seasonality in your models, sometimes the impact of a seasonal event does not impact your dependent variable until the next month, the next year, or even multiple years later. For example, coconuts take 18–24 months to develop from flower to harvestable size. Heavy rains during the flower development stage can severely reduce the number of coconuts that can be harvested months later, with significant negative economic repercussions. These patterns—where events in one time period impact our variable of interest in later time periods—are important to be able to include in our models. @@ -2825,7 +2670,7 @@ While fitting functions to time series allows you to account for seasonality in In this chapter, we introduce lagged effects into our previous discussions on interpreting time-series data (Chaps. F4.6 and F4.7). Being able to integrate lagged effects into our time-series models allows us to address many important questions. For example, streamflow can be accurately modeled by taking into account previous streamflow, rainfall, and soil moisture; this improved understanding helps predict and mitigate the impacts of drought and flood events made more likely by climate change (Sazib et al. 2020). As another example, time-series lag analysis was able to determine that decreased rainfall was associated with increases in livestock disease outbreaks one year later in India (Karthikeyan et al. 2021). -## Autocovariance and Autocorrelation +### Autocovariance and Autocorrelation Before we dive into autocovariance and autocorrelation, let’s set up an area of interest and dataset that we can use to illustrate these concepts. We will work with a detrended time series (as seen in Chap. F4.6) based on the USGS Landsat 8 Level 2, Collection 2, Tier 1 image collection. Copy and paste the code below to filter the Landsat 8 collection to a point of interest over California and specific dates, and apply the pre-processing function—to mask clouds (as seen in Chap. F4.3) and to scale and add variables of interest (as seen in Chap. F4.6). @@ -2997,7 +2842,7 @@ It's worth noting that you can do this for longer lags as well. Of course, that :::{.callout-note} Code Checkpoint F49a. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Cross-Covariance and Cross-Correlation +### Cross-Covariance and Cross-Correlation Cross-covariance is analogous to autocovariance, except instead of measuring the correspondence between a variable and itself at a lag, it measures the correspondence between a variable and a covariate at a lag. Specifically, we will define the cross-covariance and cross-correlation according to Shumway and Stoffer (2019). @@ -3065,7 +2910,7 @@ As long as there is sufficient temporal overlap between the time series, these t :::{.callout-note} Code Checkpoint F49b. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Auto-Regressive Models +### Auto-Regressive Models The discussion of autocovariance preceded this section in order to introduce the concept of lag. Now that you have a way to get previous values of a variable, it's worth considering auto-regressive models. Suppose that pixel values at time t depend in some way on previous pixel values—auto-regressive models are time series models that use observations from previous time steps as input to a regression equation to predict the value at the next time step. If you have observed significant, non-zero autocorrelations in a time series, this is a good assumption. Specifically, you may postulate a linear model such as the following, where pt is a pixel at time t, and et is a random error (Chap. F4.6): @@ -3143,15 +2988,11 @@ It may be possible to avoid this problem by substituting the output from equatio :::{.callout-note} Code Checkpoint F49d. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Synthesis {.unnumbered} - -Assignment 1. Analyze cross-correlation between NDVI and soil moisture, or precipitation and soil moisture, for example. Earth Engine contains different soil moisture datasets in its catalog (e.g., NASA-USDA SMAP, NASA-GLDAS). Try increasing the lagged time and see if it makes any difference. Alternatively, you can pick any other environmental variable/index (e.g., a different vegetation index: EVI instead of NDVI, for example) and analyze its autocorrelation. - -## Conclusion {.unnumbered} +### Conclusion {.unnumbered} In this chapter, we learned how to use autocovariance and autocorrelation to explore the relationship between elements of a time series at multiple time steps. We also explored how to use cross-covariance and cross-correlation to examine the relationship between elements of two time series at different points in time. Finally, we used auto-regressive models to regress the elements of a time series with elements of the same time series at a different point in time. With these skills, you can now examine how events in one time period impact your variable of interest in later time periods. While we have introduced the linear approach to lagged effects, these ideas can be expanded to more complex models. -## References {.unnumbered} +### References {.unnumbered} Karthikeyan R, Rupner RN, Koti SR, et al (2021) Spatio-temporal and time series analysis of bluetongue outbreaks with environmental factors extracted from Google Earth Engine (GEE) in Andhra Pradesh, India. Transbound Emerg Dis 68:3631–3642. https://doi.org/10.1111/tbed.13972 diff --git a/F5.qmd b/F5.qmd index 3124647..a88f3f4 100644 --- a/F5.qmd +++ b/F5.qmd @@ -1,4 +1,4 @@ -# Vectors and Tables +## Vectors and Tables @@ -7,29 +7,30 @@ In addition to raster data processing, Earth Engine supports a rich set of vecto -# Exploring Vectors +## Exploring Vectors :::{.callout-tip} -# Chapter Information +## Chapter Information -## Author {.unlisted .unnumbered} +#### Author {.unlisted .unnumbered} AJ Purdy, Ellen Brock, David Saah - -Overview +#### Overview {.unlisted .unnumbered} In this chapter, you will learn about features and feature collections and how to use them in conjunction with images and image collections in Earth Engine. Maps are useful for understanding spatial patterns, but scientists often need to extract statistics to answer a question. For example, you may make a false-color composite showing which areas of San Francisco are more “green”—i.e., have more healthy vegetation—than others, but you will likely not be able to directly determine which block in a neighborhood is the most green. This tutorial will demonstrate how to do just that by utilizing vectors. As described in Chap. F4.0, an important way to summarize and simplify data in Earth Engine is through the use of reducers. Reducers operating across space were used in Chap. F3.0, for example, to enable image regression between bands. More generally, chapters in Part F3 and Part F4 used reducers mostly to summarize the values across bands or images on a pixel-by-pixel basis. What if you wanted to summarize information within the confines of given spatial elements- for example, within a set of polygons? In this chapter, we will illustrate and explore Earth Engine’s method for doing that, which is through a reduceRegions call. -Learning Outcomes +#### Learning Outcomes {.unlisted .unnumbered} + + * Uploading and working with a shapefile as an asset to use in Earth Engine. * Creating a new feature using the geometry tools. @@ -38,13 +39,14 @@ Learning Outcomes * Use reduceRegions to summarize an image in irregular neighborhoods. * Exporting calculated data to tables with Tasks. -Assumes you know how to: +### Assumes you know how to: {.unlisted .unnumbered} * Import images and image collections, filter, and visualize (Part F1). * Calculate and interpret vegetation indices (Chap. F2.0). * Use drawing tools to create points, lines, and polygons (Chap. F2.1). +::: -Introduction to Theory +### Introduction {.unlisted .unnumbered} In the world of geographic information systems (GIS), data is typically thought of in one of two basic data structures: raster and vector. In previous chapters, we have principally been focused on raster data—data using the remote sensing vocabulary of pixels, spatial resolution, images, and image collections. Working within the vector framework is also a crucial skill to master. If you don’t know much about GIS, you can find any number of online explainers of the distinctions between these data types, their strengths and limitations, and analyses using both data types. Being able to move fluidly between a raster conception and a vector conception of the world is powerful, and is facilitated with specialized functions and approaches in Earth Engine. @@ -55,7 +57,7 @@ As you have seen, raster features in Earth Engine are stored as an Image or as p In the following example, you will use features and feature collections to identify which city block near the University of San Francisco (USF) campus is the most green. -## Using Geometry Tools to Create Features in Earth Engine +### Using Geometry Tools to Create Features in Earth Engine To demonstrate how geometry tools in Earth Engine work, let’s start by creating a point, and two polygons to represent different elements on the USF campus. @@ -76,7 +78,7 @@ After you create these layers, rename the geometry imports at the top of your sc :::{.callout-note} Code Checkpoint F50a. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Loading Existing Features and Feature Collections in Earth Engine +### Loading Existing Features and Feature Collections in Earth Engine If you wish to have the exact same geometry imports in this chapter for the rest of this exercise, begin this section using the code at the Code Checkpoint above. @@ -92,11 +94,11 @@ Map.addLayer(tiger, { 'color': 'black'}, 'Tiger', false); ``` You should now have the geometry for USF’s campus and a layer added to your map that is not visualized for census blocks across the United States. Next, we will use neighborhood data to spatially filter the TIGER feature collection for blocks near USF’s campus. -## Importing Features into Earth Engine +### Importing Features into Earth Engine There are many image collections loaded in Earth Engine, and they can cover a very large area that you might want to study. Borders can be quite intricate (for example, a detailed coastline), and fortunately there is no need for you to digitize the intricate boundary of a large geographic area. Instead, we will show how to find a spatial dataset online, download the data, and load this into Earth Engine as an asset for use. -### Find a Spatial Dataset of San Francisco Neighborhoods +#### Find a Spatial Dataset of San Francisco Neighborhoods Use your internet searching skills to locate the “Analysis Neighborhoods” dataset covering San Francisco. This data might be located in a number of places, including DataSF, the City of San Francisco’s public-facing data repository. @@ -107,14 +109,14 @@ After you find the Analysis Neighborhoods layer, click Export and select Shapefi Extract the folder to your computer. When you open the folder, you will see that there are actually many files. The extensions (.shp, .dbf, .shx, .prj) all provide a different piece of information to display vector-based data. The .shp file provides data on the geometry. The .dbf file provides data about the attributes. The .shx file is an index file. Lastly, the .prj file describes the map projection of the coordinate information for the shapefile. You will need to load all four files to create a new feature asset in Earth Engine. -### Upload SF Neighborhoods File as an Asset +#### Upload SF Neighborhoods File as an Asset Navigate to the Assets tab (near Scripts). Select New > Table Upload > Shape files (Fig. F5.0.4). ![Fig. F5.0.4 Import an asset as a zipped folder](F5/image52.png) -### Select Files and Name Asset +#### Select Files and Name Asset Click the Select button and then use the file navigator to select the component files of the shapefile structure (i.e., .shp, .dbf, .shx, and .prj) (Fig. F5.0.5). Assign an Asset Name so you can recognize this asset. @@ -140,9 +142,9 @@ Note that if you have any trouble with loading the FeatureCollection using the t :::{.callout-note} Code Checkpoint F50b. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Filtering Feature Collections by Attributes +### Filtering Feature Collections by Attributes -### Filter by Geometry of Another Feature +#### Filter by Geometry of Another Feature First, let’s find the neighborhood associated with USF. Use the first point you created to find the neighborhood that intersects this point; filterBounds is the tool that does that, returning a filtered feature. @@ -159,7 +161,7 @@ var usfTiger = tiger.filterBounds(usfNeighborhood); Map.addLayer(usfTiger, {}, 'usf_Tiger'); ``` -### Filter by Feature (Attribute) Properties +#### Filter by Feature (Attribute) Properties In addition to filtering a FeatureCollection by the location of another feature, you can also filter it by its properties. First, let’s print the usfTiger variable to the Console and inspect the object. @@ -187,7 +189,7 @@ Map.addLayer(housing10_g50_l250, { 'color': 'Magenta'}, 'housing'); We have combined spatial and attribute information to narrow the set to only those blocks that meet our criteria of having between 50 and 250 housing units. -### Print Feature (Attribute) Properties to Console +#### Print Feature (Attribute) Properties to Console We can print out attribute information about these features. The block of code below prints out the area of the resultant geometry in square meters. @@ -208,13 +210,13 @@ Both of the above sections of code provide meaningful information about each fea :::{.callout-note} Code Checkpoint F50c. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Reducing Images Using Feature Geometry +### Reducing Images Using Feature Geometry Now that we have identified the blocks around USF’s campus that have the right housing density, let’s find which blocks are the greenest. The Normalized Difference Vegetation Index (NDVI), presented in detail in Chap. F2.0, is often used to compare the greenness of pixels in different locations. Values on land range from 0 to 1, with values closer to 1 representing healthier and greener vegetation than values near 0. -### Create an NDVI Image +#### Create an NDVI Image The code below imports the Landsat 8 ImageCollection as landsat8. Then, the code filters for images in 2021. Lastly, the code sorts the images from 2021 to find the least cloudy day. @@ -237,7 +239,7 @@ var nir = image.select('B5'); var red = image.select('B4'); var ndvi = nir.subtract(red).divide(nir.add(red)).rename('NDVI'); -### Clip the NDVI Image to the Blocks Near USF +#### Clip the NDVI Image to the Blocks Near USF Next, you will clip the NDVI layer to only show NDVI over USF’s neighborhood. @@ -257,7 +259,7 @@ Map.centerObject(usf_point, 14); The NDVI map for all of San Francisco is interesting, and shows variability across the region. Now, let’s compute mean NDVI values for each block of the city. -### Compute NDVI Statistics by Block +#### Compute NDVI Statistics by Block The code below uses the clipped image ndviUSFblocks and computes the mean NDVI value within each boundary. The scale provides a spatial resolution for the mean values to be computed on. @@ -272,7 +274,7 @@ var ndviPerBlock = ndviUSFblocks.reduceRegions({ ``` Now we’ll use Earth Engine to find out which block is greenest. -### Export Table of NDVI Data by Block from Earth Engine to Google Drive +#### Export Table of NDVI Data by Block from Earth Engine to Google Drive Just as we loaded a feature into Earth Engine, we can export information from Earth Engine. Here, we will export the NDVI data, summarized by block, from Earth Engine to a Google Drive space so that we can interpret it in a program like Google Sheets or Excel. @@ -294,7 +296,7 @@ After you run the task, the file will be saved to your Google Drive. You have no :::{.callout-note} Code Checkpoint F50d. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Identifying the Block in the Neighborhood Surrounding USF with the Highest NDVI +### Identifying the Block in the Neighborhood Surrounding USF with the Highest NDVI You are already familiar with filtering datasets by their attributes. Now you will sort a table and select the first element of the table. @@ -318,34 +320,19 @@ Map.addLayer(GreenHousing, { 'color': 'yellow'}, 'Green Housing!'); :::{.callout-note} Code Checkpoint F50e. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Synthesis {.unnumbered} - -Now it’s your turn to use both feature classes and to reduce data using a geographic boundary. Create a new script for an area of interest and accomplish the following assignments. - -Assignment 1. Create a study area map zoomed to a certain feature class that you made. - -Assignment 2. Filter one feature collection using feature properties. - -Assignment 3. Filter one feature collection based on another feature’s location in space. - -Assignment 4. Reduce one image to the geometry of a feature in some capacity; e.g., extract a mean value or a value at a point. - -## Conclusion {.unnumbered} +### Conclusion {.unnumbered} In this chapter, you learned how to import features into Earth Engine. In Sect. 1, you created new features using the geometry tools and loaded a feature from Earth Engine’s Data Catalog. In Sect. 2, you loaded a shapefile to an Earth Engine asset. In Sect. 3, you filtered feature collections based on their properties and locations. Finally, in Sects. 4 and 5, you used a feature collection to reduce an image, then exported the data from Earth Engine. Now you have all the tools you need to load, filter, and apply features to extract meaningful information from images using vector features in Earth Engine. -# Raster/Vector Conversions - - - +## Raster/Vector Conversions :::{.callout-tip} -# Chapter Information +## Chapter Information -## Author {.unlisted .unnumbered} +#### Author {.unlisted .unnumbered} @@ -353,12 +340,12 @@ Keiko Nomura, Samuel Bowers -## Overview {.unlisted .unnumbered} +#### Overview {.unlisted .unnumbered} The purpose of this chapter is to review methods of converting between raster and vector data formats, and to understand the circumstances in which this is useful. By way of example, this chapter focuses on topographic elevation and forest cover change in Colombia, but note that these are generic methods that can be applied in a wide variety of situations. -## Learning Outcomes {.unlisted .unnumbered} +#### Learning Outcomes {.unlisted .unnumbered} * Understanding raster and vector data in Earth Engine and their differing properties. @@ -366,7 +353,7 @@ The purpose of this chapter is to review methods of converting between raster an * Knowing how and why to convert from vector to raster. * Write a function and map it over a FeatureCollection. -## Assumes you know how to:{.unlisted .unnumbered} +#### Assumes you know how to:{.unlisted .unnumbered} * Import images and image collections, filter, and visualize (Part F1). @@ -378,7 +365,7 @@ The purpose of this chapter is to review methods of converting between raster an * Use reduceRegions to summarize an image in irregular shapes (Chap. F5.0). ::: -## Introduction {.unlisted .unnumbered} +### Introduction {.unlisted .unnumbered} Raster data consists of regularly spaced pixels arranged into rows and columns, familiar as the format of satellite images. Vector data contains geometry features (i.e., points, lines, and polygons) describing locations and areas. Each data format has its advantages, and both will be encountered as part of GIS operations. @@ -388,9 +375,9 @@ Raster and vector data are commonly combined (e.g., extracting image information In this exercise, we’ll use topographic elevation and forest change images in Colombia as well as a protected area feature collection to practice the conversion between raster and vector formats, and to identify situations in which this is worthwhile. -## Raster to Vector Conversion +### Raster to Vector Conversion -### Raster to Polygons +#### Raster to Polygons In this section we will convert an elevation image (raster) to a feature collection (vector). We will start by loading the Global Multi-Resolution Terrain Elevation Data 2010 and the Global Administrative Unit Layers 2015 dataset to focus on Colombia. The elevation image is a raster at 7.5 arc-second spatial resolution containing a continuous measure of elevation in meters in each pixel. @@ -497,7 +484,7 @@ We can see now that the polygons have more distinct shapes with many fewer small ![Fig. F5.1.2 Before (left) and after (right) applying focalMode](F5/image37.png) -### Raster to Points +#### Raster to Points Lastly, we will convert a small part of this elevation image into a point vector dataset. For this exercise, we will use the same example and build on the code from the previous subsection. This might be useful when you want to use geospatial data in a tabular format in combination with other conventional datasets such as economic indicators (Fig. F5.1.3). @@ -508,6 +495,7 @@ Lastly, we will convert a small part of this elevation image into a point vector The easiest way to do this is to use sample while activating the geometries parameter. This will extract the points at the centroid of the elevation pixel. +```js var geometry = ee.Geometry.Polygon([ [-89.553, -0.929], [-89.436, -0.929], @@ -516,7 +504,7 @@ var geometry = ee.Geometry.Polygon([ [-89.553, -0.929] ]); -```js + // To zoom into the area, un-comment and run below // Map.centerObject(geometry,12); Map.addLayer(geometry, {}, 'Areas to extract points'); @@ -563,7 +551,7 @@ Map.addLayer(elevationSamplesStratified, {}, 'Stratified samples'); :::{.callout-note} Code Checkpoint F51a. The book’s repository contains a script that shows what your code should look like at this point. ::: -##3. A More Complex Example +###3. A More Complex Example In this section we’ll use two global datasets, one to represent raster formats and the other vectors: @@ -716,7 +704,7 @@ Map.addLayer(deforestationLargeOutline, { :::{.callout-note} Code Checkpoint F51b. The book’s repository contains a script that shows what your code should look like at this point. ::: -### Raster Properties to Vector Fields +#### Raster Properties to Vector Fields Sometimes we want to extract information from a raster to be included in an existing vector dataset. An example might be estimating a deforestation rate for a set of protected areas. Rather than perform this task on a case-by-case basis, we can attach information generated from an image as a property of a feature. @@ -784,11 +772,11 @@ print(wdpaSubset.reduceColumns({ :::{.callout-note} Code Checkpoint F51c. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Vector-to-Raster Conversion +### Vector-to-Raster Conversion In Sect. 1, we used the protected area feature collection as its original vector format. In this section, we will rasterize the protected area polygons to produce a mask and use this to assess rates of forest change. -### Polygons to a Mask +#### Polygons to a Mask The most common operation to convert from vector to raster is the production of binary image masks, describing whether a pixel intersects a line or falls within a polygon. To convert from vector to a raster mask, we can use the ee.FeatureCollection.reduceToImage method. Let’s continue with our example of the WDPA database and Global Forest Change data from the previous section: @@ -852,7 +840,7 @@ This output can be useful when performing large-scale raster operations, such as :::{.callout-note} Code Checkpoint F51d. The book’s repository contains a script that shows what your code should look like at this point. ::: -### A More Complex Example +#### A More Complex Example The reduceToImage method is not the only way to convert a feature collection to an image. We will create a distance image layer from the boundary of the protected area using distance. For this example, we return to the La Paya protected area explored in Sect. 1. @@ -965,37 +953,18 @@ print('Deforestation within a 1km buffer outside the protected area (ha)', :::{.callout-note} Code Checkpoint F51e. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Synthesis {.unnumbered} - -Question 1. In this lab, we quantified rates of deforestation in La Paya. There is another protected area in the Colombian Amazon named Tinigua. By modifying the existing scripts, determine how the dynamics of forest change in Tinigua compare to those in La Paya with respect to: - -* the number of deforestation events -* the year with the greatest number of change events -* the mean average area of change events -* the total area of loss - -Question 2. In Sect. 1.4, we only considered losses of tree cover, but many protected areas will also have increases in tree cover from regrowth (which is typical of shifting agriculture). Calculate growth in hectares using the Global Forest Change dataset’s gain layer for the six protected areas in Sect. 1.4 by extracting the raster properties and adding them to vector fields. Which has the greatest area of regrowth? Is this likely to be sufficient to balance out the rates of forest loss? Note: The gain layer shows locations where tree cover has increased for the period 2001–2012 (0 = no gain, 1 = tree cover increase), so for comparability use deforestation between the same time period of 2001–2012. - -Question 3. In Sect. 2.2, we considered rates of deforestation in a buffer zone around La Paya. Estimate the deforestation rates inside of La Paya using buffer zones. Is forest loss more common close to the boundary of the reserve? - -Question 4. Sometimes it’s advantageous to perform processing using raster operations, particularly at large scales. It is possible to perform many of the tasks in Sect. 1.3 and 1.4 by first converting the protected area vector to raster, and then using only raster operations. As an example, can you display only deforestation events >10 ha in La Paya using only raster data? (Hint: Consider using ee.Image.connectedPixelCount. You may also want to also look at Sect. 2.1). - -## Conclusion {.unnumbered} +### Conclusion {.unnumbered} In this chapter, you learned how to convert raster to vector and vice versa. More importantly, you now have a better understanding of why and when such conversions are useful. Our examples should give you practical applications and ideas for using these techniques. -# Zonal Statistics +## Zonal Statistics - - - - :::{.callout-tip} -# Chapter Information +## Chapter Information -## Author {.unlisted .unnumbered} +#### Author {.unlisted .unnumbered} @@ -1003,12 +972,12 @@ Sara Winsemius and Justin Braaten -## Overview {.unlisted .unnumbered} +#### Overview {.unlisted .unnumbered} The purpose of this chapter is to extract values from rasters for intersecting points or polygons. We will lay out the process and a function to calculate zonal statistics, which includes optional parameters to modify the function, and then apply the process to three examples using different raster datasets and combinations of parameters. -## Learning Outcomes {.unlisted .unnumbered} +#### Learning Outcomes {.unlisted .unnumbered} * Buffering points as square or circular regions. @@ -1017,7 +986,7 @@ The purpose of this chapter is to extract values from rasters for intersecting p * Exporting computation results to a table. * Copying properties from one image to another. -## Assumes you know how to:{.unlisted .unnumbered} +#### Assumes you know how to:{.unlisted .unnumbered} * Recognize similarities and differences among Landsat 5, 7, and 8 spectral bands (Part F1, Part F2, Part F3). @@ -1028,8 +997,10 @@ The purpose of this chapter is to extract values from rasters for intersecting p * Export calculated data to tables with Tasks (Chap. F5.0). * Understand the differences between raster and vector data (Chap. F5.0, Chap. F5.1). * Write a function and map it over a FeatureCollection (Chap. F5.1). +::: + +### Introduction {.unlisted .unnumbered} -## Introduction to Theory Anyone working with field data collected at plots will likely need to summarize raster-based data associated with those plots. For instance, they need to know the Normalized Difference Vegetation Index (NDVI), precipitation, or elevation for each plot (or surrounding region). Calculating statistics from a raster within given regions is called zonal statistics. Zonal statistics were calculated in Chaps. F5.0 and F5.1 using ee.Image.ReduceRegions. Here, we present a more general approach to calculating zonal statistics with a custom function that works for both ee.Image and ee.ImageCollection objects. In addition to its flexibility, the reduction method used here is less prone to “Computed value is too large” errors that can occur when using ReduceRegions with very large or complex ee.FeatureCollection object inputs. @@ -1040,14 +1011,14 @@ In fieldwork, researchers often work with plots, which are commonly recorded as To choose the size of your neighborhood, you will need to consider your research questions, the spatial resolution of the dataset, the size of your field plot, and the error from your GPS. For example, the raster value extracted for randomly placed 20 m diameter plots would likely merit use of a neighborhood mean when using Sentinel-2 or Landsat 8—at 10 m and 30 m spatial resolution, respectively—while using a thermal band from MODIS (Moderate Resolution Imaging Spectroradiometer) at 1000 m may not. While much of this tutorial is written with plot points and buffers in mind, a polygon asset with predefined regions will serve the same purpose. -## Functions +### Functions Two functions are provided; copy and paste them into your script: * A function to generate circular or square regions from buffered points * A function to extract image pixel neighborhood statistics for a given region -### Function: bufferPoints(radius, bounds) +#### Function: bufferPoints(radius, bounds) Our first function, bufferPoints, returns a function for adding a buffer to points and optionally transforming to rectangular bounds (see Table F5.2.1). @@ -1077,7 +1048,7 @@ function bufferPoints(radius, bounds) { return function(pt) { }; } -### Function: zonalStats(fc, params) +#### Function: zonalStats(fc, params) The second function, zonalStats, reduces images in an ImageCollection by regions defined in a FeatureCollection. Note that reductions can return null statistics that you might want to filter out of the resulting feature collection. Null statistics occur when there are no valid pixels intersecting the region being reduced. This situation can be caused by points that are outside of an image or in regions that are masked for quality or clouds. @@ -1187,7 +1158,7 @@ function zonalStats(ic, fc, params) { // Initialize internal params dictionary }); // Converts the feature collection of feature collections to a single //feature collection. }).flatten(); return results; } -## Point Collection Creation +### Point Collection Creation Below, we create a set of points that form the basis of the zonal statistics calculations. Note that a unique plot_id property is added to each point. A unique plot or point ID is important to include in your vector dataset for future filtering and joining. @@ -1202,7 +1173,7 @@ var pts = ee.FeatureCollection([ ee.Feature(ee.Geometry.Point([-118.6010, 37.0 :::{.callout-note} Code Checkpoint F52a. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Neighborhood Statistic Examples +### Neighborhood Statistic Examples The following examples demonstrate extracting raster neighborhood statistics for the following: @@ -1212,11 +1183,11 @@ The following examples demonstrate extracting raster neighborhood statistics for In each example, the points created in the previous section will be buffered and then used as regions to extract zonal statistics for each image in the image collection. -### Topographic Variables +#### Topographic Variables This example demonstrates how to calculate zonal statistics for a single multiband image. This Digital Elevation Model (DEM) contains a single topographic band representing elevation. -###Buffer the Points +####Buffer the Points Nex, we will apply a 45 m radius buffer to the points defined previously by mapping the bufferPoints function over the feature collection. The radius is set to 45 m to correspond to the 90 m pixel resolution of the DEM. In this case, circles are used instead of squares (set the second argument as false, i.e., do not use bounds). @@ -1225,7 +1196,7 @@ Nex, we will apply a 45 m radius buffer to the points defined previously by mapp var ptsTopo = pts.map(bufferPoints(45, false)); ``` -###Calculate Zonal Statistics +####Calculate Zonal Statistics There are two important things to note about the zonalStats function that this example addresses: @@ -1345,17 +1316,17 @@ slope 29.4 -### MODIS Time Series +#### MODIS Time Series A time series of MODIS eight-day surface reflectance composites demonstrates how to calculate zonal statistics for a multiband ImageCollection that requires no preprocessing, such as cloud masking or computation. Note that there is no built-in function for performing region reductions on ImageCollection objects. The zonalStats function that we are using for reduction is mapping the reduceRegions function over an ImageCollection. -###Buffer the Points +####Buffer the Points In this example, suppose the point collection represents center points for field plots that are 100 m x 100 m, and apply a 50 m radius buffer to the points to match the size of the plot. Since we want zonal statistics for square plots, set the second argument of the bufferPoints function to true, so that the bounds of the buffered points are returned. var ptsModis = pts.map(bufferPoints(50, true)); -###Calculate Zonal Statistic +####Calculate Zonal Statistic Import the MODIS 500 m global eight-day surface reflectance composite collection and filter the collection to include data for July, August, and September from 2015 through 2019. @@ -1383,13 +1354,13 @@ var ptsModisStats = zonalStats(modisCol, ptsModis, params);print('Limited MODIS ``` The result is a feature collection with a feature for all combinations of plots and images. Interpreted as a table, the result has 200 rows (5 plots times 40 images) and as many columns as there are feature properties. Feature properties include those from the plot asset and the image, and any associated non-system image properties. Note that the printed results are limited to the first 50 features for brevity. -### Landsat Time Series +#### Landsat Time Series This example combines Landsat surface reflectance imagery across three instruments: Thematic Mapper (TM) from Landsat 5, Enhanced Thematic Mapper Plus (ETM+) from Landsat 7, and Operational Land Imager (OLI) from Landsat 8. The following section prepares these collections so that band names are consistent and cloud masks are applied. Reflectance among corresponding bands are roughly congruent for the three sensors when using the surface reflectance product; therefore the processing steps that follow do not address inter-sensor harmonization. Review the current literature on inter-sensor harmonization practices if you'd like to apply a correction. -###Prepare the Landsat Image Collection +####Prepare the Landsat Image Collection First, define the function to mask cloud and shadow pixels (See Chap. F4.3 for more detail on cloud masking). @@ -1450,7 +1421,7 @@ Merge the prepared sensor collections. var landsatCol = oliCol.merge(etmCol).merge(tmCol); -###Calculate Zonal Statistics +####Calculate Zonal Statistics Reduce each image in the collection by each plot according to the following parameters. Note that this example defines the imgProps and imgPropsRename parameters to copy over and rename just two selected image properties: Landsat image ID and the satellite that collected the data. It also uses the max reducer, which, as an unweighted reducer, will return the maximum value from pixels that have their centroid within the buffer (see Sect. 4.1 below for more details). @@ -1475,7 +1446,7 @@ print('Limited Landsat zonal stats table', ptsLandsatStats.limit(50)); ``` The result is a feature collection with a feature for all combinations of plots and images. -###Dealing with Large Collections +####Dealing with Large Collections If your browser times out, try exporting the results (as described in Chap. F6.2). It’s likely that point feature collections that cover a large area or contain many points (point-image observations) will need to be exported as a batch task by either exporting the final feature collection as an asset or as a CSV/shapefile/GeoJSON to Google Drive or GCS. @@ -1496,15 +1467,15 @@ Export.table.toDrive({ :::{.callout-note} Code Checkpoint F52b. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Additional Notes +### Additional Notes -### Weighted Versus Unweighted Region Reduction +#### Weighted Versus Unweighted Region Reduction A region used for calculation of zonal statistics often bisects multiple pixels. Should partial pixels be included in zonal statistics? Earth Engine lets you decide by allowing you to define a reducer as either weighted or unweighted (or you can provide per-pixel weight specification as an image band). A weighted reducer will include partial pixels in the zonal statistic calculation by weighting each pixel's contribution according to the fraction of the area intersecting the region. An unweighted reducer, on the other hand, gives equal weight to all pixels whose cell center intersects the region; all other pixels are excluded from calculation of the statistic. For aggregate reducers like ee.Reducer.mean and ee.Reducer.median, the default mode is weighted, while identifier reducers such as ee.Reducer.min and ee.Reducer.max are unweighted. You can adjust the behavior of weighted reducers by calling unweighted on them, as in ee.Reducer.mean.unweighted. You may also specify the weights by modifying the reducer with splitWeights; however, that is beyond the scope of this book. -### Copy Properties to Computed Images +#### Copy Properties to Computed Images Derived, computed images do not retain the properties of their source image, so be sure to copy properties to computed images if you want them included in the region reduction table. For instance, consider the simple computation of unscaling Landsat SR data: @@ -1537,7 +1508,7 @@ print('Selected image properties retained', computedImg ``` Now selected properties are included. Use this technique when returning computed, derived images in a mapped function, and in single-image operations. -### Understanding Which Pixels are Included in Polygon Statistics +#### Understanding Which Pixels are Included in Polygon Statistics If you want to visualize what pixels are included in a polygon for a region reducer, you can adapt the following code to use your own region (by replacing geometry), dataset, desired scale, and CRS parameters. The important part to note is that the image data you are adding to the map is reprojected using the same scale and CRS as that used in your region reduction (see Fig. F5.2.3). @@ -1617,32 +1588,11 @@ Finally, here are some notes on CRS and scale: * Specifying the CRS and scale in both the reduceRegion and addLayer functions allows the map visualization to align with the information printed in the Console. * The Earth Engine default, WGS 84 lat long (EPSG:4326), is a generic CRS that works worldwide. The code above reprojects to EPSG:5070, North American Equal Albers, which is a CRS that preserves area for North American locations. Use the CRS that is best for your use case when adapting this to your own project, or maintain (and specify) the CRS of the image using, for example, crs: 'img.projection().crs()'. -## Synthesis {.unnumbered} - -Question 1. Look at the MODIS example (Sect. 3.2), which uses the median reducer. Try modifying the reducer to be unweighted, either by specifying unweighted or using an identifier reducer like max. What happens, and why? - -Question 2. Calculate zonal statistics for your own buffered points or polygons using a raster and reducer of interest. Be sure to consider the spatial scale of the raster and whether a weighted or unweighted reducer would be more appropriate for your interests. - -If the point or polygon file is stored in a local shapefile or CSV file, first upload the data to your Earth Engine assets. All columns in your vector file, such as the plot name, will be retained through this process. Once you have an Earth Engine table asset ready, import the asset into your script by hovering over the name of the asset and clicking the arrow at the right side, or by calling it in your script with the following code. - -var pts = ee.FeatureCollection('users/yourUsername/yourAsset'); - -If you prefer to define points or polygons dynamically rather than loading an asset, you can add them to your script using the geometry tools. See Chap. F2.1 and F5.0 for more detail on adding and creating vector data. - -Question 3. Try the code from Sect. 4.3 using the MODIS data and the first point from the pts variable. Among other modifications, you will need to create a buffer for the point, take a single MODIS image from the collection, and change visualization parameters. - -* Think about the CRS in the code: The code reprojects to EPSG:5070, but MODIS is collected in the sinusoidal projection SR-ORG:6974. Try that CRS and describe how the image changes. -* Is the count reducer weighted or unweighted? Give an example of a circumstance to use a weighted reducer and an example for an unweighted reducer. Specify the buffer size you would use and the spatial resolution of your dataset. - -Question 4. In the examples above, only a single ee.Reducer is passed to the zonalStats function, which means that only a single statistic is calculated (for example, zonal mean or median or maximum). What if you want multiple statistics—can you alter the code in Sect. 3.1 to (1) make the point buffer 500 instead of 45; (2) add the reducer parameter to the params dictionary; and (3) as its argument, supply a combined ee.Reducer that will calculate minimum, maximum, standard deviation, and mean statistics? - -To achieve this you’ll need to chain several ee.Reducer.combine functions together. Note that if you accept all the individual ee.Reducer and ee.Reducer.combine function defaults, you’ll run into two problems related to reducer weighting differences, and whether or not the image inputs are shared among the combined set of reducers. How can you manipulate the individual ee.Reducer and ee.Reducer.combine functions to achieve the goal of calculating multiple zonal statistics in one call to the zonalStats function? - -## Conclusion {.unnumbered} +### Conclusion {.unnumbered} In this chapter, you used functions containing optional parameters to extract raster values for collocated points. You also learned how to buffer points, and apply weighted and unweighted reducers to get different types of zonal statistics. These functions were applied to three examples that differed by raster dataset, reducer, spatial resolution, and scale. Lastly, you covered related topics like weighting of reducers and buffer visualization. Now you’re ready to apply these ideas to your own work! -## References {.unnumbered} +### References {.unnumbered} Cansler CA, McKenzie D (2012) How robust are burn severity indices when applied in a new region? Evaluation of alternate field-based and remote-sensing methods. Remote Sens 4:456–483. https://doi.org/10.3390/rs4020456 @@ -1650,18 +1600,12 @@ Miller JD, Thode AE (2007) Quantifying burn severity in a heterogeneous landscap -# Advanced Vector Operations +## Advanced Vector Operations - - - - - - :::{.callout-tip} -# Chapter Information +## Chapter Information -## Author {.unlisted .unnumbered} +#### Author {.unlisted .unnumbered} @@ -1669,26 +1613,25 @@ Ujaval Gandhi -## Overview {.unlisted .unnumbered} +#### Overview {.unlisted .unnumbered} This chapter covers advanced techniques for visualizing and analyzing vector data in Earth Engine. There are many ways to visualize feature collections, and you will learn how to pick the appropriate method to create visualizations, such as a choropleth map. We will also cover geoprocessing techniques involving multiple vector layers, such as selecting features in one layer by their proximity to features in another layer and performing spatial joins. -## Learning Outcomes {.unlisted .unnumbered} +#### Learning Outcomes {.unlisted .unnumbered} * Visualizing any vector dataset and creating a thematic map. * Understanding joins in Earth Engine. * Carrying out geoprocessing tasks with vector layers in Earth Engine. -## Assumes you know how to:{.unlisted .unnumbered} - +#### Assumes you know how to:{.unlisted .unnumbered} * Filter a FeatureCollection to obtain a subset (Chap. F5.0, Chap. F5.1). * Write a function and map it over a FeatureCollection (Chap. F5.1, Chap. F5.2). +::: - -## Visualizing Feature Collections +### Visualizing Feature Collections There is a distinct difference between how rasters and vectors are visualized. While images are typically visualized based on pixel values, vector layers use feature properties (i.e., attributes) to create a visualization. Vector layers are rendered on the Map by assigning a value to the red, green, and blue channels for each pixel on the screen based on the geometry and attributes of the features. The functions used for vector data visualization in Earth Engine are listed below in increasing order of complexity. @@ -1699,12 +1642,12 @@ There is a distinct difference between how rasters and vectors are visualized. W In the exercises below, we will learn how to use each of these functions and see how they can generate different types of maps. -### Creating a Choropleth Map +#### Creating a Choropleth Map We will use the TIGER: US Census Blocks layer, which stores census block boundaries and their characteristics within the United States, along with the San Francisco neighborhoods layer from Chap. F5.0 to create a population density map for the city of San Francisco. We start by loading the census blocks and San Francisco neighborhoods layers. We use ee.Filter.bounds to filter the census blocks layer to the San Francisco boundary. - +```js var blocks = ee.FeatureCollection('TIGER/2010/Blocks'); var roads = ee.FeatureCollection('TIGER/2016/Roads'); var sfNeighborhoods = ee.FeatureCollection( 'projects/gee-book/assets/F5-0/SFneighborhoods'); @@ -1712,7 +1655,6 @@ var sfNeighborhoods = ee.FeatureCollection( 'projects/gee-book/assets/F5-0/SFn var geometry = sfNeighborhoods.geometry(); Map.centerObject(geometry); -```js // Filter blocks to the San Francisco boundary. var sfBlocks = blocks.filter(ee.Filter.bounds(geometry)); @@ -1732,7 +1674,13 @@ The census blocks table has a property named 'pop10' containing the population t ```js // Add a pop_density column. -var sfBlocks = sfBlocks.map(function(f) { // Get the polygon area in square miles. var area_sqmi = f.area().divide(2.59e6); var population = f.get('pop10'); // Calculate population density. var density = ee.Number(population).divide(area_sqmi); return f.set({ 'area_sqmi': area_sqmi, 'pop_density': density +var sfBlocks = sfBlocks.map(function(f) { // Get the polygon area in square miles. +var area_sqmi = f.area().divide(2.59e6); +var population = f.get('pop10'); // Calculate population density. +var density = ee.Number(population).divide(area_sqmi); +return f.set({ + 'area_sqmi': area_sqmi, + 'pop_density': density }); }); @@ -1747,8 +1695,6 @@ print(stats); ``` You will see that the population density values have a large range. We also have values that are greater than 100,000, so we need to make sure we select a data type that can store values of this size. We create an empty image and cast it to int32, which is able to hold large integer values. -D - The result is an image with pixel values representing the population density of the polygons. We can now use the standard image visualization method to add this layer to the Map (Fig. F5.3.2). Then, we need to determine minimum and maximum values for the visualization parameters.A reliable technique to produce a good visualization is to find minimum and maximum values that are within one standard deviation. From the statistics that we calculated earlier, we can estimate good minimum and maximum values to be 0 and 50000, respectively. var palette = ['fee5d9', 'fcae91', 'fb6a4a', 'de2d26', 'a50f15']; @@ -1762,7 +1708,7 @@ Map.addLayer(sfBlocksPaint.clip(geometry), visParams, 'Population Density'); ![Fig. F5.3.2 San Francisco population density](F5/image41.png) -### Creating a Categorical Map +#### Creating a Categorical Map Continuing the exploration of styling methods, we will now learn about draw and style. These are the preferred methods of styling for points and line layers. Let’s see how we can visualize the TIGER: US Census Roads layer to create a categorical map. @@ -1923,7 +1869,7 @@ Code Checkpoint F53a. The book’s repository contains a script that shows what ::: Save your script for your own future use, as outlined in Chap. F1.0. Then, refresh the Code Editor to begin with a new script for the next section. -## Joins with Feature Collections +### Joins with Feature Collections Earth Engine was designed as a platform for processing raster data, and that is where it shines. Over the years, it has acquired advanced vector data processing capabilities, and users are now able to carry out complex geoprocessing tasks within Earth Engine. You can leverage the distributed processing power of Earth Engine to process large vector layers in parallel. @@ -1934,7 +1880,7 @@ This section shows how you can do spatial queries and spatial joins using multip Joins are one of the harder skills to master, but doing so will help you perform many complex analysis tasks within Earth Engine. We will go through practical examples that will help you understand these concepts and the workflow better. -### Selecting by Location +#### Selecting by Location In this section, we will learn how to select features from one layer that are within a specified distance from features in another layer. We will continue to work with the San Francisco census blocks and roads datasets from the previous section. We will implement a join to select all blocks in San Francisco that are within 1 km of an interstate highway. @@ -2001,7 +1947,7 @@ Map.addLayer(closeBlocksDrawn, {}, 'Blocks within 1km'); ![Fig. F5.3.6 Selected blocks within 1 km of an interstate highway](F5/image40.png) -### Spatial Joins +#### Spatial Joins A spatial join allows you to query two collections based on the spatial relationship. We will now implement a spatial join to count points in polygons. We will work with a dataset of tree locations in San Francisco and polygons of neighborhoods to produce a CSV file with the total number of trees in each neighborhood. @@ -2094,10 +2040,6 @@ The final result is a CSV file with the neighborhood names and total numbers of :::{.callout-note} Code Checkpoint F53b. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Synthesis {.unnumbered} - -Assignment 1. What join would you use if you wanted to know which neighborhood each tree belongs to? Modify the code above to do a join and post-process the result to add a neighborhood property to each tree point. Export the results as a shapefile. - -## Conclusion {.unnumbered} +### Conclusion {.unnumbered} This chapter covered visualization and analysis using vector data in Earth Engine. You should now understand different functions for FeatureCollection visualization and be able to create thematic maps with vector layers. You also learned techniques for doing spatial queries and spatial joins within Earth Engine. Earth Engine is capable of handling large feature collections and can be effectively used for many spatial analysis tasks. \ No newline at end of file diff --git a/F6.qmd b/F6.qmd index a9a91ee..2136c87 100644 --- a/F6.qmd +++ b/F6.qmd @@ -1,24 +1,14 @@ -# Advanced Topics - - - +## Advanced Topics Although you now know the most basic fundamentals of Earth Engine, there is still much more that can be done. The Part presents some advanced topics that can help expand your skill set for doing larger and more complex projects. These include tools for sharing code among users, scaling up with efficient project design, creating apps for non-expert users, and combining R with other information processing platforms. - - -# Advanced Raster Visualization - - - - - +## Advanced Raster Visualization :::{.callout-tip} -# Chapter Information +## Chapter Information -## Author {.unlisted .unnumbered} +#### Author {.unlisted .unnumbered} @@ -26,12 +16,12 @@ Gennadii Donchyts, Fedor Baart -## Overview {.unlisted .unnumbered} +#### Overview {.unlisted .unnumbered} This chapter should help users of Earth Engine to better understand raster data by applying visualization algorithms such as hillshading, hill shadows, and custom colormaps. We will also learn how image collection datasets can be explored by animating them as well as by annotating with text labels, using, for example, attributes of images or values queried from images. -## Learning Outcomes {.unlisted .unnumbered} +#### Learning Outcomes {.unlisted .unnumbered} * Understanding why perceptually uniform colormaps are better to present data and using them efficiently for raster visualization. @@ -40,7 +30,7 @@ This chapter should help users of Earth Engine to better understand raster data * Animating image collections in multiple ways (animated GIFs, exporting video clips, interactive animations with UI controls). * Adding hillshading and shadows to help visualize raster datasets. -## Assumes you know how to:{.unlisted .unnumbered} +#### Assumes you know how to:{.unlisted .unnumbered} * Import images and image collections, filter, and visualize (Part F1). @@ -48,7 +38,7 @@ This chapter should help users of Earth Engine to better understand raster data * Inspect an Image and an ImageCollection, as well as their properties (Chap. F4.1). ::: -## Introduction {.unlisted .unnumbered} +### Introduction {.unlisted .unnumbered} Visualization is the step to transform data into a visual representation. You make a visualization as soon as you add your first layer to your map in Google Earth Engine. Sometimes you just want to have a first look at a dataset during the exploration phase. But as you move towards the dissemination phase, where you want to spread your results, it is good to think about a more structured approach to visualization. A typical workflow for creating visualization consists of the following steps: @@ -66,7 +56,7 @@ A good standard work on all the choices that one can make while creating a visua In this chapter, we will cover several aspects mentioned in the Grammar of Graphics to convert (raster) data into visual elements. The accurate representation of data is essential in science communication. However, color maps that visually distort data through uneven color gradients or are unreadable to those with color-vision deficiency remain prevalent in science (Crameri, 2020). You will also learn how to add annotation text and symbology, while improving your visualizations by mixing images with hillshading as you explore some of the amazing datasets that have been collected in recent years in Earth Engine. -## Palettes +### Palettes In this section we will explore examples of colormaps to visualize raster data. Colormaps translate values to colors for display on a map. This requires a set of colors (referred to as a “palette” in Earth Engine) and a range of values to map (specified by the min and max values in the visualization parameters). @@ -182,7 +172,7 @@ If you zoom in (F6.0.3) you can see how long cracks have recently appeared near :::{.callout-note} Code Checkpoint F60b. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Remapping and Palettes +### Remapping and Palettes Classified rasters in Earth Engine have metadata attached that can help with analysis and visualization. This includes lists of the names, values, and colors associated with class. These are used as the default color palette for drawing a classification, as seen next. The USGS National Land Cover Database (NLCD) is one such example. Let’s access the NLCD dataset, name it nlcd, and view it (Fig. F6.0.4) with its built-in palette. @@ -259,7 +249,7 @@ Using this remapping approach, we can properly visualize the new color palette ( :::{.callout-note} Code Checkpoint F60c. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Annotations +### Annotations Annotations are the way to visualize data on maps to provide additional information about raster values or any other data relevant to the context. In this case, this additional information is usually shown as geometries, text labels, diagrams, or other visual elements. Some annotations in Earth Engine can be added by making use of the ui portion of the Earth Engine API, resulting in graphical user interface elements such as labels or charts added on top of the map. However, it is frequently useful to render annotations as a part of images, such as by visualizing various image properties or to highlight specific areas. @@ -410,7 +400,7 @@ Code Checkpoint F60f. The book’s repository contains a script that shows what ::: Try commenting that event handler and observe how annotation rendering changes when you zoom in or zoom out. -## Animations +### Animations Visualizing raster images as animations is a useful technique to explore changes in time-dependent datasets, but also, to render short animations to communicate how changing various parameters affects the resulting image—for example, varying thresholds of spectral indices resulting in different binary maps or the changing geometry of vector features. @@ -590,7 +580,7 @@ The main advantage of the interactive animation API is that it provides a way to :::{.callout-note} Code Checkpoint F60i. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Terrain Visualization +### Terrain Visualization This section introduces several raster visualization techniques useful to visualize terrain data such as: @@ -695,21 +685,11 @@ Steps in visualizing a topographic dataset: :::{.callout-note} Code Checkpoint F60j. The book’s repository contains a script that shows what your code should look like at this point. ::: -## Synthesis {.unnumbered} - -To synthesize what you have learned in this chapter, you can do the following assignments. - -Assignment 1. Experiment with different color palettes from the palettes library. Try combining palettes with image opacity (using ee.Image.updateMask call) to visualize different physical features (for example, hot or cold areas using temperature and elevation). - -Assignment 2. Render multiple text annotations when generating animations using image collection. For example, show other image properties in addition to date or image statistics generated using regional reducers for every image. - -Assignment 3. In addition to text annotations, try blending geometry elements (lines, polygons) to highlight specific areas of rendered images. - -## Conclusion {.unnumbered} +### Conclusion {.unnumbered} In this chapter we have learned about several techniques that can greatly improve visualization and analysis of images and image collections. Using predefined palettes can help to better comprehend and communicate Earth observation data, and combining with other visualization techniques such as hillshading and annotations can help to better understand processes studied with Earth Engine. When working with image collections, it is often very helpful to analyze their properties through time by visualizing them as animations. Usually, this step helps to better understand dynamics of the changes that are stored in image collections and to develop a proper algorithm to study these changes. -## References {.unnumbered} +### References {.unnumbered} Burrough PA, McDonnell RA, Lloyd CD (2015) Principles of Geographical Information Systems. Oxford University Press @@ -731,16 +711,16 @@ Wilkinson L (2005) The Grammar of Graphics. Springer Verlag -# Collaborating in Earth Engine with Scripts and Assets +## Collaborating in Earth Engine with Scripts and Assets :::{.callout-tip} -# Chapter Information +## Chapter Information -## Author {.unlisted .unnumbered} +#### Author {.unlisted .unnumbered} @@ -748,12 +728,12 @@ Sabrina H. Szeto -## Overview {.unlisted .unnumbered} +#### Overview {.unlisted .unnumbered} Many users find themselves needing to collaborate with others in Earth Engine at some point. Students may need to work on a group project, people from different organizations might want to collaborate on research together, or people may want to share a script or an asset they created with others. This chapter will show you how to collaborate with others and share your work. -## Learning Outcomes {.unlisted .unnumbered} +#### Learning Outcomes {.unlisted .unnumbered} * Understanding when it is important to share a script or asset. @@ -767,13 +747,13 @@ Many users find themselves needing to collaborate with others in Earth Engine at * Using the require function to load modules. * Creating a script to share as a module. -## Assumes you know how to:{.unlisted .unnumbered} +#### Assumes you know how to:{.unlisted .unnumbered} * Sign up for an Earth Engine account, open the Code Editor, and save your script (Chap. F1.0). ::: -## Introduction {.unlisted .unnumbered} +### Introduction {.unlisted .unnumbered} Many people find themselves needing to share a script when they encounter a problem; they wish to share the script with someone else so they can ask a question. When this occurs, sharing a link to the script often suffices. The other person can then make comments or changes before sending a new link to the modified script. @@ -787,7 +767,7 @@ If you or your group members find yourselves repeatedly reusing certain function Let’s get started. For this lab, you will need to work in small groups or pairs. -## Using Get Link to Share a Script +### Using Get Link to Share a Script Copy and paste the following code into the Code Editor. @@ -809,7 +789,7 @@ Question 2. What happens when you check the box for Hide code panel or Disable a Answer. Hide code panel will minimize the code panel so the person you send the script to will see the Map maximized. This is useful when you want to draw the person’s attention to the results rather than to the code. To expand the code panel, they have to click on the Show code button. Disable auto-run is helpful when you do not want the script to start running when the person you sent it to opens it. Perhaps your script takes very long to run or requires particular user inputs and you just want to share the code with the person. -## Sharing Assets from Your Asset Manager +### Sharing Assets from Your Asset Manager When you clicked the Get Link button earlier, you may have noticed a note in the popup reading: “To give others access to assets in the code snapshot, you may need to share them.” If your script uses an asset that you have uploaded into your Asset Manager, you will need to share that asset as well. If not, an error message will appear when the person you shared the script with tries to run it. @@ -841,7 +821,7 @@ Question 4. Share an asset with a group member and give them writer access. Send Answer: You can view details about the asset and import it for use in a script in the Code Editor. You can also share the asset with others and delete the asset. -## Working with Shared Repositories +### Working with Shared Repositories Now that you know how to share assets and scripts, let’s move on to sharing repositories. In this section, you will learn about different types of repositories and how to add a repository that someone else shared with you. You will also learn how to view previous versions of a script and how to revert back to an earlier version. @@ -937,7 +917,7 @@ Answer: The script reverts to the previous version, in which the only line of co print('The owner of this repository is GroupMemberName.'); -## Using the Require Function to Load a Module +### Using the Require Function to Load a Module In earlier chapters, you may have noticed that the require function allows you to reuse code that has already been written without having to copy and paste it into your current script. For example, you might have written a function for cloud masking that you would like to use in multiple scripts. Saving this function as a module enables you to share the code across your own scripts and with other people. Or you might discover a new module with capabilities you need written by other authors. This section will show you how to use the require function to create and share your own module or to load a module that someone else has shared. @@ -1093,15 +1073,11 @@ Map.addLayer(composite, { }); ``` -## Synthesis {.unnumbered} - -Apply what you learned in this chapter by setting up a shared repository for your project, lab group, or organization. What scripts would you share? What permissions should different users have? Are there any scripts that you would turn into modules? - -## Conclusion {.unnumbered} +### Conclusion {.unnumbered} In this chapter, you learned how to collaborate with others in the Earth Engine Code Editor through sharing scripts, assets, and repositories. You learned about different roles and permissions available for sharing and when it is appropriate to use each. In addition, you are now able to see what changes have been made to a script and revert to a previous version. Lastly, you loaded and used a module that was shared with you and created your own module for sharing. You are now ready to start collaborating and developing scripts with others. -## References {.unnumbered} +### References {.unnumbered} Chang A (2017) Making it easier to reuse code with Earth Engine script modules. In: Google Earth and Earth Engine. https://medium.com/google-earth/making-it-easier-to-reuse-code-with-earth-engine-script-modules-2e93f49abb13. Accessed 24 Feb 2022 @@ -1109,16 +1085,16 @@ Donchyts G, Baart F, Braaten J (2019) ee-palettes. https://github.com/gee-commun -# Scaling Up in Earth Engine +## Scaling Up in Earth Engine :::{.callout-tip} -# Chapter Information +## Chapter Information -## Author {.unlisted .unnumbered} +#### Author {.unlisted .unnumbered} @@ -1126,12 +1102,12 @@ Jillian M. Deines, Stefania Di Tommaso, Nicholas Clinton, Noel Gorelick -## Overview {.unlisted .unnumbered} +#### Overview {.unlisted .unnumbered} Commonly, when Earth Engine users move from tutorials to developing their own processing scripts, they encounter the dreaded error messages, “computation timed out” or “user memory limit exceeded.” Computational resources are never unlimited, and the team at Earth Engine has designed a robust system with built-in checks to ensure server capacity is available to everyone. This chapter will introduce general tips for creating efficient Earth Engine workflows that accomplish users’ ambitious research objectives within the constraints of the Earth Engine ecosystem. We use two example case studies: 1) extracting a daily climate time series for many locations across two decades, and 2) generating a regional, cloud-free median composite from Sentinel-2 imagery. -## Learning Outcomes {.unlisted .unnumbered} +#### Learning Outcomes {.unlisted .unnumbered} * Understanding constraints on Earth Engine resource use. @@ -1139,7 +1115,7 @@ Commonly, when Earth Engine users move from tutorials to developing their own pr * Managing large projects and multistage workflows. * Recognizing when using the Python API may be advantageous to execute large batches of tasks. -## Assumes you know how to:{.unlisted .unnumbered} +#### Assumes you know how to:{.unlisted .unnumbered} * Import images and image collections, filter, and visualize (Part F1). @@ -1149,7 +1125,7 @@ Commonly, when Earth Engine users move from tutorials to developing their own pr * Use the require function to load code from existing modules (Chap. F6.1). ::: -## Introduction {.unlisted .unnumbered} +### Introduction {.unlisted .unnumbered} Parts F1–F5 of this book have covered key remote sensing concepts and demonstrated how to implement them in Earth Engine. Most exercises have used local-scale examples to enhance understanding and complete tasks within a class-length time period. But Earth Engine’s power comes from its scalability—the ability to apply geospatial processing across large areas and many years. @@ -1161,7 +1137,7 @@ How we go from small to large scales is influenced by Earth Engine’s design. E 3. Breaking up jobs across space. 4. Building a multipart workflow and exporting intermediate assets. -### Earth Engine: Under the Hood +#### Earth Engine: Under the Hood As you use Earth Engine, you may begin to have questions about how it works and how you can use that knowledge to optimize your workflow. In general, the inner workings are opaque to users. Typical fixes and approaches that data scientists use to manage memory constraints often don’t apply. It’s helpful to know what users can and cannot control, and how your scripts translate to Earth Engine’s server operations. @@ -1221,7 +1197,7 @@ Table assets Maximum of 100 million features, 1000 properties (columns), and 100,000 vertices for a geometry. -### The Importance of Coding Best Practices +#### The Importance of Coding Best Practices Good code scales better than bad code. But what is good code? Generally, for Earth Engine, good code means 1) using Earth Engine’s server-side operators; 2) avoiding multiple passes through the same image collection; 3) avoiding unnecessary conversions; and 4) setting the processing scale or sample numbers appropriate for your use case, i.e., avoid using very fine scales or large samples without reason. @@ -1230,11 +1206,11 @@ We encourage readers to become familiar with the “Coding Best Practices” pag In addition, some Earth Engine functions are more efficient than others. For example, Image.reduceRegions is more efficient than Image.sampleRegions, because sampleRegions regenerates the geometries under the hood. These types of best practices are trickier to enumerate and somewhat idiosyncratic. We encourage users to learn about and make use of the Profiler tab, which will track and display the resources used for each operation within your script. This can help identify areas to focus efficiency improvements. Note that the profiler itself increases resource use, so only use it when necessary to develop a script and remove it for production-level execution. Other ways to discover best practices include following/posting questions to GIS StackExchange or the Earth Engine Developer’s Discussion Group, swapping code with others, and experimentation. -## Scaling Across Time +### Scaling Across Time In this section we use an example of extracting climate data at features (points or polygons) to demonstrate how to scale an operation across many features (Sect. 1.1) and how to break up large jobs by time units when necessary (e.g, by years; Sect. 1.2). -### 1.1. Scaling Up with Earth Engine Operators: Annual Daily Climate Data +#### 1.1. Scaling Up with Earth Engine Operators: Annual Daily Climate Data Earth Engine’s operators are designed to parallelize queries on the backend without user intervention. In many cases, they are sufficient to accomplish a scaling operation. @@ -1355,7 +1331,7 @@ Perhaps even more important if you are seeking to scale up an analysis, includin For example, when we ran the same job as above but did not use the selectors argument, the output dataset was 5.7 GB (versus 6.8 MB!) and the runtime was slower. This is a cumbersomely large file, with no real benefit. We generally recommend dropping the .geo column and other unnecessary properties. To retain spatial information, a unique identifier for each feature can be used for downstream joins with the spatial data or other properties. If working with point data, latitude and longitude columns can be added prior to export to maintain easily accessible geographic information, although the .geo column for point data is far smaller than for irregularly shaped polygon features. -### 1.2. Scaling Across Time by Batching: Get 20 Years of Daily Climate Data +#### 1.2. Scaling Across Time by Batching: Get 20 Years of Daily Climate Data Above, we extracted one year of daily data for our 293 counties. Let’s say we want to do the same thing, but for 2001–2020. We have already written our script to flexibly specify years, so it’s fairly adaptable to this new use case: @@ -1447,7 +1423,7 @@ Note: If at any time you have submitted several tasks to the server but want to In order to auto-execute jobs in batch mode, we’d need to use the Python API. Interested users can see the Earth Engine User Guide Python API tutorial for further details about the Python API. -## Scaling Across Space via Spatial Tiling +### Scaling Across Space via Spatial Tiling Breaking up jobs in space is another key strategy for scaling operations in Earth Engine. Here, we will focus on making a cloud-free composite from the Sentinel-2 Level 2A Surface Reflectance product. The approach is similar to that in Chap. F4.3, which explores cloud-free compositing. The main difference is that Landsat scenes come with a reliable quality band for each scene, whereas the process for Sentinel-2 is a bit more complicated and computationally intense (see below). @@ -1467,7 +1443,7 @@ These three image collections all contain 10 m resolution data for every Sentine To sum up, the cloud masking approach is computationally costly, thus requiring some thought when applying it at scale. -### 2.1. Generate a Cloud-Free Satellite Composite: Limits to On-the-Fly Computing +#### 2.1. Generate a Cloud-Free Satellite Composite: Limits to On-the-Fly Computing Note: Our focus here is on code structure for implementing spatial tiling. Below, we import existing tested functions for cloud masking using the require command. @@ -1559,7 +1535,7 @@ Map.addLayer(median, viz, 'median'); As you can see, this is an excellent candidate for an export task rather than running in “on-the-fly” interactive mode, as above. -### 2.2. Generate a Regional Composite Through Spatial Tiling +#### 2.2. Generate a Regional Composite Through Spatial Tiling Our goal is to apply the cloud masking method in Sect. 2.1 to the state of Washington, United States. In our testing, we successfully exported one Sentinel-2 composite for this area in about nine hours, but for this tutorial, let’s presume we need to split the task up to be successful. @@ -1593,8 +1569,8 @@ Map.addLayer(roi, {}, 'roi'); Map.addLayer(grid, {}, 'grid'); ``` -![Fig. F6.2.7 Visualization of the regular spatial grid generated for use in spatial batch processing](F6/image12.png) +![Fig. F6.2.7 Visualization of the regular spatial grid generated for use in spatial batch processing](F6/image12.png) Next, create a new, empty ImageCollection asset to use as our export destination (Assets > New > Image Collection; Fig. F6.2.8). Name the image collection 'S2_composite_WA' and specify the asset location in your user folder (e.g., "path/to/your/asset/s2_composite_WA"). @@ -1675,7 +1651,7 @@ Code Checkpoint F62e. The book’s repository contains a script that shows what Note the ease, speed, and joy of panning and zooming to explore the pre-computed composite asset (Fig. F6.2.10) compared to the on-the-fly version discussed in Sect. 2.1. -## Multistep Workflows and Intermediate Assets +### Multistep Workflows and Intermediate Assets Often, our goals require several processing steps that cannot be completed within one Earth Engine computational chain. In these cases, the best strategy becomes breaking down tasks into individual pieces that are created, stored in assets, and used across several scripts. Each sequential script creates an intermediate output, and this intermediate output becomes the input to the next script. @@ -1734,463 +1710,11 @@ Tip 4. Consider Google Cloud Platform for hosting larger intermediates If you are working with very large or very many files, you can link Earth Engine with Cloud Projects on Google Cloud Platform. See the Earth Engine documentation on “Setting Up Earth Engine Enabled Cloud Projects” for more information. -## Synthesis {.unnumbered} - -Earth Engine is built to be scaled. Scaling up working scripts, however, can present challenges when the computations take too long or return results that are too large or numerous. We have covered some key strategies to use when you encounter memory or computational limits. Generally, they involve 1) optimizing your code based on Earth Engine’s functions and infrastructure; 2) working at scales appropriate for your data, question, and region of interest, and not at higher resolutions than necessary; and 3) breaking up tasks into discrete units. - -## References {.unnumbered} - -Abatzoglou JT (2013) Development of gridded surface meteorological data for ecological applications and modelling. Int J Climatol 33:121–131. https://doi.org/10.1002/joc.3413 - -Frantz D, Haß E, Uhl A, et al (2018) Improvement of the Fmask algorithm for Sentinel-2 images: Separating clouds from bright surfaces based on parallax effects. Remote Sens Environ 215:471–481. https://doi.org/10.1016/j.rse.2018.04.046 - -Gorelick N, Hancher M, Dixon M, et al (2017) Google Earth Engine: Planetary-scale geospatial analysis for everyone. Remote Sens Environ 202:18–27. https://doi.org/10.1016/j.rse.2017.06.031 - -Wilson G, Bryan J, Cranston K, et al (2017) Good enough practices in scientific computing. PLoS Comput Biol 13:e1005510. https://doi.org/10.1371/journal.pcbi.1005510 - - - -# Sharing Work in Earth Engine: Basic UI and Apps - - - - - -:::{.callout-tip} -# Chapter Information - -## Author {.unlisted .unnumbered} - - - -Qiusheng Wu - - - -## Overview {.unlisted .unnumbered} - - -The purpose of this chapter is to demonstrate how to design and publish Earth Engine Apps using both JavaScript and Python. You will be introduced to the Earth Engine User Interface JavaScript API and the geemap Python package. Upon completion of this chapter, you will be able to publish an Earth Engine App with a split-panel map for visualizing land cover change. - -## Learning Outcomes {.unlisted .unnumbered} - - -* Designing a user interface for an Earth Engine App using JavaScript. -* Publishing an Earth Engine App for visualizing land cover change. -* Developing an Earth Engine App using Python and geemap. -* Deploying an Earth Engine App using a local computer as a web server. -* Publishing an Earth Engine App using Python and free cloud platforms. -* Creating a conda environment using Anaconda/Miniconda. -* Installing Python packages and using Jupyter Notebook. -* Commiting changes to a GitHub repository. - -## Assumes you know how to:{.unlisted .unnumbered} - - -* Import images and image collections, filter, and visualize (Part F1). -* Use the basic functions and logic of Python. - -::: -## Introduction {.unlisted .unnumbered} - - -Earth Engine has a user interface API that allows users to build and publish interactive web apps directly from the JavaScript Code Editor. Many readers will have encountered a call to ui.Chart in other chapters, but much more interface functionality is available. In particular, users can utilize the ui functions to construct an entire graphical user interface (GUI) for their Earth Engine script. The GUI may include simple widgets (e.g., labels, buttons, checkboxes, sliders, text boxes) as well as more complex widgets (e.g., charts, maps, panels) for controlling the GUI layout. A complete list of the ui widgets and more information about panels can be found at the links below. Once a GUI is constructed, users can publish the App from the JavaScript Code Editor by clicking the Apps button above the script panel in the Code Editor. - -* Widgets: [https://developers.google.com/earth-engine/guides/ui_widgets](https://www.google.com/url?q=https://developers.google.com/earth-engine/guides/ui_widgets&sa=D&source=editors&ust=1671458841273029&usg=AOvVaw10aLP4KU7kHJTwcnM5Pr4-) -* Panels: [https://developers.google.com/earth-engine/guides/ui_panels](https://www.google.com/url?q=https://developers.google.com/earth-engine/guides/ui_panels&sa=D&source=editors&ust=1671458841273501&usg=AOvVaw1XwxHFyCULBnago_-VpDUq) - -Unlike the Earth Engine JavaScript API, the Earth Engine Python API does not provide functionality for building interactive user interfaces. Fortunately, the Jupyter ecosystem has ipywidgets, an architecture for creating interactive user interface controls (e.g., buttons, sliders, checkboxes, text boxes, dropdown lists) in Jupyter notebooks that communicate with Python code. The integration of graphical widgets into the Jupyter Notebook workflow allows users to configure ad hoc control panels to interactively sweep over parameters using graphical widget controls. One very powerful widget is the output widget, which can be used to display rich output generated by IPython, such as text, images, charts, and videos. A complete list of widgets and more information about the output widget can be found at the links below. By integrating ipyleaflet (for creating interactive maps) and ipywidgets (for designing interactive user interfaces), the geemap Python package ([https://geemap.org](https://www.google.com/url?q=https://geemap.org&sa=D&source=editors&ust=1671458841274245&usg=AOvVaw3rVUS8lKoS9clb8r42oxe0)) makes it much easier to explore and analyze massive Earth Engine datasets via a web browser in a Jupyter environment suitable for interactive exploration, teaching, and sharing. Users can build interactive Earth Engine Apps using geemap with minimal coding (Fig. F6.3.1). - -* Widgets: [https://ipywidgets.readthedocs.io/en/latest/examples/Widget%20List.html](https://www.google.com/url?q=https://ipywidgets.readthedocs.io/en/latest/examples/Widget%2520List.html&sa=D&source=editors&ust=1671458841274780&usg=AOvVaw2SzMGW9F4r3E-llbiCksCA) -* Output: [https://ipywidgets.readthedocs.io/en/latest/examples/Output%20Widget.html](https://www.google.com/url?q=https://ipywidgets.readthedocs.io/en/latest/examples/Output%2520Widget.html&sa=D&source=editors&ust=1671458841275255&usg=AOvVaw13a3qsly2f0mg0egoTAMqw) - -![Fig. F6.3.1 The GUI of geemap in a Jupyter environment](F6/image1.png) - - - -## Building an Earth Engine App Using JavaScript - -In this section, you will learn how to design a user interface for an Earth Engine App using JavaScript and the Earth Engine User Interface API. Upon completion of this section, you will have an Earth Engine App with a split-panel map for visualizing land cover change using the Landsat-based United States Geological Survey National Land Cover Database (NLCD). - -First, let’s define a function for filtering the NLCD ImageCollection by year and select the landcover band. The function returns an Earth Engine ui.Map.Layer of the landcover band of the selected NLCD image. Note that as of this writing, NLCD spans nine epochs: 1992, 2001, 2004, 2006, 2008, 2011, 2013, 2016, and 2019. The 1992 data are primarily based on unsupervised classification of Landsat data, while the rest of the images rely on the imperviousness data layer for the urban classes and on a decision-tree classification for the rest. The 1992 image is not directly comparable to any later editions of NLCD (see the Earth Engine Data Catalog for more details, if needed). Therefore, we will use only the eight epochs after 2000 in this lab. - -```js -// Get an NLCD image by year. -var getNLCD = function(year) { // Import the NLCD collection. var dataset = ee.ImageCollection( 'USGS/NLCD_RELEASES/2019_REL/NLCD'); // Filter the collection by year. var nlcd = dataset.filter(ee.Filter.eq('system:index', year)) - .first(); // Select the land cover band. var landcover = nlcd.select('landcover'); return ui.Map.Layer(landcover, {}, year); -}; - -``` -Our intention is to create a dropdown list so that when a particular epoch is selected, the corresponding NLCD image layer will be displayed on the map. We’ll define a dictionary with each NLCD epoch as the key and its corresponding NLCD image layer as the value. The keys of the dictionary (i.e., the eight NLCD epochs) will be used as the input to the dropdown lists (ui.Select) on the split-level map. - -```js -// Create a dictionary with each year as the key -// and its corresponding NLCD image layer as the value. -var images = { '2001': getNLCD('2001'), '2004': getNLCD('2004'), '2006': getNLCD('2006'), '2008': getNLCD('2008'), '2011': getNLCD('2011'), '2013': getNLCD('2013'), '2016': getNLCD('2016'), '2019': getNLCD('2019'), -}; - -``` -The split-panel map is composed of two individual maps, leftMap and rightMap. The map controls (e.g., zoomControl, scaleControl, mapTypeControl) will be shown only on rightMap. A control panel (ui.Panel) composed of a label (ui.Label) and a dropdown list (ui.Select) is added to each map. When an NLCD epoch is selected from a dropdown list, the function updateMap will be called to show the corresponding image layer of the selected epoch. - -```js -// Create the left map, and have it display the first layer. -var leftMap = ui.Map(); -leftMap.setControlVisibility(false); -var leftSelector = addLayerSelector(leftMap, 0, 'top-left'); - -// Create the right map, and have it display the last layer. -var rightMap = ui.Map(); -rightMap.setControlVisibility(true); -var rightSelector = addLayerSelector(rightMap, 7, 'top-right'); - -// Adds a layer selection widget to the given map, to allow users to -// change which image is displayed in the associated map. -function addLayerSelector(mapToChange, defaultValue, position) { var label = ui.Label('Select a year:'); // This function changes the given map to show the selected image. function updateMap(selection) { - mapToChange.layers().set(0, images[selection]); - } // Configure a selection dropdown to allow the user to choose // between images, and set the map to update when a user // makes a selection. var select = ui.Select({ - items: Object.keys(images), - onChange: updateMap - }); - select.setValue(Object.keys(images)[defaultValue], true); var controlPanel = ui.Panel({ - widgets: [label, select], - style: { - position: position - } - }); - - mapToChange.add(controlPanel); -} - -When displaying a land cover classification image on the Map, it would be useful to add a legend to make it easier for users to interpret the land cover type associated with each color. Let’s define a dictionary that will be used to construct the legend. The dictionary contains two keys: names (a list of land cover types) and colors (a list of colors associated with each land cover type). The legend will be placed in the bottom right of the Map. - -// Set the legend title. -var title = 'NLCD Land Cover Classification'; - -// Set the legend position. -var position = 'bottom-right'; - -// Define a dictionary that will be used to make a legend -var dict = { 'names': [ '11 Open Water', '12 Perennial Ice/Snow', '21 Developed, Open Space', '22 Developed, Low Intensity', '23 Developed, Medium Intensity', '24 Developed, High Intensity', '31 Barren Land (Rock/Sand/Clay)', '41 Deciduous Forest', '42 Evergreen Forest', '43 Mixed Forest', '51 Dwarf Scrub', '52 Shrub/Scrub', '71 Grassland/Herbaceous', '72 Sedge/Herbaceous', '73 Lichens', '74 Moss', '81 Pasture/Hay', '82 Cultivated Crops', '90 Woody Wetlands', '95 Emergent Herbaceous Wetlands', - ], 'colors': [ '#466b9f', '#d1def8', '#dec5c5', '#d99282', '#eb0000', '#ab0000', '#b3ac9f', '#68ab5f', '#1c5f2c', '#b5c58f', '#af963c', '#ccb879', '#dfdfc2', '#d1d182', '#a3cc51', '#82ba9e', '#dcd939', '#ab6c28', '#b8d9eb', '#6c9fb8', - ] -}; - -``` -With the legend dictionary defined above, we can now create a panel to hold the legend widget and add it to the Map. Each row on the legend widget is composed of a color box followed by its corresponding land cover type. - -```js -// Create a panel to hold the legend widget. -var legend = ui.Panel({ - style: { - position: position, - padding: '8px 15px' } -});// Function to generate the legend. -function addCategoricalLegend(panel, dict, title) { // Create and add the legend title. var legendTitle = ui.Label({ - value: title, - style: { - fontWeight: 'bold', - fontSize: '18px', - margin: '0 0 4px 0', - padding: '0' } - }); - panel.add(legendTitle); var loading = ui.Label('Loading legend...', { - margin: '2px 0 4px 0' }); - panel.add(loading); // Creates and styles 1 row of the legend. var makeRow = function(color, name) { // Create the label that is actually the colored box. var colorBox = ui.Label({ - style: { - backgroundColor: color, // Use padding to give the box height and width. padding: '8px', - margin: '0 0 4px 0' } - }); // Create the label filled with the description text. var description = ui.Label({ - value: name, - style: { - margin: '0 0 4px 6px' } - }); return ui.Panel({ - widgets: [colorBox, description], - layout: ui.Panel.Layout.Flow('horizontal') - }); - }; // Get the list of palette colors and class names from the image. var palette = dict.colors; var names = dict.names; - loading.style().set('shown', false); for (var i = 0; i < names.length; i++) { - panel.add(makeRow(palette[i], names[i])); - } - - rightMap.add(panel); - -} - -The last step is to create a split-panel map to hold the linked maps (leftMap and rightMap) and tie everything together. When users pan and zoom one map, the other map will also be panned and zoomed to the same extent automatically. When users select a year from a dropdown list, the image layer will be updated accordingly. Users can use the slider to swipe through and visualize land cover change easily (Fig. F6.3.2). Please make sure you minimize the Code Editor and maximize the Map so that you can see the dropdown widget in the upper-right corner of the map. - -addCategoricalLegend(legend, dict, title); - -// Create a SplitPanel to hold the adjacent, linked maps. -var splitPanel = ui.SplitPanel({ - firstPanel: leftMap, - secondPanel: rightMap, - wipe: true, - style: { - stretch: 'both' } -});// Set the SplitPanel as the only thing in the UI root. -ui.root.widgets().reset([splitPanel]); -var linker = ui.Map.Linker([leftMap, rightMap]); -leftMap.setCenter(-100, 40, 4); - -``` -:::{.callout-note} -Code Checkpoint F63a. The book’s repository contains a script that shows what your code should look like at this point. -::: -![Fig. F6.3.2 A split-panel map for visualizing land cover change using NLCD](F6/image32.png) - - -## Publishing an Earth Engine App from the Code Editor - -The goal of this section is to publish the Earth Engine App that we created in Sect. 1. The look and feel of interfaces changes often; if the exact windows described below change over time, the concepts should remain stable to help you to publish your App. First, load the script (see the Code Checkpoint in Sect. 1) into the Code Editor. Then, open the Manage Apps panel by clicking the Apps button above the script panel in the Code Editor (Fig. F6.3.3). - -![Fig. F6.3.3 The Apps button in the JavaScript Code Editor](F6/image14.png) - - -Now click on the New App button (Fig. F6.3.4). - -![Fig. F6.3.4 The New App button](F6/image61.png) - - -In the Publish New App dialog (Fig. F6.3.5), choose a name for the App (e.g., NLCD Land Cover Change), select a Google Cloud Project, provide a thumbnail to be shown in the Public Apps Gallery, and specify the location of the App’s source code. You may restrict access to the App to a particular Google Group or make it publicly accessible. Check Feature this app in your Public Apps Gallery if you would like this App to appear in your public gallery of Apps available at https://USERNAME.users.earthengine.app. When all fields are filled out and validated, the Publish button will be enabled; click it to complete publishing the App. - -![Fig. F6.3.5 The Publish New App dialog](F6/image37.png) - - -To manage an App from the Code Editor, open the Manage Apps panel (Fig. F6.3.6) by clicking the Apps button above the script panel in the Code Editor (Fig. F6.3.3). There, you can update your App’s configuration or delete the App. - -![Fig. F6.3.6 The Manage Apps panel](F6/image65.png) - - -## Developing an Earth Engine App Using geemap - -In this section, you will learn how to develop an Earth Engine App using the geemap Python package and Jupyter Notebook. The geemap package is available on both [PyPI](https://www.google.com/url?q=https://pypi.org/project/leafmap&sa=D&source=editors&ust=1671458841303040&usg=AOvVaw25Lg4yxDI4b31VJovOcUOc) (pip) and [conda-forge](https://www.google.com/url?q=https://anaconda.org/conda-forge/leafmap&sa=D&source=editors&ust=1671458841303462&usg=AOvVaw267tbdQ3NyyHL53Jmk4W4k). It is highly recommended that you create a fresh conda environment to install geemap. - -:::{.callout-note} -Code Checkpoint F63b. The book’s repository contains information about setting up a conda environment and installing geemap. -::: -Once you have launched a Jupyter notebook in your browser, you can continue with the next steps of the lab. - -On the Jupyter Notebook interface, click the New button in the upper-right corner and select Python 3 to create a new notebook. Change the name of the notebook from “Untitled” to something meaningful (e.g., “nlcd_app”). With the newly created notebook, we can start writing and executing Python code. - -First, let’s import the Earth Engine and geemap libraries. Press Alt + Enter to execute the code and create a new cell below. - -import eeimport geemap - -Create an interactive map by specifying the map center (latitude, longitude) and zoom level (1–18). If this is the first time you use geemap and the Earth Engine Python API, a new tab will open in your browser asking you to authenticate Earth Engine. Follow the on-screen instructions to authenticate Earth Engine. - -Map = geemap.Map(center=[40, -100], zoom=4) -Map - -Retrieve the NLCD 2019 image by filtering the NLCD ImageCollection and selecting the landcover band. Display the NLCD 2019 image on the interactive Map by using Map.addLayer. - -# Import the NLCD collection. -dataset = ee.ImageCollection('USGS/NLCD_RELEASES/2019_REL/NLCD') - -# Filter the collection to the 2019 product. -nlcd2019 = dataset.filter(ee.Filter.eq('system:index', '2019')).first() - -# Select the land cover band. -landcover = nlcd2019.select('landcover') - -# Display land cover on the map.Map.addLayer(landcover, {}, 'NLCD 2019') -Map - -Next, add the NLCD legend to the Map. The geemap package has several built-in legends, including the NLCD legend. Therefore, you can add the NLCD legend to the Map by using just one line of code (Map.add_legend). - -title = 'NLCD Land Cover Classification' -Map.add_legend(legend_title=title, builtin_legend='NLCD') - -Alternatively, if you want to add a custom legend, you can define a legend dictionary with labels as keys and colors as values, then you can use Map.add_legend to add the custom legend to the Map. - -legend_dict = { '11 Open Water': '466b9f', '12 Perennial Ice/Snow': 'd1def8', '21 Developed, Open Space': 'dec5c5', '22 Developed, Low Intensity': 'd99282', '23 Developed, Medium Intensity': 'eb0000', '24 Developed High Intensity': 'ab0000', '31 Barren Land (Rock/Sand/Clay)': 'b3ac9f', '41 Deciduous Forest': '68ab5f', '42 Evergreen Forest': '1c5f2c', '43 Mixed Forest': 'b5c58f', '51 Dwarf Scrub': 'af963c', '52 Shrub/Scrub': 'ccb879', '71 Grassland/Herbaceous': 'dfdfc2', '72 Sedge/Herbaceous': 'd1d182', '73 Lichens': 'a3cc51', '74 Moss': '82ba9e', '81 Pasture/Hay': 'dcd939', '82 Cultivated Crops': 'ab6c28', '90 Woody Wetlands': 'b8d9eb', '95 Emergent Herbaceous Wetlands': '6c9fb8'} -title = 'NLCD Land Cover Classification' -Map.add_legend(legend_title=title, legend_dict=legend_dict) - -The Map with the NLCD 2019 image and legend should look like Fig. F6.3.7. - -![Fig. F6.3.7 The NLCD 2019 image layer displayed in geemap](F6/image60.png) - - -The Map above shows only the NLCD 2019 image. To create an Earth Engine App for visualizing land cover change, we need a stack of NLCD images. Let’s print the list of system IDs of all available NLCD images. - -dataset.aggregate_array('system:id').getInfo() - -The output should look like this. - -['USGS/NLCD_RELEASES/2019_REL/NLCD/2001', -'USGS/NLCD_RELEASES/2019_REL/NLCD/2004', -'USGS/NLCD_RELEASES/2019_REL/NLCD/2006', -'USGS/NLCD_RELEASES/2019_REL/NLCD/2008', -'USGS/NLCD_RELEASES/2019_REL/NLCD/2011', -'USGS/NLCD_RELEASES/2019_REL/NLCD/2013', -'USGS/NLCD_RELEASES/2019_REL/NLCD/2016', -'USGS/NLCD_RELEASES/2019_REL/NLCD/2019'] - -Select the eight NLCD epochs after 2000. - -years = ['2001', '2004', '2006', '2008', '2011', '2013', '2016', '2019'] - -Define a function for filtering the NLCD ImageCollection by year and select the 'landcover' band. - -# Get an NLCD image by year. -def getNLCD(year): # Import the NLCD collection. dataset = ee.ImageCollection('USGS/NLCD_RELEASES/2019_REL/NLCD') # Filter the collection by year. nlcd = dataset.filter(ee.Filter.eq('system:index', year)).first() # Select the land cover band. landcover = nlcd.select('landcover'); return landcover - -Create an NLCD ImageCollection to be used in the split-panel map. - -# Create an NLCD image collection for the selected years. -collection = ee.ImageCollection(ee.List(years).map(lambda year: getNLCD(year))) - -Print out the list of system IDs of the selected NLCD images covering the contiguous United States. - -collection.aggregate_array('system:id').getInfo() - -The output should look like this. - -['USGS/NLCD_RELEASES/2019_REL/NLCD/2001', -'USGS/NLCD_RELEASES/2019_REL/NLCD/2004', -'USGS/NLCD_RELEASES/2019_REL/NLCD/2006', -'USGS/NLCD_RELEASES/2019_REL/NLCD/2008', -'USGS/NLCD_RELEASES/2019_REL/NLCD/2011', - - 'USGS/NLCD_RELEASES/2019_REL/NLCD/2013', -'USGS/NLCD_RELEASES/2019_REL/NLCD/2016', -'USGS/NLCD_RELEASES/2019_REL/NLCD/2019'] - -Next, create a list of labels to populate the dropdown list. - -labels = [f'NLCD {year}' for year in years] -labels - -The output should look like this. - -['NLCD 2001', -'NLCD 2004', -'NLCD 2006', -'NLCD 2008', -'NLCD 2011', -'NLCD 2013', - - 'NLCD 2016', -'NLCD 2019'] - -The last step is to create a split-panel map by passing the NLCD ImageCollection and list of labels to Map.ts_inspector. - -Map.ts_inspector(left_ts=collection, right_ts=collection, left_names=labels, right_names=labels) -Map - -The split-panel map should look like Fig. F6.3.8. - -![Fig. F6.3.8 A split-panel map for visualizing land cover change with geemap](F6/image34.png) - - -To visualize land cover change, choose one NLCD image from the left dropdown list and another image from the right dropdown list, then use the slider to swipe through to visualize land cover change interactively. Click the close button in the bottom-right corner to close the split-panel map and return to the NLCD 2019 image shown in Fig. F6.3.7. - -:::{.callout-note} -Code Checkpoint F63c. The book’s repository contains information about what your code should look like at this point. -::: -## Publishing an Earth Engine App Using a Local Web Server - -In this section, you will learn how to deploy an Earth Engine App using a local computer as a web server. Assume that you have completed Sect. 3 and created a Jupyter notebook named nlcd_app.ipynb. First, you need to download ngrok, a program that can turn your computer into a secure web server and connect it to the ngrok cloud service, which accepts traffic on a public address. Download ngrok from [https://ngrok.com](https://www.google.com/url?q=https://ngrok.com&sa=D&source=editors&ust=1671458841326366&usg=AOvVaw20VcUaADMd_4ptInlJt5X1) and unzip it to a directory on your computer, then copy nlcd_app.ipynb to the same directory. Open the Anaconda Prompt (on Windows) or the Terminal (on macOS/Linux) and enter the following commands. Make sure you change /path/to/ngrok/dir to your computer directory where the ngrok executable is located, e.g., ~/Downloads. - -cd /path/to/ngrok/dir - -conda activate gee -voila --no-browser nlcd_app.ipynb - -The output of the terminal should look like this. - -![Fig. F6.3.9 The output of the terminal running Voilà](F6/image38.png) - - -[Voilà](https://www.google.com/url?q=https://voila.readthedocs.io/en/stable/&sa=D&source=editors&ust=1671458841329454&usg=AOvVaw16GCXDSkxu4fUl6f6LFafy) can be used to run, convert, and serve a Jupyter notebook as a standalone app. Click the link (e.g., [http://localhost:8866](https://www.google.com/url?q=http://localhost:8866&sa=D&source=editors&ust=1671458841329837&usg=AOvVaw052EvlmJeJTl9zIHvFobLl)) shown in the terminal window to launch the interactive dashboard. Note that the port number is 8866, which is needed in the next step to launch ngrok. Open another terminal and enter the following command. - -cd /path/to/ngrok/dir - -ngrok http 8866 - -The output of the terminal should look like Fig. F6.3.10. Click the link shown in the terminal window to launch the interactive dashboard. The link should look like https://random-string.ngrok.io, which is publicly accessible. Anyone with the link will be able to launch the interactive dashboard and use the split-panel map to visualize NLCD land cover change. Keep in mind that the dashboard might take several seconds to load, so please be patient. - -![Fig. F6.3.10 The output of the terminal running ngrok](F6/image8.png) - - -To stop the web server, press Ctrl + C on both terminal windows. See below for some optional settings for running Voilà and ngrok. - -To show code cells from your App, run the following from the terminal. - -voila --no-browser --strip_sources=False nlcd_app.ipynb - -To protect your App with a password, run the following. - -ngrok http -auth="username:password" 8866 - -## Publish an Earth Engine App Using Cloud Platforms - -In this section, you will learn how to deploy an Earth Engine App on cloud platforms, such as Heroku and Google App Engine. Heroku is a “platform as a service” that enables developers to build, run, and operate applications entirely in the cloud. It has a free tier with limited computing hours, which would be sufficient for this lab. Follow the steps below to deploy the Earth Engine App on Heroku. - -First, go to [https://github.com/signup](https://www.google.com/url?q=https://github.com/signup&sa=D&source=editors&ust=1671458841334585&usg=AOvVaw3ZgJomY_82B7WiuJEZym5g) to sign up for a GitHub account if you don’t have one already. Once your GitHub account has been created, log into your account and navigate to the sample app repository: [https://github.com/giswqs/earthengine-apps](https://www.google.com/url?q=https://github.com/giswqs/earthengine-apps&sa=D&source=editors&ust=1671458841335025&usg=AOvVaw1IQJyv6EFLU5bCmsOsYkfX). Click the Fork button in the top-right corner to fork this repository into your account. Two important files in the repository are worth mentioning here: [requirements.txt](https://www.google.com/url?q=https://github.com/giswqs/earthengine-apps/blob/master/requirements.txt&sa=D&source=editors&ust=1671458841335397&usg=AOvVaw23DFUPWYa0pVU-kSzjGeHr) lists the required packages (e.g., geemap) to run the App, while [Procfile](https://www.google.com/url?q=https://github.com/giswqs/earthengine-apps/blob/master/Procfile&sa=D&source=editors&ust=1671458841335755&usg=AOvVaw2_QbFT8c3ahut0eHTdTmUh) specifies the commands that are executed by the App on startup. The content of [Procfile](https://www.google.com/url?q=https://github.com/giswqs/earthengine-apps/blob/master/Procfile&sa=D&source=editors&ust=1671458841336035&usg=AOvVaw2MdhcFIHFTRIbZBL-XZFCb) should look like this. - -web: voila --port=$PORT --no-browser --strip_sources=True --enable_nbextensions=True --MappingKernelManager.cull_interval=60 --MappingKernelManager.cull_idle_timeout=120 notebooks/ - -The above command instructs the server to hide the source code and periodically check for idle kernels—in this example, every 60 seconds—and cull them if they have been idle for more than 120 seconds. This can avoid idle kernels using up the server computing resources. The page served by Voilà will contain a list of any notebooks in the [notebooks](https://www.google.com/url?q=https://github.com/giswqs/earthengine-apps/tree/master/notebooks&sa=D&source=editors&ust=1671458841337601&usg=AOvVaw1pHGG9qzGbdG7h7tVr7C9e) directory. - -Next, go to [https://signup.heroku.com](https://www.google.com/url?q=https://signup.heroku.com&sa=D&source=editors&ust=1671458841338152&usg=AOvVaw03SbXy3optdTU9A2yG-bBq) to sign up for a Heroku account if you don’t have one already. Log into your account and click the New button in the top-right corner, then choose Create new app from the dropdown list (Fig. F6.3.11). - -![Fig. F6.3.11 Creating a new app in Heroku](F6/image5.png) - - -Choose a name for your App (Fig. F6.3.12). Note that the App name must be unique; if an App name has already been taken, you won’t be able to use it. - -![Fig. F6.3.12 Choosing an App name](F6/image73.png) - - -Once the Heroku App has been created, click the Deploy tab and choose GitHub as the deployment method. Connect to your GitHub account and enter earthengine-apps in the search box. The repository should be listed beneath the search box. Click the Connect button to connect the repository to Heroku (Fig. F6.3.13). - -![Fig. F6.3.13 Connecting a GitHub account to Heroku](F6/image27.png) - - -Under the same Deploy tab, scroll down and click Enable Automatic Deploys (Fig. F6.3.14). This will enable Heroku to deploy a new version of the App whenever the GitHub repository gets updated. - -![Fig. F6.3.14 Enabling Automatic Deploys on Heroku](F6/image31.png) - - -Since using Earth Engine requires authentication, we need to set the Earth Engine token as an environment variable so that the web App can pass the Earth Engine authentication. If you have completed Sect. 3, you have successfully authenticated Earth Engine on your computer and the token can be found in the following file path, depending on the operating system you are using. Note that you might need to show the hidden directories on your computer in order to see the .config folder under the home directory. - -Windows: C:UsersUSERNAME.configearthenginecredentialsLinux: /home/USERNAME/.config/earthengine/credentials -MacOS: /Users/USERNAME/.config/earthengine/credentials - -Open the credentials file using a text editor. Select and copy the long string wrapped by the double quotes (Fig. F6.3.15). Do not include the double quotes in the copied string. - -![Fig. F6.3.15 The Earth Engine authentication token](F6/image68.png) - - -Next, navigate to your web App on Heroku. Click the Settings tab, then click Reveal Config Vars on the page. Enter EARTHENGINE_TOKEN as the key and paste the string copied above as the value. Click the Add button to set EARTHENGINE_TOKEN as an environment variable (Fig. F6.3.16) that will be used by geemap to authenticate Earth Engine. - -![Fig. F6.3.16 Setting the Earth Engine token as an environment variable on Heroku](F6/image4.png) - - -The last step is to commit some changes to your forked GitHub repository to trigger Heroku to build and deploy the App. If you are familiar with [Git](https://www.google.com/url?q=https://git-scm.com/book/en/v2/Getting-Started-Installing-Git&sa=D&source=editors&ust=1671458841348079&usg=AOvVaw0sqY5Rw_8XGKlIcdI82G_f) commands, you can push changes to the repository from your local computer to GitHub. If you have not used Git before, you can navigate to your repository on GitHub and make changes directly using the browser. For example, you can navigate to [README.md](https://www.google.com/url?q=https://github.com/giswqs/earthengine-apps/blob/master/README.md&sa=D&source=editors&ust=1671458841348610&usg=AOvVaw1JFdaBCpAIcitrc-P7DHhR) and click the Edit icon on the page to start editing the file. Simply place the cursor at the end of the file and press Enter to add an empty line to the file, then click the Commit changes button at the bottom of the page to save changes. This should trigger Heroku to build the App. Check the build status under the latest activities of the Overview tab. Once the App has been built and deployed successfully, you can click the Open app button in the top-right corner to launch the web App (Fig. F6.3.17). When the App is open in your browser, click nlcd_app.ipynb to launch the split-panel map for visualizing land cover change. This is the same notebook that we developed in Sect. 3. The App URL should look like this: [https://APP-NAME.herokuapp.com/voila/render/nlcd_app.ipynb](https://www.google.com/url?q=https://app-name.herokuapp.com/voila/render/nlcd_app.ipynb&sa=D&source=editors&ust=1671458841349532&usg=AOvVaw2GAtxe9X-ymmUvBhjV_PBD). - -Congratulations! You have successfully deployed the Earth Engine App on Heroku. - -![Fig. F6.3.17 Setting the Earth Engine token as an environment variable on Heroku](F6/image51.png) - - -:::{.callout-note} -Code Checkpoint F63d. The book’s repository contains information about what your code should look like at this point. -::: -Question 1. What are the pros and cons of designing Earth Engine Apps using geemap and ipywidgets, compared to the JavaScript Earth Engine User Interface API? - -Question 2. What are the pros and cons of deploying Earth Engine Apps on Heroku, compared to a local web server and the Earth Engine Code Editor? - -## Synthesis {.unnumbered} - -Assignment 1. Replace the NLCD datasets with other multitemporal land cover datasets (e.g., United States Department of Agriculture National Agricultural Statistics Service Cropland Data Layers, visible in the Earth Engine Data Catalog) and modify the web App for visualizing land cover change using the chosen land cover datasets. Deploy the web App using multiple platforms (i.e., JavaScript Code Editor, ngrok, and Heroku). More land cover datasets can be found at [https://developers.google.com/earth-engine/datasets/tags/landcover](https://www.google.com/url?q=https://developers.google.com/earth-engine/datasets/tags/landcover&sa=D&source=editors&ust=1671458841352828&usg=AOvVaw0rYGxg8JD8U1d92jnx8V1i). - -## Conclusion {.unnumbered} +### Conclusion {.unnumbered} In this chapter, you learned how to design Earth Engine Apps using both the Earth Engine User Interface API (JavaScript) and geemap (Python). You also learned how to deploy Earth Engine Apps on multiple platforms, such as the JavaScript Code Editor, a local web server, and Heroku. The skill of designing and deploying interactive Earth Engine Apps is essential for making your research and data products more accessible to the scientific community and the general public. Anyone with the link to your web App can analyze and visualize Earth Engine datasets without needing an Earth Engine account. -## References {.unnumbered} +### References {.unnumbered} * Earth Engine User Interface API: [https://developers.google.com/earth-engine/guides/ui](https://www.google.com/url?q=https://developers.google.com/earth-engine/guides/ui&sa=D&source=editors&ust=1671458841355402&usg=AOvVaw1ZgOA5XvSnn5Uo83ccgHPT) * Earth Engine Apps: [https://developers.google.com/earth-engine/guides/apps](https://www.google.com/url?q=https://developers.google.com/earth-engine/guides/apps&sa=D&source=editors&ust=1671458841355890&usg=AOvVaw1ZYwRsMtxZ2RPAn-dC-E0n) @@ -2212,9 +1736,9 @@ Chapter F6.4: Combining R and Earth Engine :::{.callout-tip} -# Chapter Information +## Chapter Information -## Author {.unlisted .unnumbered} +#### Author {.unlisted .unnumbered} @@ -2222,12 +1746,12 @@ Cesar Aybar, David Montero, Antony Barja, Fernando Herrera, Andrea Gonzales, and -## Overview {.unlisted .unnumbered} +#### Overview {.unlisted .unnumbered} The purpose of this chapter is to introduce rgee, a non-official API for Earth Engine. You will explore the main features available in rgee and how to set up an environment that integrates rgee with third-party R and Python packages. After this chapter, you will be able to combine R, Python, and JavaScript in the same workflow. -## Learning Outcomes {.unlisted .unnumbered} +#### Learning Outcomes {.unlisted .unnumbered} * Becoming familiar with rgee, the Earth Engine R API interface. @@ -2236,7 +1760,7 @@ The purpose of this chapter is to introduce rgee, a non-official API for Earth E * Integrating Python and R packages using reticulate. * Combining Earth Engine JavaScript and Python APIs with R. -## Assumes you know how to:{.unlisted .unnumbered} +#### Assumes you know how to:{.unlisted .unnumbered} * Install the Python environment (Chap. F6.3). @@ -2246,7 +1770,7 @@ The purpose of this chapter is to introduce rgee, a non-official API for Earth E * Create Python virtual environments. ::: -## Introduction {.unlisted .unnumbered} +### Introduction {.unlisted .unnumbered} R is a popular programming language established in statistical science with large support in reproducible research, geospatial analysis, data visualization, and much more. To get started with R, you will need to satisfy some extra software requirements. First, install an up-to-date R version (at least 4.0) in your work environment. The installation procedure will vary depending on your operating system (i.e., Windows, Mac, or Linux). Hands-On Programming with R (Garrett Grolemund 2014, Appendix A) explains step by step how to proceed. We strongly recommend that Windows users install Rtools to avoid compiling external static libraries. @@ -2255,7 +1779,7 @@ If you are new to R, a good starting point is the book Geocomputation with R (Lo The following R packages must be installed (find more information in the R manual) in order to go through the practicum section. -# Use install.packages to install R packages from the CRAN repository. +## Use install.packages to install R packages from the CRAN repository. install.packages('reticulate') # Connect Python with R. install.packages('rayshader') # 2D and 3D data visualizations in R. install.packages('remotes') # Install R packages from remote repositories. @@ -2277,7 +1801,7 @@ Earth Engine officially supports client libraries only for the JavaScript and Py Figure F6.4.1 illustrates how rgee bridges the Earth Engine platform with the R ecosystem. When an Earth Engine request is created in R, rgee transforms this piece of code into Python. Next, the Earth Engine Python API converts the Python code into JSON. Finally, the JSON file request is received by the server through a Web REST API. Users could get the response using the getInfo method by following the same path in reverse. -## Installing rgee +### Installing rgee To run, rgee needs a Python environment with two packages: NumPy and earthengine-api. Because instructions change frequently, installation is explained at the following checkpoint: @@ -2288,17 +1812,17 @@ After installing both the R and Python requirements, you can now initialize your The Google Drive and GCS APIs will allow you to transfer your Earth Engine completed task exports to a local environment automatically. In these practice sessions, we will use only the Earth Engine and Google Drive APIs. Users that are interested in GCS can look up and explore the GCS vignette. To initialize your Earth Engine account alongside Google Drive, use the following commands. -# Initialize just Earth Engine +## Initialize just Earth Engine ee_Initialize() -# Initialize Earth Engine and GDee_Initialize(drive = TRUE) +## Initialize Earth Engine and GDee_Initialize(drive = TRUE) If the Google account is verified and the permission is granted, you will be directed to an authentication token. Copy and paste this token into your R console. Consider that the GCS API requires setting up credentials manually; look up and explore the rgee vignette for more information. The verification step is only required once; after that, rgee saves the credentials in your system. :::{.callout-note} Code Checkpoint F64b. The book’s repository contains information about what your code should look like at this point. ::: -## Creating a 3D Population Density Map with rgee and rayshader +### Creating a 3D Population Density Map with rgee and rayshader First, import the rgee, rayshader, and raster packages. @@ -2373,7 +1897,7 @@ render_snapshot( :::{.callout-note} Code Checkpoint F64c. The book’s repository contains information about what your code should look like at this point. ::: -## Displaying Maps Interactively +### Displaying Maps Interactively Similar to the Code Editor, rgee supports the interactive visualization of spatial Earth Engine objects by Map$addLayer. First, let’s import the rgee and cptcity packages. The cptcity R package is a wrapper to the cpt-city color gradients web archive. @@ -2433,17 +1957,17 @@ Map$addLayer(building) # eeObject is a ee$FeatureCollection The rgee functionality also supports the display of ee$ImageCollection via Map$addLayers (note the extra “s” at the end). Map$addLayers will use the same visualization parameters for all the images (Fig. F6.4.6). If you need different visualization parameters per image, use a Map$addLayer within a for loop. -# Define a ImageCollection +## Define a ImageCollection etp <- ee$ImageCollection$Dataset$MODIS_NTSG_MOD16A2_105 %>% ee$ImageCollection$select("ET") %>% ee$ImageCollection$filterDate('2014-04-01', '2014-06-01') -# Set viz paramsviz <- list( +## Set viz paramsviz <- list( min = 0, max = 300, palette = cpt(pal = "grass_bcyr", rev = TRUE) ) -# Print map results interactively +## Print map results interactively Map$setCenter(0, 0, 1) etpmap <- Map$addLayers(etp, visParams = viz) etpmap @@ -2492,7 +2016,7 @@ m2 | m1 :::{.callout-note} Code Checkpoint F64d. The book’s repository contains information about what your code should look like at this point. ::: -## Integrating rgee with Other Python Packages +### Integrating rgee with Other Python Packages As noted in Sect. 1, rgee set up a Python environment with NumPy and earthengine-api in your system. However, there is no need to limit it to just two Python packages. In this section, you will learn how to use the Python package ndvi2gif to perform a Normalized Difference Vegetation Index (NDVI) multi-seasonal analysis in the Ocoña Valley without leaving R. @@ -2549,10 +2073,10 @@ myclass <- ngif$NdviSeasonality( sat = 'Sentinel', # 'Sentinel', 'Landsat', 'MODIS', 'sar' key = 'max' # 'max', 'median', 'perc_90' ) -# Estimate the median of the yearly composites from 2016 to 2020. +## Estimate the median of the yearly composites from 2016 to 2020. median <- myclass$get_year_composite()$median() -# Estimate the median of the winter season composites from 2016 to 2020. +## Estimate the median of the winter season composites from 2016 to 2020. wintermax <- myclass$get_year_composite()$select('winter')$max() We can display maps interactively using the Map$addLayer (Fig. F6.4.9), and use the leaflet layers control to switch between basemaps. @@ -2572,7 +2096,7 @@ To get more information about the ndvi2gif package, visit its [GitHub](https://w :::{.callout-note} Code Checkpoint F64e. The book’s repository contains information about what your code should look like at this point. ::: -## Converting JavaScript Modules to R +### Converting JavaScript Modules to R In recent years, the Earth Engine community has developed a lot of valuable third-party modules. Some incredible ones are geeSharp (Zuspan 2020), ee-palettes (Donchyts et al. 2020), spectral (Montero 2021), and LandsatLST (Ermida et al. 2020). While some of these modules have been implemented in Python and JavaScript (e.g., geeSharp and spectral), most are available only for JavaScript. This is a critical drawback, because it divides the Earth Engine community by programming languages. For example, if an R user wants to use tagee (Safanelli et al. 2020), the user will have to first translate the entire module to R. @@ -2676,15 +2200,11 @@ Question 1. When and why might users prefer to use R instead of Python to connec Question 2. What are the advantages and disadvantages of using rgee instead of the Earth Engine JavaScript Code Editor? -## Synthesis {.unnumbered} - -Assignment 1. Estimate the Gaussian curvature map from a digital elevation model using rgee and rgeeExtra. Hint: Use the module 'users/joselucassafanelli/TAGEE:TAGEE-functions'. - -## Conclusion {.unnumbered} +### Conclusion {.unnumbered} In this chapter, you learned how to use Earth Engine and R in the same workflow. Since rgee uses reticulate, rgee also permits integration with third-party Earth Engine Python packages. You also learned how to use Map$addLayer, which works similarly to the Earth Engine User Interface API (Code Editor). Finally, we also introduced rgeeExtra, a new R package that extends the Earth Engine API and supports JavaScript module execution. -## References {.unnumbered} +### References {.unnumbered} Aybar C, Wu Q, Bautista L, et al (2020) rgee: An R package for interacting with Google Earth Engine. J Open Source Softw 5:2272. https://doi.org/10.21105/joss.02272 diff --git a/_quarto.yml b/_quarto.yml index 44775f1..488944b 100644 --- a/_quarto.yml +++ b/_quarto.yml @@ -22,6 +22,7 @@ book: - lights.qmd - refineries.qmd - ships.qmd + - blast.qmd repo-url: https://github.com/oballinger/GEE_OSINT/ google-analytics: G-RK9ZLZQ6GL repo-actions: [edit] @@ -40,6 +41,9 @@ format: dark: darkly light: cosmo code-copy: true + code-overflow: wrap + highlight-style: monokai.theme + linkcolor: "#34a832" pdf: documentclass: scrreprt diff --git a/blast.qmd b/blast.qmd new file mode 100644 index 0000000..f89cbe1 --- /dev/null +++ b/blast.qmd @@ -0,0 +1,429 @@ +# Blast Damage Assessment {.unnumbered} + +On August 4th, 2020, a warehouse full of fertilizer in the port of Beirut exploded. The blast killed over 200 people, injured thousands, and caused widespread damage to the city: + + +
+ +
+ + +Assessing blast damage is a common task for open source investigators, and satellite imagery analysis is one of the best tools we have at our disposal to analyze this sort of phenomenon. [NASA](https://earthobservatory.nasa.gov/images/147098/scientists-map-beirut-blast-damage) used Sentinel-1 imagery in the aftermath of the explosion to generate an estimated damage map. They explain that Sentinel-1 Synthetic Aperture Radar (SAR) imagery is good for this sort of task: + +"SAR instruments send pulses of microwaves toward Earth’s surface and listen for the reflections of those waves. The radar waves can penetrate cloud cover, vegetation, and the dark of night to detect changes that might not show up in visible light imagery. When Earth’s crust moves due to an earthquake, when dry land is suddenly covered by flood water, or when buildings have been damaged or toppled, the amplitude and phase of radar wave reflections changes in those areas and indicates to the satellite that something on the ground has changed." + +The NASA team produced this estimate the day after the explosion, which is very impressive. However, due to the quick turnaround, were pretty light on the description of their methodology (they didnt provide any code, or even explicity say how exactly they generated the estimate), making their analysis hard to replicate. They also failed to validate their results, which is a critical step in any analysis. + +In this case study we'll be developing our own change detection algorithm from scratch, applying to Sentinel-1 imagery of Beirut before and after the blast, and validating our results using building footprints and U.N. damage estimates as the ground truth. Below is the final result of the analysis, which shows building footprints colored according to the predicted level of damage they sustained from the blast: + +:::{.column-page} + + + +::: + + +## Change Detection + +There are many ways to detect change between two images. The simplest way would be to take an image taken before an event, and subtract it from an image taken after. This is a good way to get a general sense of where change has occurred, but if you only use two images (one before an event and another after), it would be difficult to differentiate between areas that have changed as a result of the event in question, and areas that have changed for other reasons. Things in Beirut (and cities in general) are constantly changing: construction, cars/planes/ships moving, vegetation growing, etc. So we wouldn't know if the change we're seeing is a result of the explosion or whether that area is generally prone to change. We can overcome this by comparing a bunch of pre-event images to a bunch of post-event images. This way we can see if the change we're seeing is consistent across all of the images. If it is, then we can be fairly confident that the change is a result of the event in question. The mean is simply the sum of all the values ($x_i$) in a set divided by the number of values ($n$): + +$$\large \overline{x} = \frac{1}{n} \sum_{i=1}^n x_i$$ + +But if we just take the average pixel value before and subtract the average pixel value after, we're not accounting for the *variability* of that pixel's values. For example, if we have a pixel that has had an average value of 1 for the month before the event, and a value of 2 in the month after the event, the difference is 1. If that pixel's value is extremely consistent (it never varies by more than 0.1), such a change would be very significant. But if that pixel's value is very variable (it varies by 2 or even 3 on a regular basis), then the change is not significant. So we need to account for the variability of the pixel's values using the standard deviation. It is calculated as the square root of the variance, which is the average of the squared differences from the mean: + +$$\large s = \sqrt{\frac{1}{n-1} \sum_{i=1}^n (x_i - \overline{x})^2}$$ + +With just the mean and the standard deviations of two sets of numbers, we can use a statistical test to determine whether the change in means is significant. The simplest way to do this is to use a pixelwise t-test, which is basically just a signal-to-noise ratio: it calculates the difference between two sample means (signal), and divides it by the standard deviations of both samples (noise). In this case, the two samples are the pre- and post-event images. The t-test is applied to each pixel in the image, allowing us to determine whether the change is statistically significant. Given two groups, $x_1$ before the event, and $x_2$ after the event, the $t$ statistic is calculated as: + +$$ \Large t = {\frac{\overline{x_1}-\overline{x_2}} {\sqrt{\frac{s^2_1}{n_1} + \frac{s^2_2}{n_2}}}} $$ + +Where: + +* $\overline{x}$: Sample Mean +* $s^2$: Sample Standard Deviation +* $n$: Number of observations + +This procedure gives us a number called a t-value, which is a measure of how many standard deviations the difference between the two means is. We're not going to get into the details here, but a rule of thumb is that if the t-value is greater than 2, then the difference between the two means is significant. If the t-value is less than 2, then the difference is not significant. We're going to calculate the t-value for each pixel in the image to determine whether that pixel has changed significantly following the event in question. You don't need to know the details of the t-test to understand the results (but hopefully you've got an intuition for what it's doing). If you're interested in learning more about statistical tests of this sort, I teach a course on Data Science at the University College London, and have made all of the lectures and courseware open-source. The T-test lecture is here. + +## Implementing a t-test in Earth Engine + +Now lets go about implementing this in Earth Engine. We'll start by centering the map on the port of Beirut, and setting the map to satellite view, and defining an area of interest (AOI) as a 3km buffer around the port: + +```js +Map.setCenter(35.51898, 33.90153, 15); + +Map.setOptions("satellite"); + +var aoi = ee.Geometry.Point(35.51898, 33.90153).buffer(3000); +``` +Next, let's define a function in earth engine that will perform the T-Test. The block of code below defines a function to implement a t-test for every pixel in a set of images. The function will be called 'ttest', and takes four arguments: + +* s1: the image collection +* shock: the date of the event +* pre_interval: the number of months before the event +* post_interval: the number of months after the event + +The function will return a t-value image, which we can use to determine whether a pixel has changed significantly following the event in question. + +```js +function ttest(s1, shock, pre_interval, post_interval) { + + // Convert the shock date to a date object + var shock = ee.Date(shock); + // Filter the image collection to the pre-event period + var pre = s1.filterDate( + shock.advance(ee.Number(pre_interval).multiply(-1), "month"), + shock + ); + // Filter the image collection to the post-event period + var post = s1.filterDate(shock, shock.advance(post_interval, "month")); + + // Calculate the mean, standard deviation, and number of images for the pre-event period + var pre_mean = pre.mean(); + var pre_sd = pre.reduce(ee.Reducer.stdDev()); + var pre_n = ee.Number(pre.filterBounds(aoi).size()); + + // Calculate the mean, standard deviation, and number of images for the pre-event period + var post_mean = post.mean(); + var post_sd = post.reduce(ee.Reducer.stdDev()); + var post_n = ee.Number(post.filterBounds(aoi).size()); + + // Calculate the pooled standard deviation + var pooled_sd = pre_sd + .multiply(pre_sd) + .multiply(pre_n.subtract(1)) + .add(post_sd.multiply(post_sd).multiply(post_n.subtract(1))) + .divide(pre_n.add(post_n).subtract(2)) + .sqrt(); + + // Calculate the denominator of the t-test + var denom = pooled_sd.multiply( + ee.Number(1).divide(pre_n).add(ee.Number(1).divide(post_n)).sqrt() + ); + + // Calculate the Degrees of Freedom, which is the number of observations minus 2 + var df = pre_n.add(post_n).subtract(2); + + print("Number of Images: ", df); + + // Calculate the t-test using the: + // mean of the pre-event period, + // the mean of the post-event period, + // and the pooled standard deviation + var change = post_mean + .abs() + .subtract(pre_mean.abs()) + .divide(denom) + .abs() + .subtract(2); + + // return the t-values for each pixel + return change; +} +``` +An important detail in the code above is that we've actually tweaked the t-test slightly, in two ways. + +First, the algorithm above returns tha *absolute* value of t (i.e. the absolute value of the difference between the two means). This is because we're interested in whether the pixel has changed at all, not whether it's changed in a particular direction. Second, we've subtracted 2 from the t-value. + +The t-value is a measure of how many standard deviations the difference between the two means is. Generally speaking, if the t-value is greater than 2, then the difference between the two means is considered statistically significant. 2 is a fairly abitrary cutoff, but it's the most commonly used one since it corresponds to the 95% confidence interval (i.e., theres less than a 5% chance of observing a difference that big due to random chance). Now we've got a function that can take an image collection and return a t-value image, where a value greater than 0 corresponds to a statistically significant change between the pre-event and post-event periods. + +## Filtering the Sentinel-1 Imagery + +We can't just blindly apply this algorithm to the entire image collection, because the image collection contains images from both ascending and descending orbits. We need to filter the image collection to the ascending and descending orbits, and then calculate the t-value for each orbit separately: this is because the satellite is viewing the scene from a completely different angle when it's in ascending and descending orbits, which will generate a lot of noise in our data. In fact, even when the satellite is either ascending or descending, we can have multiple images of the same place taken from slightly different orbital tracks because these overlap (see [this visualization of orbits](/ch1#orbits)). We need to filter the image collection to the relative orbit number that is most common within the image collection. For that, we define a new function called 'filter_s1', which takes a single argument: the path (either 'ASCENDING' or 'DESCENDING'). + +```js +function filter_s1(path) { + + // Filter the image collection to the ascending or descending orbit + var s1 = ee + .ImageCollection("COPERNICUS/S1_GRD") + .filter(ee.Filter.listContains("transmitterReceiverPolarisation", "VH")) + .filter(ee.Filter.eq("instrumentMode", "IW")) + .filter(ee.Filter.eq("orbitProperties_pass", path)) + .filterBounds(aoi) + .select("VH"); + + // Find the most common relative orbit number + var orbit = s1 + .aggregate_array("relativeOrbitNumber_start") + .reduce(ee.Reducer.mode()); + + // Filter the image collection to the most common relative orbit number + var s1 = s1.filter(ee.Filter.eq("relativeOrbitNumber_start", orbit)); + + // Calculate the t-test for the filtered image collection using the function we defined earlier + var change = ttest(s1, "2020-08-04", 12, 2); + + // Return the t-values for each pixel + return change; +} +``` +You'll notice that we've called the ttest function we defined earlier with four arguments: + +* s1: the Sentinel-1 image collection filtered to the ascending or descending orbit +* 2020-08-04: the date of the explosion +* 12: the number of months before the explosion to use for the pre-event period; I'm choosing to include the full year prior to the explosion to get a good baseline +* 2: the number of months after the explosion to use for the post-event period; I'm including 2 months after the explosion. Much less than that, and we risk missing the effects of the explosion. Much more than that, and we risk including the effects other changes that happened after the explosion, including the reconstruction effort. + +Now we want to apply this function to the image collection twice (once for each orbit) and then combine the two images into a single image. After that, we can clip it to the area of interest and display it on the map: + +```js + +// Call the filter_s1 function twice, once for each orbit, and then combine the two images into a single image +var composite = ee + .ImageCollection([filter_s1("ASCENDING"), filter_s1("DESCENDING")]) + .mean() + .clip(aoi); + +// Define a color palette +var palette = ["440154", "3b528b", "21918c", "5ec962", "fde725"]; + +// Add the composite to the map +Map.addLayer( + composite, + { min: 0, max: 4, opacity: 0.8, palette: palette }, + "change" +); +``` +The visualization parameters correspond the statitical significance of the change in pixel values. Using the Viridis color palette which ranges from purple to yellow, dark purple pixels indicate no significant change, and yellow pixels indicate a significant change with with 95% confidence. The brighter the yellow of a pixel, the more significant the change. + +![Pixelwise T-Test, 2020](./images/beirut/beirut_change_2020.jpg) + +This seems to be working quite well; but remember, ports are generally prone to change. The t-test is accounting for this by calculating each pixel's variance over the entire time period, but it's still possible that the change we're seeing is due to the port rather than the explosion. To test this, we can run the same algorithm on the same area, using the same date cutoff (August 4th), but in a different year; i've chosen 2018. This is what's known as a placebo test: if it's still showing loads of statistically significant change around the cutoff, our algorithm is probably picking up on port activity rather than the explosion. + +![Pixelwise T-Test, 2018](./images/beirut/beirut_change_2018.jpg) + +Compared to the 2020 image, there's a lot less yellow (significant change). That being said there are a few yellow areas. This could be due to a number of reasons: ships coming and going, cranes moving, and containers being loaded and unloaded would all register in the change detection algorithm. There are also a number of yellow specks throughout the city, which is also to be expected since cities are also generally in a state of flux. Construction, demolition, and even the growth of vegetation can all be detected by the algorithm. + +However, the scale and quantity of the change is nowhere near what it was for the 2020 image. This is a good sign that the algorithm detecting change resulting from the explosion. + +## Validation + +Great. We've developed our very own change detection algorithm in earth engine, applied it to the Beirut explosion, and it seems to be working using a basic placebo test. But how do we know that it's correctly predicting the *extent* of the damage, and not wildly over/underestimating? + +Given that this was a few years ago, we have the benefit of hindsight. In particular, the United Nations and the Municipality of Beirut have [published a report](https://unhabitat.org/sites/default/files/2020/10/municipality_of_beirut_-_beirut_explosion_rapid_assessment_report.pdf) on the damage caused by the explosion. This report includes estimates of the number of buildings damaged or destroyed by the explosion, as well as the number of people displaced. The report states that approximately 10,000 buildings were damaged within a 3km radius of the port. If our algorithm suggests that only 1,000 buildings were damaged, it's undershooting. If it suggests that 100,000 buildings were damaged, it's overshooting. + +Using building footprint data and the t-test image we just generated, we can generate an estimate of the number of damaged buildings according to our model. First, we want to generate a thresholded image, where pixels with a value greater than 0 are set to 1, and all other pixels are set to 0. We can then use this mask to reduce the building footprints to a single value for each building, where the value is the mean of the t-test image within the footprint. If the mean value is greater than 0, the building is damaged. If it's less than 0, the building is not damaged. + +```js +// Create a mask of the t-test image, where pixels with a value greater than 0 are set to 1, and all other pixels are set to 0 +var threshold = composite.updateMask(composite.gt(0)); + +// Load the building footprints +var buildings = ee + .FeatureCollection("projects/sat-io/open-datasets/MSBuildings/Lebanon") + .filterBounds(aoi); + +// Calculate the mean value of the t-test image within each building footprint +var damaged_buildings = threshold.reduceRegions({ + collection: buildings, + reducer: ee.Reducer.mean(), + scale: 1, +}); + +// Print the number of buildings with a mean value greater than 0 +// i.e., those displaying statistically significant change +print(damaged_buildings.filter(ee.Filter.gt("mean", 0)).size()); +``` + +The result is 9,256, which is pretty damn close to 10,000. We can also visualize the building footprints on the map, colored according the mean value of the t-test image within the footprint, where: + +* Blue = no damage +* Green = low damage +* Yellow/Orange = medium damage +* Red = high levels of damage + + +```js + +// Create an empty image +var empty = ee.Image().byte(); + +// Paint the building footprints onto the empty image +var outline = empty.paint({ + featureCollection: damaged_buildings, + color: "mean", + width: 5, +}); + +// Define a color palette +var building_palette = [ + "0034f5", + "1e7d83", + "4da910", + "b3c120", + "fcc228", + "ff8410", + "fd3000", +]; + +// Add the image to the map +Map.addLayer( + outline, + { palette: building_palette, min: 0, max: 2 }, + "Damaged Buildings" +); +``` + +The result naturally resembles the underlying t-test image, with high levels of damage concetrated around the port, and progressively decreasing damage with distance: + +![Building Footprints colored according to estimated blast damage](./images/beirut/beirut_footprints.jpg) + +To get a better sense of how these predicitons correspond to actual damage, we can zoom in and turn on the Google satellite basemap, which has imagery taken just after the explosion; you can still see capsized boats in the port. Zooming in to the epicentre, we can see several warehouses that were effectively vaporized. Our change detection algorithm picks up on a high degree of change, as indicated by the red outlines of the building footprints: + +![Predicted damage and optical satellite imagery in the Port of Beirut, August 2020](./images/beirut/beirut_footprints_port.jpg) + +This is pretty low-hanging fruit. Let's look at a different area, around 1.3km east from the epicentre with a mix of warehouses and residential buildings: + +![Area east of the port: 35.533194, 33.9024](./images/beirut/beirut_footprints_zoomed.jpg) + +Here, there's greater variation in the predictions. I've highlighted three areas. + +In Area A, we see a warehouse with a highily deformed roof; panels of corrugated iron are missing, and much of the roof is warped. The building footprint for this warehouse is red, suggesting that our algorithm correctly predicts a significant amount of blast damage. + +In Area B, we see a medium-rise building. If you look closely at the southern edge of the building, you'll see the siding has been completely torn off and is laying on the sidewalk. The bulding footprint is orange, suggesting a medium amount of change. We may be underestimating a bit here. + +In Area C, there are a bunch of high rise buildings clustered together. The building footprints are all blue, suggesting little to no damage. This is a bit of a surprise given how damaged areas A and B are. If you squint at the satellite image, it is indeed hard to tell if these buildings are damaged because we're looking at them from the top down, when much of the damage (e.g., the windows being blown out) would only be visible from the side. Indeed, our own estimate of the number of damaged buildings based on the algorithm we developed is about 8% shy of the U.N.'s estimate. This may be why. + +## Conclusion + +In this practical, we created a custom change detection algorithm that conducts a pixelwise t-test to detect change resulting from the 2020 explosion in the port of Beirut. By defining our own functions to do most of this analysis, we can apply the same workflow quite easily to a different context by simply moving the AOI and inputting the date of the shock. A placebo test showed that it's not just detecting general change in the area, but specifically change resulting from the explosion: when we keep everythgin the same but change the year of the shock, we see very little significant change being detected. Finally, by joining the predicted damage map to building footprints, we come up with an estimate of 9,256 damaged buildings, which is pretty close to the U.N.'s estimate of 10,000. That concludes the portion of this case study that deals with Earth Engine, but if you're interested in learning more about why we're coming up a bit short on the damage estimate (and some different ways of looking at the problem), read on. + +## Extension: Satellite Imagery and its Limits + +Though satellite imagery analysis is undoubtedly one of the best tools we have at our disposal to analyze this sort of phenomenon, it appears to systematically underestimate the extent of damage in Beirut. I outline an alternative approach using Open Street Map data to create a 3D model of Beirut and the explosion to analyze directional blast damage. Again, we're now leaving Earth Engine and moving to Blender, so if you're not interested in that feel free to skip ahead to the next case study. + +Below is one of the most viewed videos of the explosion: + +

Stunning video shows explosions just minutes ago at Beirut port pic.twitter.com/ZjltF0VcTr

— Borzou Daragahi 🖊🗒 (@borzou) August 4, 2020
+ +Geolocating this video was pretty simple thanks to the Greek Orthodox church (highlighted in green below) and the road leading to it (highlighted in blue). The red box indicates the likely location (33.889061, 35.515909) from which the person was filming: + +![](./images/beirut/IMG_2.png) + +The video shows heavy damage being sustained by areas well outside the zones classified as damaged in the maps above (both my own and NASA's). Indeed, substantial damage was reported several kilometers away. + +Why are satellite images underestimating damage in Beirut? Satellite images are taken from above, and are two-dimensional. Much of the damage caused by the blast, however, was directional; the pressure wave hit the sides of buildings, as shown in this diagram from a FEMA manual: + +![](./images/beirut/IMG_3.png) + +Areas close to the explosion suffered so much damage that it could be seen from above, but even if an apartment building had all of its windows blown out, this would not necessarily be visible in a top-down view. Even for radar, which does technically collect data in three dimensions, the angle problem remains; a high resolution radar might be able to tell you how tall an apartment complex is, but it won't give you a clear image of all sides. Case in point: the NASA damage map was created using Sentinel-1 SAR data. In a nutshell, damage assessment in this case is a three-dimensional problem, and remote sensing is a two-dimensional solution. + + + +### Creating a 3D model of Beirut + +To create a more accurate rendering of directional blast damage, three dimensional data are required. Data from Open Street Maps (OSM) contains information on both the "footprints" (i.e., the location and shape) as well as the height of buildings, which is enough to create a three dimensional model of Beirut. 3D rendering was done in Blender using the Blender-OSM add-on to import a satellite basemap, terrain raster, and OSM data. + +Geolocated videos of the blast can be used to verify and adjust the model. Below is a side-by-side comparison of the twitter video and a 3D rendition of OSM data: + +![](./images/beirut/IMG_4.png) + +Some slight adjustments to the raw OSM data were made to achieve the image on the right. The building footprints are generally very accurate and comprehensive in coverage, but the building height data does occasionally have to be adjusted manually. A simple and reliable way of doing this is to look at the shadows cast by the building on the satellite base map and scale accordingly. I also added a rough texture to the buildings to help differentiate them, and added the domed roof of the Greek Orthodox church for reference. + +For good measure, a second video is geolocated following the same procedure: + +

Another view of the explosions in Beirut pic.twitter.com/efT5VlpMkj

— Borzou Daragahi 🖊🗒 (@borzou) August 4, 2020
+ +The second pier (highlighted in green) and the angle (in blue) serve as references: + +![](./images/beirut/IMG_5.png) + +The video was taken from the rooftop of a japanese restaurant called Clap Beirut (in red above). This is confirmed by a picture of the rooftop bar on google images, which matches the bar that can be seen at 0:02 in the twitter video. Below is a comparison of the video view and the 3D OSM model: + +![](./images/beirut/IMG_6.png) + +Though somewhat grainy, the basemap on the OSM rendering shows the same parking lot in the foreground, the second pier, and the same two buildings highlighted in yellow. Having created a 3D model of Beirut using OSM data, we can now simulate how the explosion would interact with the cityscape. + + + +### Using a Viewshed Analysis to Assess Blast Exposure + +As the pressure wave moved through the Beirut, some buildings bore the full force of the explosion, while others were partially shielded by taller structures. A viewshed analysis can be conducted to identify surfaces that were directly exposed to the explosion by creating a lighting object at ground zero; areas that are lit up experienced unobstructed exposure to the blast: + +![](./images/beirut/GIF_1.gif) + +Pressure waves, like sound, are capable of diffraction (beding around small obstructions). To roughly simluate this, the lighting object is gradually raised, allowing the light to pass "around" obstructions. Warehouses on the Eastern side of the docks, as well as the first row of apartment buildings facing the docks are immediately affected. As the lighting object rises above the warehouse, more areas suffer direct exposure. + +Using two lighting objects-- a red one at 10 meters and a blue one at 20 meters above the warehouse at ground zero-- the intensity of the blast in different areas is highlighted; red areas suffered direct exposure, blue areas suffered partially obstructed exposure, and black areas were indirectly exposed. + +![](./images/beirut/IMG_7.png) + +In the immediate vicinity of the explosion the large "L" shaped building (Lebanon's strategic grain reserve) is bright red, and was barely left standing. It absorbed a large amount of the blast, shielding areas behind it and thereby casting a long blue shadow to the West. If one refers back to the satellite damage maps above, there appears to be significantly less damage in the area just West of ("behind") the grain silo, roughly corresponding to the blue shadow above. While these areas were still heavily damaged, they seem to have suffered less damage than areas of equal distance to the East. + + + +### Accounting for Diffraction + +The viewshed analysis tells us which sides of a building are exposed to the blast, but it's a pretty rough approximation of the way the pressure wave would respond to obstacles in its path. As previously mentioned, pressure waves behave much like sound waves or waves in water: they bounce off of objects, move around obstructions, and gradually fade. + +To get a more precise idea of the way in which the blast interacted with the urban environment, we can model the blast as an actual wave using the "dynamic wave" feature in Blender. This effectively involves creating a two-dimensional plane, telling it to behave like water, and simulating an object being dropped into the water. By putting an obstruction in this plane, we can see how the wave responds to it. As an example, the grain silo has been isolated below: + +![](./images/beirut/GIF_2.gif) + +As the blast hits the side of the silo, it is reflected. Two large waves can be seen traveling to the right: the initial blast wave, and the reflection from the silo which rivals the initial wave in magnitude. To the left, the wave travels around the silo but is significantly weakened. + +Broadening the focus and adding the rest of the OSM data back in, we can observe how the pressure wave interacted with buildings on the waterfront: + +![](./images/beirut/GIF_3.gif) + +The warehouses on the docks were omitted to emphasize the interaction between the pressure wave and the waterfront buildings; their light metal structure and low height means they would have caused little reflection anyway. The general pattern of the dynamic wave is consistent with the viewshed, but adds a layer of detail. The blast is reflected off of the silo towards the East, leading to a double hit. Though the wave still moves around the silo to the West, the pressure is diminished. Once the wave hits the highrises, the pattern becomes noisy as the wave both presses forward into the mainland and is reflected back towards the pier. + + + +### Modeling the Pressure Wave + +Now that we've accounted for the directionality of the blast and the influence of buildings, we can model the pressure wave itself. An expanding sphere centered at ground zero is used to model the progression of the pressure wave through the city. To get a visual sense of the blast's force, the color of the sphere will be a function of the pressure exerted by pressure wave. + +The pressure exerted by the explosion in kilopascals (kPa) at various distances can be calculated using the DoD's Blast Effects Computer, which allows users to input variables such as the TNT equivalent of the ordnance, storage method, and elevation. Though there are several estimates, the blast was likely equivalent to around 300 tons of TNT. The direct "incident pressure" of the pressure wave is shown in blue. However, pressure waves from explosions that occur on the ground are reflected upwards, amplifying the total pressure exerted by the blast. This "reflected pressure" is shown in orange: + + + +For reference, 137 kPa results in 99% fatalities, 68 kPa is enough to cause structural damage to most buildings, and 20 kPa results in serious injuries. 1-6 kPa is enough to break an average window. At 1km, the reflected pressure of the blast (18 kPa) was still enough to seriously injure. Precisely calculating the force exerted by an explosion is exceptionally complicated, however, so these numbers should be treated as rough estimates. Further analysis of the damage caused by blasts blast can be derived from the UN's Explosion Consequences Analysis calculator which provides distance values for different types of damage and injuries. + +Linking the values in this graph to the color of the pressure wave sphere provides a visual representation of the blast's force as it expands. An RGB color scale corresponds to the blast's overpressure at three threshold values. + +![](./images/beirut/beirut.gif) + + +By keeping the lighting object from the viewshed analysis and placing it within the expanding sphere of the pressure wave, we combine two key pieces of information: the pressure exerted by the blast (the color of the sphere), and the level of directional exposure (brightness). + +Now, referring back to the two geolocated twitter videos from earlier, we can recreate the blast in our 3D model and get some new insights. Below is a side-by-side comparison of the first video and the 3D model: + +![](./images/beirut/GIF_5.gif) + +Judging by the twitter video alone, it would be very hard to tell the fate of the person filming or the damage caused to the building that they were in. However, the 3D model shows that despite having an unobstructed view of the explosion, the incident pressure of the pressure wave had decreased significtantly by the time it reached the viewing point. The blue-green color corresponds to roughly 15 kPa-- enough to injure and break windows, but not enough to cause structural damage to the building. + +The second twitter video was taken slightly closer to ground zero, but the view was partially obstructed by the grain silo: + +![](./images/beirut/GIF_6.gif) + +Though the pressure wave probably exerted more pressure compared to the first angle, the partial obstruction of the grain silo likely tempered the force of the blast. + +### Assessing Damage to the Skyline Tower + +As a concrete example of how this approach can be used to assess damage (or predict it, if one had the foresight), let us consider the Skyline Tower, pictured below following the explosion: + +![](./images/beirut/IMG_8.png) + +This partial side view shows two faces of the building, labelled "A" and "B" above. Side A was nearly perpendicular to the blast, and just 600 m from ground zero. Based on the previous modeling, the pressure wave exerted roughly 40 kPa on this side of the building. The corner where sides A and B meet, highlighted in green, shows total destruction of windows, removal of most siding panels, and structural damage. The back corner, highlighted in red, shows many windows still intact, indicating that the maximum overpressure on this side of the building likely didn't exeed 10 kPa. In other words, standing on the front balcony would likely have led to serious injury but standing on the back balcony would have been relatively safe. + +The animation below shows the Skyline Tower as it is hit by the pressure wave, with sides A and B labeled: + +![](./images/beirut/GIF_7.gif) + +The bright green color of the pressure wave indicates a strong likelihood of structural damage. Side A can be seen taking a direct hit, while side B is angled slighly away. Despite not being directly exposed to the blast, it likely still took reflective damage from some of the neighbouring buildings. Both the incident overpressure indicated by the color of the sphere, as well as the relative brightness of sides A and B both correspond closely to the observed damage taken by the Skyline Tower. + +### Further Research + +Though satellite imagery analysis is an indispensable tool in disaster response, it has limitations. Urban blast damage in particular is difficult to assess accurately because it is highly directional and much of it cannot be seen from a bird's eye view. Using free and open source tools, an interactive 3D model of an urban explosion can be generated, allowing for a highly detailed investigation of directional blast damage. This can be achieved in three steps: + +First, creating a 3D model of the urban area using Blender and Open Street Maps data. Second, conducting a viewshed analysis using lighting objects to gauge levels of unobstructed exposure to the pressure wave. Third, modeling the explosion using geolocated videos of the event and ordnance calculators. For added detail, a dynamic wave analysis can be used to more precisely model how the pressure wave interacts with buildings. + +Once properly modeled, the explosion can be viewed from any angle in the city. The viewshed analysis can be calibrated more finely by ground-truthing various damage levels (e.g. broken windows) at different locations. In the absence of an official address registry in Beirut, OSM is already being used by the Lebanese Red Cross (donate here) to conduct neighborhood surveys assessing blast damage. As such, this type of damage analysis can quickly be integrated into relief efforts, adapted to model disasters in different cities, and can even be used to simulate the destructive potential of hypothetical explosions to promote readiness. + +Speacial thanks to my nuclear physicist brother, Sean, for making sure I didn't commit too many crimes against Physics. \ No newline at end of file diff --git a/docs/F1.html b/docs/F1.html index e771483..00bbcc1 100644 --- a/docs/F1.html +++ b/docs/F1.html @@ -290,69 +290,44 @@ gtag('config', 'G-RK9ZLZQ6GL', { 'anonymize_ip': true});

Table of contents

-
-

Synthesis

-

Question 1. Look at the MODIS example (Sect. 3.2), which uses the median reducer. Try modifying the reducer to be unweighted, either by specifying unweighted or using an identifier reducer like max. What happens, and why?

-

Question 2. Calculate zonal statistics for your own buffered points or polygons using a raster and reducer of interest. Be sure to consider the spatial scale of the raster and whether a weighted or unweighted reducer would be more appropriate for your interests.

-

If the point or polygon file is stored in a local shapefile or CSV file, first upload the data to your Earth Engine assets. All columns in your vector file, such as the plot name, will be retained through this process. Once you have an Earth Engine table asset ready, import the asset into your script by hovering over the name of the asset and clicking the arrow at the right side, or by calling it in your script with the following code.

-

var pts = ee.FeatureCollection(‘users/yourUsername/yourAsset’);

-

If you prefer to define points or polygons dynamically rather than loading an asset, you can add them to your script using the geometry tools. See Chap. F2.1 and F5.0 for more detail on adding and creating vector data.

-

Question 3. Try the code from Sect. 4.3 using the MODIS data and the first point from the pts variable. Among other modifications, you will need to create a buffer for the point, take a single MODIS image from the collection, and change visualization parameters.

-
    -
  • Think about the CRS in the code: The code reprojects to EPSG:5070, but MODIS is collected in the sinusoidal projection SR-ORG:6974. Try that CRS and describe how the image changes.
  • -
  • Is the count reducer weighted or unweighted? Give an example of a circumstance to use a weighted reducer and an example for an unweighted reducer. Specify the buffer size you would use and the spatial resolution of your dataset.
  • -
-

Question 4. In the examples above, only a single ee.Reducer is passed to the zonalStats function, which means that only a single statistic is calculated (for example, zonal mean or median or maximum). What if you want multiple statistics—can you alter the code in Sect. 3.1 to (1) make the point buffer 500 instead of 45; (2) add the reducer parameter to the params dictionary; and (3) as its argument, supply a combined ee.Reducer that will calculate minimum, maximum, standard deviation, and mean statistics?

-

To achieve this you’ll need to chain several ee.Reducer.combine functions together. Note that if you accept all the individual ee.Reducer and ee.Reducer.combine function defaults, you’ll run into two problems related to reducer weighting differences, and whether or not the image inputs are shared among the combined set of reducers. How can you manipulate the individual ee.Reducer and ee.Reducer.combine functions to achieve the goal of calculating multiple zonal statistics in one call to the zonalStats function?

-
-
-

Conclusion

+
+

Conclusion

In this chapter, you used functions containing optional parameters to extract raster values for collocated points. You also learned how to buffer points, and apply weighted and unweighted reducers to get different types of zonal statistics. These functions were applied to three examples that differed by raster dataset, reducer, spatial resolution, and scale. Lastly, you covered related topics like weighting of reducers and buffer visualization. Now you’re ready to apply these ideas to your own work!

-
-

References

+
+

References

Cansler CA, McKenzie D (2012) How robust are burn severity indices when applied in a new region? Evaluation of alternate field-based and remote-sensing methods. Remote Sens 4:456–483. https://doi.org/10.3390/rs4020456

Miller JD, Thode AE (2007) Quantifying burn severity in a heterogeneous landscape with a relative version of the delta Normalized Burn Ratio (dNBR). Remote Sens Environ 109:66–80. https://doi.org/10.1016/j.rse.2006.12.006

-
-

10 Advanced Vector Operations

-

:::{.callout-tip} # Chapter Information

-
-

Author

+
+

6.4 Advanced Vector Operations

+
+
+
+ +
+
+Chapter Information +
+
+
+
+

Author

Ujaval Gandhi

-
-

Overview

+
+

Overview

This chapter covers advanced techniques for visualizing and analyzing vector data in Earth Engine. There are many ways to visualize feature collections, and you will learn how to pick the appropriate method to create visualizations, such as a choropleth map. We will also cover geoprocessing techniques involving multiple vector layers, such as selecting features in one layer by their proximity to features in another layer and performing spatial joins.

-
-

Learning Outcomes

+
+

Learning Outcomes

  • Visualizing any vector dataset and creating a thematic map.
  • Understanding joins in Earth Engine.
  • Carrying out geoprocessing tasks with vector layers in Earth Engine.
-
-

Assumes you know how to:

+
+

Assumes you know how to:

  • Filter a FeatureCollection to obtain a subset (Chap. F5.0, Chap. F5.1).
  • Write a function and map it over a FeatureCollection (Chap. F5.1, Chap. F5.2).
-
-

10.1 Visualizing Feature Collections

+
+
+
+

6.4.1 Visualizing Feature Collections

There is a distinct difference between how rasters and vectors are visualized. While images are typically visualized based on pixel values, vector layers use feature properties (i.e., attributes) to create a visualization. Vector layers are rendered on the Map by assigning a value to the red, green, and blue channels for each pixel on the screen based on the geometry and attributes of the features. The functions used for vector data visualization in Earth Engine are listed below in increasing order of complexity.

  • Map.addLayer: As with raster layers, you can add a FeatureCollection to the Map by specifying visualization parameters. This method supports only one visualization parameter: color. All features are rendered with the specified color.
  • @@ -1909,17 +1864,19 @@ Note
  • style: This is the most versatile function. It can apply a different style to each feature, including color, pointSize, pointShape, width, fillColor, and lineType.

In the exercises below, we will learn how to use each of these functions and see how they can generate different types of maps.

-
-

10.1.1 Creating a Choropleth Map

+
+

6.4.1.1 Creating a Choropleth Map

We will use the TIGER: US Census Blocks layer, which stores census block boundaries and their characteristics within the United States, along with the San Francisco neighborhoods layer from Chap. F5.0 to create a population density map for the city of San Francisco.

We start by loading the census blocks and San Francisco neighborhoods layers. We use ee.Filter.bounds to filter the census blocks layer to the San Francisco boundary.

-

var blocks = ee.FeatureCollection(‘TIGER/2010/Blocks’);
-var roads = ee.FeatureCollection(‘TIGER/2016/Roads’);
-var sfNeighborhoods = ee.FeatureCollection( ‘projects/gee-book/assets/F5-0/SFneighborhoods’);

-

var geometry = sfNeighborhoods.geometry();
-Map.centerObject(geometry);

-
// Filter blocks to the San Francisco boundary.  
-var sfBlocks = blocks.filter(ee.Filter.bounds(geometry));
+
var blocks = ee.FeatureCollection('TIGER/2010/Blocks');  
+var roads = ee.FeatureCollection('TIGER/2016/Roads');  
+var sfNeighborhoods = ee.FeatureCollection(   'projects/gee-book/assets/F5-0/SFneighborhoods');  
+  
+var geometry = sfNeighborhoods.geometry();  
+Map.centerObject(geometry);  
+  
+// Filter blocks to the San Francisco boundary.  
+var sfBlocks = blocks.filter(ee.Filter.bounds(geometry));

The simplest way to visualize this layer is to use Map.addLayer (Fig. F5.3.1). We can specify a color value in the visParams parameter of the function. Each census block polygon will be rendered with stroke and fill of the specified color. The fill color is the same as the stroke color but has a 66% opacity.

// Visualize with a single color.  
 Map.addLayer(sfBlocks, {  
@@ -1932,15 +1889,20 @@ Map.centerObject(geometry);

The census blocks table has a property named ‘pop10’ containing the population totals as of the 2010 census. We can use this to create a choropleth map showing population density. We first need to compute the population density for each feature and add it as a property. To add a new property to each feature, we can map a function over the FeatureCollection and calculate the new property called ‘pop_density’. Earth Engine provides the area function, which can calculate the area of a feature in square meters. We convert it to square miles and calculate the population density per square mile.

// Add a pop_density column.  
-var sfBlocks = sfBlocks.map(function(f) {   // Get the polygon area in square miles.   var area_sqmi = f.area().divide(2.59e6);   var population = f.get('pop10');   // Calculate population density.   var density = ee.Number(population).divide(area_sqmi);   return f.set({       'area_sqmi': area_sqmi,       'pop_density': density  
-   });  
-});
+var sfBlocks = sfBlocks.map(function(f) { // Get the polygon area in square miles. +var area_sqmi = f.area().divide(2.59e6); +var population = f.get('pop10'); // Calculate population density. +var density = ee.Number(population).divide(area_sqmi); +return f.set({ + 'area_sqmi': area_sqmi, + 'pop_density': density + }); +});

Now we can use the paint function to create an image from this FeatureCollection using the pop_density property. The paint function needs an empty image that needs to be cast to the appropriate data type. Let’s use the aggregate_stats function to calculate basic statistics for the given column of a FeatureCollection.

// Calculate the statistics of the newly computed column.  
 var stats = sfBlocks.aggregate_stats('pop_density');  
 print(stats);

You will see that the population density values have a large range. We also have values that are greater than 100,000, so we need to make sure we select a data type that can store values of this size. We create an empty image and cast it to int32, which is able to hold large integer values.

-

D

The result is an image with pixel values representing the population density of the polygons. We can now use the standard image visualization method to add this layer to the Map (Fig. F5.3.2). Then, we need to determine minimum and maximum values for the visualization parameters.A reliable technique to produce a good visualization is to find minimum and maximum values that are within one standard deviation. From the statistics that we calculated earlier, we can estimate good minimum and maximum values to be 0 and 50000, respectively.

var palette = [‘fee5d9’, ‘fcae91’, ‘fb6a4a’, ‘de2d26’, ‘a50f15’];
var visParams = {
@@ -1956,8 +1918,8 @@ Map.addLayer(sfBlocksPaint.clip(geometry), visParams, ‘Population Density’);

-
-

10.1.2 Creating a Categorical Map

+
+

6.4.1.2 Creating a Categorical Map

Continuing the exploration of styling methods, we will now learn about draw and style. These are the preferred methods of styling for points and line layers. Let’s see how we can visualize the TIGER: US Census Roads layer to create a categorical map.

We start by filtering the roads layer to the San Francisco boundary and using Map.addLayer to visualize it.

// Filter roads to San Francisco boundary.  
@@ -2066,8 +2028,8 @@ Note
 

Save your script for your own future use, as outlined in Chap. F1.0. Then, refresh the Code Editor to begin with a new script for the next section.

-
-

10.2 Joins with Feature Collections

+
+

6.4.2 Joins with Feature Collections

Earth Engine was designed as a platform for processing raster data, and that is where it shines. Over the years, it has acquired advanced vector data processing capabilities, and users are now able to carry out complex geoprocessing tasks within Earth Engine. You can leverage the distributed processing power of Earth Engine to process large vector layers in parallel.

This section shows how you can do spatial queries and spatial joins using multiple large feature collections. This requires the use of joins. As described for Image Collections in Chap. F4.9, a join allows you to match every item in a collection with items in another collection based on certain conditions. While you can achieve similar results using map and filter, joins perform better and give you more flexibility. We need to define the following items to perform a join on two collections.

    @@ -2075,8 +2037,8 @@ Note
  1. Join type: While the filter determines which features will be joined, the join type determines how they will be joined. There are many join types, including simple join, inner join, and save-all join.

Joins are one of the harder skills to master, but doing so will help you perform many complex analysis tasks within Earth Engine. We will go through practical examples that will help you understand these concepts and the workflow better.

-
-

10.2.1 Selecting by Location

+
+

6.4.2.1 Selecting by Location

In this section, we will learn how to select features from one layer that are within a specified distance from features in another layer. We will continue to work with the San Francisco census blocks and roads datasets from the previous section. We will implement a join to select all blocks in San Francisco that are within 1 km of an interstate highway.

We start by loading the census blocks and roads collections and filtering the roads layer to the San Francisco boundary.

var blocks = ee.FeatureCollection(‘TIGER/2010/Blocks’);
@@ -2132,8 +2094,8 @@ Map.addLayer(closeBlocksDrawn, {}, ‘Blocks within 1km’);

-
-

10.2.2 Spatial Joins

+
+

6.4.2.2 Spatial Joins

A spatial join allows you to query two collections based on the spatial relationship. We will now implement a spatial join to count points in polygons. We will work with a dataset of tree locations in San Francisco and polygons of neighborhoods to produce a CSV file with the total number of trees in each neighborhood.

The San Francisco Open Data Portal maintains a street tree map dataset that has a list of street trees with their latitude and longitude. We will also use the San Francisco neighborhood dataset from the same portal. We downloaded, processed, and uploaded these layers as Earth Engine assets for use in this exercise. We start by loading both layers and using the paint and style functions, covered in Sect. 1, to visualize them (Fig. F5.3.7).

var sfNeighborhoods = ee.FeatureCollection( ‘projects/gee-book/assets/F5-0/SFneighborhoods’);
@@ -2227,12 +2189,8 @@ Note

-
-

Synthesis

-

Assignment 1. What join would you use if you wanted to know which neighborhood each tree belongs to? Modify the code above to do a join and post-process the result to add a neighborhood property to each tree point. Export the results as a shapefile.

-
-
-

Conclusion

+
+

Conclusion

This chapter covered visualization and analysis using vector data in Earth Engine. You should now understand different functions for FeatureCollection visualization and be able to create thematic maps with vector layers. You also learned techniques for doing spatial queries and spatial joins within Earth Engine. Earth Engine is capable of handling large feature collections and can be effectively used for many spatial analysis tasks.

diff --git a/docs/F6.html b/docs/F6.html index 3012ae9..c26980f 100644 --- a/docs/F6.html +++ b/docs/F6.html @@ -291,81 +291,54 @@ gtag('config', 'G-RK9ZLZQ6GL', { 'anonymize_ip': true});

Table of contents

@@ -391,8 +364,8 @@ gtag('config', 'G-RK9ZLZQ6GL', { 'anonymize_ip': true});

Although you now know the most basic fundamentals of Earth Engine, there is still much more that can be done. The Part presents some advanced topics that can help expand your skill set for doing larger and more complex projects. These include tools for sharing code among users, scaling up with efficient project design, creating apps for non-expert users, and combining R with other information processing platforms.

-
-

8 Advanced Raster Visualization

+
+

7.1 Advanced Raster Visualization

@@ -403,16 +376,16 @@ Chapter Information
-
-

Author

+
+

Author

Gennadii Donchyts, Fedor Baart

-
-

Overview

+
+

Overview

This chapter should help users of Earth Engine to better understand raster data by applying visualization algorithms such as hillshading, hill shadows, and custom colormaps. We will also learn how image collection datasets can be explored by animating them as well as by annotating with text labels, using, for example, attributes of images or values queried from images.

-
-

Learning Outcomes

+
+

Learning Outcomes

  • Understanding why perceptually uniform colormaps are better to present data and using them efficiently for raster visualization.
  • Using palettes with images before and after remapping values.
  • @@ -421,8 +394,8 @@ Chapter Information
  • Adding hillshading and shadows to help visualize raster datasets.
-
-

Assumes you know how to:

+
+

Assumes you know how to:

  • Import images and image collections, filter, and visualize (Part F1).
  • Write a function and map it over an ImageCollection (Chap. F4.0).
  • @@ -431,8 +404,8 @@ Chapter Information
-
-

Introduction

+
+

Introduction

Visualization is the step to transform data into a visual representation. You make a visualization as soon as you add your first layer to your map in Google Earth Engine. Sometimes you just want to have a first look at a dataset during the exploration phase. But as you move towards the dissemination phase, where you want to spread your results, it is good to think about a more structured approach to visualization. A typical workflow for creating visualization consists of the following steps:

  • Defining the story (what is the message?)
  • @@ -446,8 +419,8 @@ Chapter Information

    A good standard work on all the choices that one can make while creating a visualization is provided by the Grammar of Graphics (GoG) by Wilkinson (1999). It was the inspiration behind many modern visualization libraries (ggplot, vega). The main concept is that you can subdivide your visualization into several aspects.

    In this chapter, we will cover several aspects mentioned in the Grammar of Graphics to convert (raster) data into visual elements. The accurate representation of data is essential in science communication. However, color maps that visually distort data through uneven color gradients or are unreadable to those with color-vision deficiency remain prevalent in science (Crameri, 2020). You will also learn how to add annotation text and symbology, while improving your visualizations by mixing images with hillshading as you explore some of the amazing datasets that have been collected in recent years in Earth Engine.

-
-

8.1 Palettes

+
+

7.1.1 Palettes

In this section we will explore examples of colormaps to visualize raster data. Colormaps translate values to colors for display on a map. This requires a set of colors (referred to as a “palette” in Earth Engine) and a range of values to map (specified by the min and max values in the visualization parameters).

There are multiple types of colormaps, each used for a different purpose. These include the following:

Sequential: These are probably the most commonly used colormaps, and are useful for ordinal, interval, and ratio data. Also referred to as a linear colormap, a sequential colormap looks like the viridis colormap (Fig. F6.0.1) from matplotlib. It is popular because it is a perceptual uniform colormap, where an equal interval in values is mapped to an equal interval in the perceptual colorspace. If you have a ratio variable where zero means nothing, you can use a sequential colormap starting at white, transparent, or, when you have a black background, at black—for example, the turku colormap from Crameri (Fig. F6.0.1). You can use this for variables like population count or gross domestic product.

@@ -563,8 +536,8 @@ Note
-
-

8.2 Remapping and Palettes

+
+

7.1.2 Remapping and Palettes

Classified rasters in Earth Engine have metadata attached that can help with analysis and visualization. This includes lists of the names, values, and colors associated with class. These are used as the default color palette for drawing a classification, as seen next. The USGS National Land Cover Database (NLCD) is one such example. Let’s access the NLCD dataset, name it nlcd, and view it (Fig. F6.0.4) with its built-in palette.

// Advanced remapping using NLCD.  
 // Import NLCD.  
@@ -646,8 +619,8 @@ Note
 
-
-

8.3 Annotations

+
+

7.1.3 Annotations

Annotations are the way to visualize data on maps to provide additional information about raster values or any other data relevant to the context. In this case, this additional information is usually shown as geometries, text labels, diagrams, or other visual elements. Some annotations in Earth Engine can be added by making use of the ui portion of the Earth Engine API, resulting in graphical user interface elements such as labels or charts added on top of the map. However, it is frequently useful to render annotations as a part of images, such as by visualizing various image properties or to highlight specific areas.

In many cases, these annotations can be mixed with output images generated outside of Earth Engine, for example, by post-processing exported images using Python libraries or by annotating using GIS applications such as QGIS or ArcGIS. However, annotations could also be also very useful to highlight and/or label specific areas directly within the Code Editor. Earth Engine provides a sufficiently rich API to turn vector features and geometries into raster images which can serve as annotations. We recommend checking the ee.FeatureCollection.style function in the Earth Engine documentation to learn how geometries can be rendered.

For textual annotation, we will make use of an external package ‘users/gena/packages:text’ that provides a way to render strings into raster images directly using the Earth Engine raster API. It is beyond the scope of the current tutorials to explain the implementation of this package, but internally this package makes use of bitmap fonts which are ingested into Earth Engine as raster assets and are used to turn every character of a provided string into image glyphs, which are then translated to desired coordinates.

@@ -802,8 +775,8 @@ Note

Try commenting that event handler and observe how annotation rendering changes when you zoom in or zoom out.

-
-

8.4 Animations

+
+

7.1.4 Animations

Visualizing raster images as animations is a useful technique to explore changes in time-dependent datasets, but also, to render short animations to communicate how changing various parameters affects the resulting image—for example, varying thresholds of spectral indices resulting in different binary maps or the changing geometry of vector features.

Animations are very useful when exploring satellite imagery, as they allow viewers to quickly comprehend dynamics of changes of earth surface or atmospheric properties. Animations can also help to decide what steps should be taken next to designing a robust algorithm to extract useful information from satellite image time series. Earth Engine provides two standard ways to generate animations: as animated GIFs, and as AVI video clips. Animation can also be rendered from a sequence of images exported from Earth Engine, using numerous tools such as ffmpeg or moviepy. However, in many cases it is useful to have a way to quickly explore image collections as animation without requiring extra steps.

In this section, we will generate animations in three different ways:

@@ -993,8 +966,8 @@ Note
-
-

8.5 Terrain Visualization

+
+

7.1.5 Terrain Visualization

This section introduces several raster visualization techniques useful to visualize terrain data such as:

  • Basic hillshading and parameters (light azimuth, elevation)
  • @@ -1079,19 +1052,12 @@ Note
-
-

Synthesis

-

To synthesize what you have learned in this chapter, you can do the following assignments.

-

Assignment 1. Experiment with different color palettes from the palettes library. Try combining palettes with image opacity (using ee.Image.updateMask call) to visualize different physical features (for example, hot or cold areas using temperature and elevation).

-

Assignment 2. Render multiple text annotations when generating animations using image collection. For example, show other image properties in addition to date or image statistics generated using regional reducers for every image.

-

Assignment 3. In addition to text annotations, try blending geometry elements (lines, polygons) to highlight specific areas of rendered images.

-
-
-

Conclusion

+
+

Conclusion

In this chapter we have learned about several techniques that can greatly improve visualization and analysis of images and image collections. Using predefined palettes can help to better comprehend and communicate Earth observation data, and combining with other visualization techniques such as hillshading and annotations can help to better understand processes studied with Earth Engine. When working with image collections, it is often very helpful to analyze their properties through time by visualizing them as animations. Usually, this step helps to better understand dynamics of the changes that are stored in image collections and to develop a proper algorithm to study these changes.

-
-

References

+
+

References

Burrough PA, McDonnell RA, Lloyd CD (2015) Principles of Geographical Information Systems. Oxford University Press

Crameri F, Shephard GE, Heron PJ (2020) The misuse of colour in science communication. Nat Commun 11:1–10. https://doi.org/10.1038/s41467-020-19160-7

Imhof E (2015) Cartographic Relief Presentation. Walter de Gruyter GmbH & Co KG

@@ -1103,8 +1069,8 @@ Note

Wilkinson L (2005) The Grammar of Graphics. Springer Verlag

-
-

9 Collaborating in Earth Engine with Scripts and Assets

+
+

7.2 Collaborating in Earth Engine with Scripts and Assets

@@ -1115,16 +1081,16 @@ Chapter Information
-
-

Author

+
+

Author

Sabrina H. Szeto

-
-

Overview

+
+

Overview

Many users find themselves needing to collaborate with others in Earth Engine at some point. Students may need to work on a group project, people from different organizations might want to collaborate on research together, or people may want to share a script or an asset they created with others. This chapter will show you how to collaborate with others and share your work.

-
-

Learning Outcomes

+
+

Learning Outcomes

  • Understanding when it is important to share a script or asset.
  • Understanding what roles and permission options are available.
  • @@ -1138,24 +1104,24 @@ Chapter Information
  • Creating a script to share as a module.
-
-

Assumes you know how to:

+
+

Assumes you know how to:

  • Sign up for an Earth Engine account, open the Code Editor, and save your script (Chap. F1.0).
-
-

Introduction

+
+

Introduction

Many people find themselves needing to share a script when they encounter a problem; they wish to share the script with someone else so they can ask a question. When this occurs, sharing a link to the script often suffices. The other person can then make comments or changes before sending a new link to the modified script.

If you have included any assets from the Asset Manager in your script, you will also need to share these assets in order for your script to work for your colleague. The same goes for sharing assets to be displayed in an app.

Another common situation involves collaborating with others on a project. They may have some scripts they have written that they want to reuse or modify for the new project. Alternatively, several people might want to work on the same script together. For this situation, sharing a repository would be the best way forward; team members will be able to see who made what changes to a script and even revert to a previous version.

If you or your group members find yourselves repeatedly reusing certain functions for visualization or for part of your analysis, you could use the require module to call that function instead of having to copy and paste it into a new script each time. You could even make this function or module available to others to use via require.

Let’s get started. For this lab, you will need to work in small groups or pairs.

-