1Introduction
In this module we will discuss the following concepts:
- How to rename bands of an image in GEE.
- How to use existing remote sensing indices.
- How to use band math to generate your own remote sensing index.
Satellite image of the number of years a field has been irrigated. The most likely source of irrigation water is the Ogallala Aquifer. Image from near Holyoke, Colorado. Source: NASA
2Background
Much of the world’s most productive agricultural land is extensively developed for farming. Efforts to increase overall agricultural production often rely on improving the viability of marginal lands by importing water through diversion projects or extracting groundwater from underground reservoirs called aquifers. The Ogallala Aquifer lies beneath 10 states in the western United States and provides water for about 30 percent of the country’s irrigated agriculture. Extensive use of the aquifer coupled with slow recharge rates means water levels have dropped significantly in many areas of the aquifer. The process of extracting water from these aquifers is so dramatic that the NASA GRACE satellite has been able to detect changes in Earth’s gravitational pull on these aquifers over time. The Ogallala’s decline is small compared to other major aquifers around the world. However, there is a wealth of data available showing Ogallala’s agricultural land yields. In this module we will learn which remote sensing data are well suited for identifying irrigated farmland within the Ogallala. We aim to apply these same methods to other regions of the world that have less publicly available agricultural data but rely on groundwater extraction to promote agriculture on marginal lands.
Trends in groundwater storage in Earth’s 37 largest aquifers. Read more about this image here. Image credit: University of California, Irvine/NASA/JPL-Caltech.
2.1 About data
Average precipitation
We know that irrigation is more likely to occur where crops do not receive enough moisture from rainfall. Therefore, there should be a relationship between total annual precipitation and the probability of irrigated agriculture within a region. To see this relationship in Ogallala, we will import an annual precipitation raster from the WorldClim dataset, a global resource for climate data. There are 19 climate variables in WorldClim version 1, but currently we are only interested in one, bio12: annual precipitation.
WorldClim version 1 metadata example forGEE.
Farmland Data Layer
We are interested in determining what remote sensing data can be used to distinguish irrigated farmland from the natural landscape at a location. In order to access the appropriateness of different remote sensing data, we need a validation dataset to compare our values. For this we will rely on the USDA NASS Cropland data layer. This dataset is developed annually with a spatial resolution of 30 m for the 48 US states. It contains 255 unique crop categories.
Note: We used a different classification dataset for land cover in Module 3.
TheUSDA Global Farmland Data Layer spans the continental United States. Yellow represents corn, green represents soybeans, and red represents cotton. It’s also interesting to observe the diversity of California’s Central Valley.
3 Use remote sensing images to identify irrigated land
You need to open a new Google Earth Engine script for this module. Do this at code.earthengine.google.com/
3.1 Define a region of interest
Our first step is to define the region of interest.
Using the Ogallala Aquifer image below, draw the bounding box in GEE script using the Draw Rectangle tool in the upper left corner of the map.
Image showing the extent of the Ogallala Aquifer.
This map and a brief description of the aquifer can be found here.
The bounding box doesn’t need to be perfect. We’re just trying to narrow our focus to the general area where this water source might come from. That’s all the land above the aquifer.
3.2 Characterizing precipitation
First, we will introduce the annual precipitation layer by calling the BioClim dataset with a unique ID. This data is stored as a multiband image, so we will call it using the ee.Image() function.
// Call in image by unique ID.
var bioClim = ee.Image("WORLDCLIM/V1/BIO");
print(bioClim, "bioClim");
From the print statement we can see that each variable is stored as a band. ee.Image.select()We can use this function to pull the specific band we are interested in. Once selected, we will clip the results to our area of interest and rename it so we remember what “bio12” stands for.
// Select specific band, clip to region of interest, and rename.
var precip = bioClim
.select("bio12")
.clip(geometry)
.rename("precip");
Map.addLayer(precip, {<!-- -->min:240, max:1300}, "precip", false);
We can see from the resulting images that there is a predictable decrease in precipitation as you travel across the Great Plains from east to west. It makes sense that we would see more irrigated agriculture in the western part of the region. We’ll add a Cropland data layer to test this hypothesis.
3.3 Characterizing Agriculture
In this case, the farmland data layer is stored as a collection of images. We’re still calling it with a unique ID, but we’re using the ee.ImageCollection() function.
// Call in collection base on unique ID and print.
var crops = ee.ImageCollection("USDA/NASS/CDL");
print(crops, "crops");
An image collection contains multiple images with the same properties. Therefore, processing image collections relies heavily on filtering functions. Filtering is useful, but applying a filter does not change the image’s data type. Filters always return an image collection object. This is fine, but different GEE objects have different functions. Often you need to convert a collection of images into an image. One of the most common methods is to apply a reducer. Reducers are one of the elements that make GEE stand out compared to other remote sensing software. Although they are very powerful, we will not go into detail about their capabilities in this module.
// Filter collection by a date range and pull a specific image from the collection.
var crop2016 = crops
.filterDate("2016-01-01","2016-12-30")
.select(['cropland']);
print(crop2016, "crop2016");
As we expected, the crop2016 object is still an image collection object, even though it only contains a single image.
Print statement for the2016 farmland dataset.
Instead of applying a reducer, we can simply call this image by the unique ID it contains and then use the Image.select() function to pull the specific band we are interested in. We can view the image ID by printing the object and reading the metadata.
// Call the specific image by the unique ID, select band of interest, and clip to region of interest.
var cropType2016 = ee.Image("USDA/NASS/CDL/2016")
.select("cropland")
.clip(geometry);
print(cropType2016, "cropType2016");
Map.addLayer(cropType2016, {<!-- -->}, "cropType2016", false);
Here we load the image, select the band of interest, and clip the image to our output region.
We can use these two datasets to see big picture patterns in land cropland and average precipitation in the area. If you really want to dig into the data, you can look at the types of crops grown. You need to extract metadata from the USDA layer by searching the dataset. As mentioned above, there are 255 different land cover types in the dataset, so it can become overwhelming very quickly. Below is a table of common crops in our area of interest.
| Value | Color | Crop |
|---|---|---|
| 1 | Yellow | Corn |
| 2 | Red | Cotton |
| 24 | Brown | Winter wheat |
The circular area of crops represents a field irrigated using a center pivot. Given how arid the area is and the number of major rivers that flow through it, the water used by Center Pivot likely comes from below the surface, known as the Ogallala Aquifer. Fields that are not round can still be irrigated. Although it may be less efficient, flood irrigation is used in many areas to water crops. It is difficult to determine whether a field is flood-irrigated based on these two layers without extensive knowledge of the crop and growing conditions on the site. We hope that the difference in spectral reflectance values between cropland and non-cropland will help us make the distinction.
You can turn on both layers and adjust the transparency to see the overlap between precipitation and farmland.
By increasing the transparency of precipitation images, we can visually assess the relationship between crop types and annual precipitation over large geographic areas.
3.4 Visualizing existing indices
One of the main advantages of GEE is that it is a data repository. In this section, we’ll access some pre-existing datasets that have been generated and see if we can rely on them to answer our questions without worrying about spending time generating our own datasets. No need to reinvent the wheel. We’ll look at three options:
Evaporation Transpiration: The process by which water is transferred from land to the atmosphere through evaporation from soil and other surfaces and transpiration from plants.
Primary productivity: The rate at which photosynthesis or chemical synthesis occurs.
Normalized Difference Vegetation Index: A dimensionless index that describes the difference between visible and near-infrared reflectance of vegetation cover and can be used to estimate the density and health of vegetation within a land area.
These indices give us information about how much plant life is present at a given location. Since irrigation supports plant growth by providing more water than would otherwise be available, we might suspect that irrigated land would stand out in the exponential image. You can search the details of all these datasets in GEE using the following names.
Evaporation: Terra Net evaporation 500m globally in 8 days
Gross primary productivity: Landsat gross primary productivity CONUS
NDVI : Landsat 8 32-day NDVI composite
// Evapotranspiration: call in image collection by unique ID and print to understand data structure.
var et = ee.ImageCollection("MODIS/006/MOD16A2")
.filterBounds(geometry);
print(et, "modisET");
It’s a good practice to call up a new dataset and print it so you can see the data structure and image availability before trying to filter. GEE applies filters to the dataset even if the request does not match information present in the dataset, so you need to understand whether your request is feasible.
From the metadata, we build a date range that makes sense, and we know which bands we’re interested in. Now we will filter the collection.
// Filter by data, select band, apply reducer to get single image and clip image.
var et2 = et
.filterDate("2016-07-01","2016-09-30")
.select(['ET'])
.median()
.clip(geometry);
Map.addLayer(et2, {<!-- -->max: 357.18 , min: 29.82}, "evapotransparation", false);
The evapotranspiration layer is provided as an image set. We apply a median reducer here because we want generalized values that represent the growing season in the region. The median reducer will take a bunch of images, calculate the median of each cell, and compile all of those cells back into a single image. After applying the median reducer, your image no longer represents a single point in time, but rather a collection of information over time.
We will repeat this process for the “Gross Primary Productivity Product” layer as well.
// Gross Primary Productivity.
var GPP = ee.ImageCollection("UMT/NTSG/v2/MODIS/GPP");
print(GPP, "modisGPP");
var GPP2 = GPP
.filterDate("2016-07-01","2016-09-30")
.select(['GPP'])
.median()
.clip(geometry);
Map.addLayer(GPP2, {<!-- -->max: 1038.5 , min: 174.5}, "gross Primary Productivity", false);
The two previous datasets were from MODIS images. The NDVI layer we use is from a different sensor, Landsat. While there are many differences between the two sensors, the main factors are spatial resolution (pixel area) and temporal resolution (image re-capture period). The spatial resolution of Landsat is much better (30 m compared to 500 m for MODIS), but the temporal resolution is much lower (return time is 16 days compared to 1 day for MODIS). It is always recommended to compare multiple sensors as each sensor has advantages and disadvantages. You need to know which tool is best for your application.
// Pulling in NDVI data.
var LS8NDVI = ee.ImageCollection("LANDSAT/LC8_L1T_32DAY_NDVI")
print(LS8NDVI, "ls8NDVI");
// Filter, reduce, then clip.
var LS8NDVI2 = LS8NDVI
.filterDate("2016-07-01","2016-09-30")
.median()
.clip(geometry);
Map.addLayer(LS8NDVI2, {<!-- -->max: 0.66 , min: -0.23}, "NDVI", false);
Loading the three images on the map will take some time to explore the visual patterns identified across the landscape. A good place to observe the difference in reflectivity between images is nearJoseph, Colorado. You can find it by searching for “Joes, Colorado” in the search bar just like you would in Google Maps. Within a few miles of Jos there are some distinct central hub irrigation areas, undeveloped land and riparian zones that represent natural landscapes.
We add all images to the map using the defined visualization parameters min and max. These values are created by applying a stretch and recording the minimum and maximum values in the script. This is a way to make differences in the data easier to see without changing the underlying values in the data.
NDVI image of the area surroundingJos, Colorado. Bright areas represent locations with high NDVI values and are associated with irrigated lands and riparian areas.
3.5 Ideas and Limitations
What index you decide to use is entirely up to you. The three used here seem to capture the differences between irrigated and non-irrigated land, but all images have some unique considerations.
-
Evapotranspiration: Cell size is a bit large for distinguishing center pivot locations.
-
Total primary productivity: Similar cell size issues to evaporation data. Additionally, this dataset only applies to the contiguous United States, so it doesn’t help us solve global problems.
-
NDVI: Looks like the winner, but the data is only available until 2017, so it doesn’t help us after that.
Based on these observed limitations, we will use Landsat 8 data to generate our own set of indices that can be applied to any location at any point in time of image collection.
4Create your own index
We will create the NDVI index on the Landsat 8 image as it highlights the vegetated areas really well.
The first step is to bring in Landsat imagery. Because we are working over a large area and compiling images from multiple time periods, we will use the preprocessed USGS Landsat 8 Surface Reflectance Tier 1 dataset. This is the highest level of preprocessing available for Landsat imagery, effectively allowing you to directly compare values from different images by mitigating atmospheric effects. You can find detailed instructions for the various preprocessing levels in Module 5.
// Call in Landsat 8 surface reflectance Image Collection, filter by region and cloud cover, and reduce to single image.
var LS8_SR2 = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
.filterDate("2017-07-01","2017-09-30")
.filterBounds(geometry)
.filterMetadata('CLOUD_COVER', 'less_than', 20)
.median();
print(LS8_SR2,'LS8_SR2');
Filter by metadata is a new command. You can find out how it works by checking the image’s documentation and the function itself.
There are several ways to define indexes. We’ll cover three in this section. Two of the methods require defining individual bands and one requires selecting bands from the image. Our first step is to define a single image from the bands.
// Pull bands from the LS8 image and save them as images.
var red = LS8_SR2.select('B4').rename("red");
var nir = LS8_SR2.select('B5').rename("nir");
print(red);
The Normalized Difference Vegetation Index (NDVI) is a very common vegetation measurement method. Calculate using the following equation:
NDVI = (nir – red)/(nir + red)
The first way to generate an exponent is to use a mathematical function to reproduce the NDVI equation. You cannot use symbols ( + ,-,/,*) alone as mathematical expressions in GEE. Instead, you need to write out the operation.
// Generate NDVI based on predefined bands.
var ndvi1 = nir.subtract(red).divide(nir.add(red)).rename('ndvi');
Map.addLayer(ndvi1, {<!-- -->max: 0.66 , min: -0.23},'NDVI1',false);
Calling math functions directly can quickly become cumbersome. Another option is to rely on built-in functions. In our example, there is a feature normalizedDifference() function ee.Image. We applied a reducer on a collection of Landsat 8 images, so we got images that we could use. So calling this function is as simple as pointing it to the correct band.
// Define LS8 NDVI using built in normalized difference function.
var ndvi2 = LS8_SR2.normalizedDifference(['B5', 'B4']);
Map.addLayer(ndvi2, {<!-- -->max: 0.66 , min: -0.23},'NDVI2',false);
The final method you can use to generate an index is to build an expression. Given that the normalized difference function is built into GEE, it does not make sense to use this approach with NDVI. We use it as an example here. This process is necessary for more complex indices such as fringed cap conversion or enhanced vegetation index (EVI).
// Use the expression function to generate NDVI.
var ndvi3 = LS8_SR2.expression("(nir - red) / (nir + red)", {<!-- -->
"nir" : nir,
"red" : red }
)
Map.addLayer(ndvi3, {<!-- -->max: 0.66 , min: -0.23},'NDVI3',false);
The expression is much more complex. We define mathematical expressions in strings. We then use the dictionary to reference what the variable in the expression refers to. We have saved the approximate inference (B5’) and the red band (B4’) as variables so we can call them directly.
Identical NDVI values for three different methods of generating indexes This image is centered on the town of Ogallala, Nebraska. The dark feature in the center is Lake McConaughey.
If you compare the three NDVI images, you’ll see that they are all the same; different routes to the same end point. What method you use depends on the dataset you have and the index you are trying to create. There are hundreds of options out there. The index database is a good reference to understand the possibilities of indexing the database for each specific sensor.
4.1 Change location
Looking back at the changes to the groundwater map created using the GRACE sensor (see image below), let’s quickly apply our index to new areas with severely depleted groundwater resources and see if we can uncover some irrigated agriculture.
Thered areas represent measured decreases in the aquifer over the past 20 years.
To change the area of interest, you can too.
- Delete the current geometric feature and draw a new one.
or
- Adds a second feature to the existing geometry layer.
We’ll choose the second one.
Select the geometry layer from the Geometry Import drop-down box. Once selected, the feature name should appear in bold. You can then draw a new box on Arabian Peninsula and click run.
New geometric features have been added to the Arabian Peninsula, a region where groundwater supplies are declining significantly.
The results are clear. Much of the landscape has been converted to central hub irrigated agriculture. In desert ecosystems, there is no doubt that this agricultural water comes from groundwater. Unfortunately, in this case, the groundwater was not recharged. As with all non-renewable resources, this period of abundance will not last.
Identified large areas of irrigated agriculture using remote sensing data. This image is located near the city of Al Huwayd in Saudi Arabia.
5Conclusion
Through this module, we performed a comparison of multiple data sets related to irrigated land. Concepts involving filtering and manipulation of various images and image collections can be applied to a wide range of processes. The process of utilizing or generating your own indices is a first step toward more quantitative analysis within GEE.
5.1 module final code
// Call in image by unique ID.
var bioClim = ee.Image("WORLDCLIM/V1/BIO");
print(bioClim, "bioClim");
// Select specific band, clip to region of interest, and rename.
var precip = bioClim
.select("bio12")
.clip(geometry)
.rename("precip");
Map.addLayer(precip, {<!-- -->min:240, max:1300}, "precip", false);
// Call in collection base on unique ID and print.
var crops = ee.ImageCollection("USDA/NASS/CDL");
print(crops, "crops");
// Filter collection by a date range and pull a specific image from the collection.
var crop2016 = crops
.filterDate("2016-01-01","2016-12-30")
.select(['cropland']);
print(crop2016, "crop2016");
// Call the specific image by the unique ID, select band of interest, and clip to region of interest.
var cropType2016 = ee.Image("USDA/NASS/CDL/2016")
.select("cropland")
.clip(geometry);
print(cropType2016, "cropType2016");
Map.addLayer(cropType2016, {<!-- -->}, "cropType2016", false);
// Evapotranspiration: call in image collection by unique ID and print to understand data structure.
var et = ee.ImageCollection("MODIS/006/MOD16A2")
.filterBounds(geometry);
print(et, "modisET");
// Filter by data, select band, apply reducer to get single image and clip image.
var et2 = et
.filterDate("2016-07-01","2016-09-30")
.select(['ET'])
.median()
.clip(geometry);
Map.addLayer(et2, {<!-- -->max: 357.18 , min: 29.82}, "evapotransparation", false);
// Gross Primary Productivity.
var GPP = ee.ImageCollection("UMT/NTSG/v2/MODIS/GPP");
print(GPP, "modisGPP");
var GPP2 = GPP
.filterDate("2016-07-01","2016-09-30")
.select(['GPP'])
.median()
.clip(geometry);
Map.addLayer(GPP2, {<!-- -->max: 1038.5 , min: 174.5}, "gross Primary Productivity", false);
// Pulling in NDVI data.
var LS8NDVI = ee.ImageCollection("LANDSAT/LC8_L1T_32DAY_NDVI");
print(LS8NDVI, "ls8NDVI");
// Filter, reduce, then clip.
var LS8NDVI2 = LS8NDVI
.filterDate("2016-07-01","2016-09-30")
.median()
.clip(geometry);
Map.addLayer(LS8NDVI2, {max: 0.66 , min: -0.23}, "NDVI", false);
// Call in Landsat 8 surface reflectance Image Collection, filter by region and cloud cover, and reduce to single image.
var LS8_SR2 = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
.filterDate("2017-07-01","2017-09-30")
.filterBounds(geometry)
.filterMetadata('CLOUD_COVER', 'less_than', 20)
.median();
print(LS8_SR2,'LS8_SR2');
// Pull bands from the LS8 image and save them as images.
var red = LS8_SR2.select('B4').rename("red");
var nir = LS8_SR2.select('B5').rename("nir");
print(red);
// Generate NDVI based on predefined bands.
var ndvi1 = nir.subtract(red).divide(nir.add(red)).rename('ndvi');
Map.addLayer(ndvi1, {<!-- -->max: 0.66 , min: -0.23},'NDVI1',false);
// Define LS8 NDVI using built in normalized difference function.
var ndvi2 = LS8_SR2.normalizedDifference(['B5', 'B4']);
Map.addLayer(ndvi2, {<!-- -->max: 0.66 , min: -0.23},'NDVI2',false);
// Use the expression function to generate NDVI.
var ndvi3 = LS8_SR2.expression("(nir - red) / (nir + red)", {<!-- -->
"nir" : nir,
"red" : red }
)
Map.addLayer(ndvi3, {<!-- -->max: 0.66 , min: -0.23},'NDVI3',false);