1 Introduction

In this module, we will discuss the following concepts:

  1. How to rename bands of an image in GEE.
  2. How to work with premade indices.
  3. How to generate your own indices using band math.


Satellite derived image of the number of years a field has been irrigated. The most likely source of irrigation water is the Ogallala aquifer. The image is from near Holyoke, Colorado. Source: NASA

2 Background

The majority of the most productive agricultural land in the world has been extensively developed for cultivation. Efforts to increase overall agricultural production often rely on increasing the viability of marginal lands by importing water through water diversion projects or pulling ground water from subsurface reservoirs called aquifers. The Ogallala aquifer, which is found below 10 western states in the United States, provides water for approximately 30% of all irrigated agriculture in the country. The extensive use of the aquifer combined with a slow recharge rate means that many regions of the aquifer have seen significant decreases in water levels. The extraction of water from these aquifers can be so dramatic that NASA GRACE satellites have been able to detect changes in the gravitational pull of the earth over these aquifers over time. The decrease in the Ogallala is minor compared to other major aquifers around the world. Yet, there is a large amount of data available that shows agricultural land production over the Ogallala. In this module, we are going to see what remote sensing data are well suited for identifying irrigated crop lands within the Ogallala. Our goal is to apply these same methods to other regions of the world that have less publicly available agricultural data but rely on ground water extraction to promote marginal land agriculture.



Trends in the changes of ground water storage in the Earth’s 37 largest aquifers. Read more about this image here. Image credit: UC Irvine/NASA/JPL-Caltech.

2.1 About the Data

Average Precipitation

We know that irrigation is much more likely to occur in places where the crops cannot get enough water from rain-fall alone. Therefore, a relationship between the total annual precipitation and the probability of irrigated agriculture within a region should be present. To view this relationship within the Ogallala, we are going to bring in an annual precipitation raster from the WorldClim dataset, a global resource for climate data. There are 19 climate variables within WorldClim version 1 but at the moment we are only concerned with one, bio12 : annual precipitation.

Example of the WorldClim version 1 metadata from GEE.

Croplands Data Layer

We are interested in identifying what remote sensing data, can be used to distinguish irrigated cropland from the natural landscape at a location. To access the appropriateness of different remote sensing data we need a validation dataset to compare our values against. For this, we are going to rely on the USDA NASS Cropland Data layer. This dataset is developed yearly at a 30 m spatial resolution for the lower 48 of the United States. It contains 255 unique crop classes.

Note: We used an different categorical dataset for landcover in Module 3.



The USDA global cropland data layer across the contiguous United States. Yellow represents corn, green for soybeans, and red is cotton. It is also interesting to observe the diversity in California’s central valley.

3 Identifying Irrigated Lands with Remotely Sensed Imagery

You will need to open a new Google Earth Engine script for this module. Do so at code.earthengine.google.com/

3.1 Defining an Area of Interest

Our first step is going to be defining the area of interest.

Use the image of the Ogallala aquifer below to draw a bounding box within your GEE script using the draw a rectangle tool in the upper left-hand corner of the map.



An image showing the extent of the Ogallala aquifer.

This map and short description of the aquifer can be found here.

The bounding box does not need to be perfect. We just want to narrow our focus to the general area, which could be pulling from this water source. That is all the land that is above the aquifer.

3.2 Characterizing Precipitation

First, we will bring in the annual precipitation layer by calling in the BioClim data set using the unique ID. These data are stored as a multiband image so we will use the ee.Image() function to call it in.

// 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. We can pull the specific band we are interested in by using the ee.Image.select() function. Once selected, we will clip the result 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 tell from the resulting image that there is a predictable decrease in precipitation as you move from east to west across the Great Plains. It would make sense that we will see more irrigated agriculture on the western side of this region. We will add the Cropland data layer to test this assumption.

3.3 Characterizing Agriculture

In this case, the cropland data layer is stored as an image collection. We are still calling it in using a unique ID but we are using the ee.ImageCollection() function.

// Call in collection base on unique ID and print.
var crops = ee.ImageCollection("USDA/NASS/CDL");
print(crops, "crops");

Image collections contain multiple images with the same attributes. As a result, working with image collections relies heavily on filtering functions. Filtering is useful but applying a filter does not change the data type of the imagery. Filters always return an image collection object. This is fine, but different GEE objects have different functions. Often you will need to convert an image collection to an image. One of the most common ways to do so is by applying a reducer. Reducers are one of the elements that make GEE stand out when compared to other remote sensing software. While they are very powerful, we will not be covering their functionality in detail 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 the 2016 cropland dataset.

Rather than 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 called the image, selected the band of interest, and clipped the image to our output area.

We can use these two datasets to look at big picture patterns of where land is cultivated and the average precipitation of that area. If you really want to dig into these data, you can look at the types of crops grown. You will need to pull up the 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 get overwhelming in a hurry. Below is a table to the common crops in our area of interest.

Value Color Crop
1 Yellow Corn
2 Red Cotton
24 Brown Winter Wheat

Circular areas of a crop indicate fields using center pivot irrigation. Considering how arid this region is and how few major rivers run through it, it is most likely that water used by the center pivots is coming from below the surface, namely, the Ogallala aquifer. Fields that are not round can still be irrigated land. While potentially less efficient, many areas use flood irrigation to water crops. Determining if a field is flood irrigated or not based on these two layers would be difficult without extensive knowledge of the crops and the growing conditions of the place. We are hoping that the difference in the spectral reflectance values between the crop land and the non-crop lands will help us make that distinction.

You can turn both layers on and adjust the transparency to see the overlap between precipitation and crop lands.

By increasing the transparency of the precipitation image, we can visually assess relationships between crop type and annual precipitation over a large geographic area.

3.4 Visualize Existing Indices

One of the major benefits of GEE is that it is a data repository. In this section we are going to access a few preexisting datasets that have been generated to see if we can rely on them to answer our question and not have to worry about taking the time to generate our own. No need to reinvent the wheel. We are going to look at three options:

Evapotranspiration: the process by which water is transferred from the land to the atmosphere by evaporation from the soil and other surfaces and through transpiration from plants.

Gross Primary Productivity: the rate at which photosynthesis or chemosynthesis 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 in an area of land.

These indices give us information about how much plant life is present at a given location. Because irrigation is supporting plant growth by providing more water than is otherwise available, we may suspect that irrigated land will stand out in the images derived from the indices. You can find details on all these datasets by searching for them in GEE using the following names.

Evapotranspiration: Terra Net Evapotranspiration 8-Day Global 500m

Gross Primary Productivity: Landsat Gross Primary Production 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 is good practice to call in a new dataset and print it so you can look at the data structures and image availability before trying to filter. GEE will apply filters to datasets even if the request does not match the information present in the dataset, so it is up to you to understand if what you are asking is possible.

From the metadata we established a date range that makes sense and we know what bands we are 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 collection. We are applying a median reducer here because we want generalized values that are representative of the growing season in this region. The median reducer will take a stack of images, calculate the median value of each cell, and compile all those cells back into a single image. Once you apply a median reducer, your image no longer represents a single point in time but a collection of information over a period of time.

We will repeat this process for the Gross Primary Productivity 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 derived from MODIS imagery. The NDVI layer that we are using is from a different sensor, Landsat. While there are many differences between the two sensors, the big factors are spatial resolution (pixel area) and temporal resolution (image recapture period). Landsat has a much finer spatial resolution (30 m compared to MODIS 500 m) but has a much lower temporal resolution (16 days return time compared to the 1 day of MODIS). Comparing multiple sensors is always recommended because each one has benefits and pitfalls. It is up to you to understand what tool will work 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);

With the three images loaded on the map take some time to explore the visual patterns identified across the landscape. A great place to see the difference in reflectance between the images is around Joes, Colorado. You can find it by searching for “Joes, Colorado” in the search bar as you would in Google Maps. There are some obvious center pivot irrigation areas, undeveloped land that is representative of the natural landscape and riparian zones all within a few mile radius of Joes.

We added all images to the map with defined visualizations parameters of min and max. These values were created by applying a stretch and recording the min and max values in the script. This is a method for making the differences in the data easier to see and does not alter the underlying values within the data.

An NDVI image of the area around Joes, Colorado. Bright areas represent locations with high NDVI values and are associated with irrigated lands and riparian areas.

3.5 Thoughts and Limitations

What indices you decide to use is really up to you. The three used here seem to capture difference in the irrigated and non-irrigated lands but there are some unique considerations for all the images.

  1. Evapotranspiration: Cell size is a little big for distinguishing center pivot locations.

  2. Gross Primary Productivity: Similar issue with cell size as the evapotranspiration data. Also, this dataset is only available for the contiguous United States, so it will not help us with global questions.

  3. NDVI: Seems like a winner but data are only available until 2017, so it does not help us with anything after that time.

Based on these observed limitations, we are going to generate our own set of indices using Landsat 8 data that can be applied to any location at any point in time that the imagery has been collected.

4 Building Your Own Indices

We are going to create NDVI indices on Landsat 8 imagery because it really highlights vegetated areas well.

The first step is bringing 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 with Landsat imagery and effectively allows you to directly compare values from different images by mitigating the atmospheric effects. You can find a detailed description of 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');

The filter by metadata is a new command. You can see how it works by checking the documentation of the imagery and the function itself.

There are multiple ways to define indices. We will cover three in this section. Two of these methods require defining individual bands and one requires selecting bands from the image. Our first step is defining individual images 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 measure of vegetation. It is calculated using the following equation:

NDVI = (nir - red)/(nir + red)

The first method to generate an index involves using the math functions to reproduce the NDVI equation. You cannot use the symbols (+,-,/,*) as mathematic expressions on their own in GEE. Instead, you will need to write out the operations.

// 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 relying on built-in functions. In our case there is a normalizedDifference() function for ee.Image features. We applied a reducer to our Landsat 8 image collection so we are left with an image that we can use. Therefore, calling this function is as simple as pointing it to the correct bands.

// 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 last method available for generating indices is building an expression. Given that there is a normalized difference function built into GEE, it does not make sense to use this method for NDVI. We use it here as an example. The process is necessary for more complex indices such as Tassel Cap Transformations or the 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);
Expressions are much more complicated. We define the mathematical expression within a string. Then we use a dictionary to reference what the variables in the expression are referring to. We saved the near inferred (‘B5’) and red band (‘B4’) as variables already so we can call them directly.

Identical NDVI values from the three different methods for generating the index The image is centered over the town Ogallala, Nebraska. The dark feature in the center is McConaughey Lake.

If you compare the three NDVI images you will realize that they are all the same; different routes to the same end point. What method you use will depend on the datasets you have and what index you are attempting to create. There are hundreds of options out there. Index database is a great reference for what is possible with each specific sensor indexdatabase.

4.1 Changing Location

Looking back at the change in ground water map that was created using the GRACE sensors (see image below), let us quickly apply our indices to a new region with significant depletion of groundwater resources and see if we can spot some irrigated agriculture.



Areas of red represent measured decrease in aquifers over the last 20 years.

To change the area of interest you can either.

  1. Delete the current geometry feature and draw a new one.

or

  1. Add a second feature to the existing geometry layer.

We will go with the second.

Select the geometry layer in the geometry imports drop down box. Once selected the features name should appear bold. You can then draw a new box over the Arabian Peninsula and click run.



New geometry feature added over the Arabian Peninsula, an area that is seeing substantial declines in groundwater availability.

The result is pretty obvious. Large sections of the landscape have been converted to center pivot irrigation agriculture. In a desert ecosystem there is little doubt that the water for this agriculture is coming from ground water. Unfortunately, in this situation, it is ground water that is not being replenished. As with all non-renewable resources, this time of abundance will not last.



A very large area of irrigated agriculture identified used remotely sensed data. This image is near the city Al Huwayd, Saudi Arabia.

5 Conclusion

Through this module, we conducted a comparison of multiple datasets that relate to irrigated lands. The concepts involved with the filtering and manipulation of the various images and images collections can be applied to a wide verity of processes. The process of utilizing or generating your own indices is the first step to conducting more quantitative analysis within GEE.

5.1 Final Code for module

// 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);