Lab: Computing Index Rasters

Introduction

There is an ever increasing availability of remotely sensed (RS) data for spatial analysis. Most recently, UAVs are being used to collect data at lower cost than ever before. However, to effectively use RS data, we need to know to georeference it, correct for sensor and atmospheric effects, and convert it into numerical values we can use.

Preparation

To save time, I've provided the following data sets for this lab. Make sure to download them ASAP.

• LandSat 7 Scene

Creating Index Rasters

Creating index rasters is done using the "Raster Math" tool. The tutorial below will walk you through how to create an NDVI raster from LandSat images. Some of the material here is reproduced from a tutorial created by Colorado State University students Grant J. Firl and Lane with Paul Evangelista and Jim Graham as advisors. The original tutorial is available at the bottom of this page and you can complete other indexes in the original tutorial for extra credit.

1. Background

Landsat data acquired from the USGS Global Visualization Viewer as in the tutorial “Acquiring Landsat Data” has been terrain-, radiometrically-, and geographically- corrected, but has been formatted to fit in an 8-bit number (ranges from 0-255). Data in such a format is referred to as “digital number,” or DN data. Before it can be used to calculate vegetation indices, the data must be converted to reflectance, a physical measurement.

2. The Data

Take a look at the data and you'll see that there are a number of files, each one ending in "BX" is a Landsat image where the X indicates the band of the spectrum in the file.

2.1 Use the provided metadata and the web to determine what area of the electromagnetic spectrum each band represents. For NDVI we want to use the "red" and "near-infrared" bands.

Question 1: What do each of the eight bands represent?

3. Preparing the data

We do not need to project the data but we do need to setup the no data values.

3.1 Use the "SetNull()" operator in raster calculator to convert the 0 values to NoData. This operator is similar to Con() but is required to set pixels in a raster to NoData (yep, the tool should be called SetNoData()). If you need help with this tool, try using the web and search for "ArcGIS SetNull". You may need to dig a bit to find a good example of how to use SetNull().

4. Convert DN to radiance data

Before converting to reflectance data, one must convert the DN data to radiance. This is done using the following expression:

L   = (gain  *  DN7) +  bias

Where L is the calculated radiance (in Watts / (sq. meter * micrometers * ster) or Watts per unit area, per solid angle, per wavelength), DN7 is the Landsat 7 ETM+ DN data, and the gain and bias are band-specific numbers.

The latest gain and bias numbers for the Landsat 7 ETM+ sensor are given in Chander et al. (2009) and are shown in the following table.

Band	Gain	Bias
1	0.778740	-6.98
2	0.798819	-7.20
3	0.621654	-5.62
4	0.639764	-5.74
5	0.126220	-1.13
7	0.043898	-0.39

 

4.1 Use the “Raster Calculator” tool in the Spatial Analyst to perform the calculation, saving the newly created radiance raster as something memorable. Be sure to use the correct gain and bias values for the band you’re working with.

Note that the calculations convert the integer data into floating point because the values in the equation have floating point numbers in them.

5. Convert radiance data to reflectance data

While radiance is the quantity actually measured by the Landsat sensors, a conversion to reflectance facilitates better comparison among different scenes. It does this by removing differences caused by the position of the sun and the differing amounts of energy output by the sun in each band. The reflectance can be thought of as a “planetary albedo,” or fraction of the sun’s energy that is reflected by the surface. It can be calculated using the following expression:

where R is the reflectance (unit-less ratio), L is the radiance calculated in step 3, d is the earth-sun distance (in astronomical units), E is the band-specific radiance emitted by the sun, and theta is the solar elevation angle. One needs three pieces of information (in addition to the radiance calculated in step 3) in order to calculate the reflectance. The first is the band-specific radiance emitted by the sun. These values are given in Chander et al. (2009) and are repeated in the following table.

 

The second and third pieces of information are d, the earth-sun distance, and theta, the solar elevation angle. These two values are dependent on the individual scene, specifically the day of the year and the time of day when the scene was captured. The solar elevation angle and the day of year are listed in the header file for each scene. This file is included with the data and ends with “_MTL.txt”. Search the file for the solar elevation angle labeled “SUN_ELEVATION”. The solar elevation angle is given in degrees.

Now take a look at the "DATE_ACQUIRED" field in the metadata. Then, take a look at the file names for the Landsat bands. The file names should contain the year the data was acquired and just after the year the "day of the year", or DOY, the data was acquired (i.e. the file names contain "YYYYDDD" where YYYY is the year and DDD is the day of the year). Check this against the date acquired and see if it makes sense (i.e. you should be able to take the number of months before the data was acquired times 30, add the day of the month, and find an approximation of the day of the year.

Then, you can find the value for from the table below.

5.1 Compute the sun elevation in radians

Keep in mind that the sine function within ArcMap requires the solar elevation angle to be in radians instead of degrees. Convert the sun angle from degrees to radians using

radians = (degrees * PI)/180°

5.3 Compute a factor

Notice that you can take the radiance out of the equation and compute just one value to multiply times the radiance. This will make your raster calculation much simpler. Compute a factor to multiply times each of the Landsat bands.

Question 2: What were your factors for bands 3 and 4?

5.4 Computing the rasters

Use the “Raster Calculator” tool in the Spatial Analyst toolbar to perform the calculation, saving the newly created reflectance raster as something memorable.

6. Enforce positive reflectance

During the conversion from DN data to reflectance, it is possible to create small negative reflectance. These values are not physical and should be set to zero. It should be noted that only very small negative numbers should be produced with this procedure. If large negative numbers are calculated, this may signify a problem with the implementation of this procedure.

6.1 Use raster calculator to check for negative values and replace them with zero.

Reflectance_NoNegatives = Con( [reflectance] < 0.0, 0.0, [reflectance])

7. Calculate NDVI

With top-of-atmosphere reflectance calculated for each band, one can now calculate vegetation indices. NDVI is given by

NDVI = (band 4 – band 3) / (band 4 + band 3)

where the reflectances are used for each band.

7.1 Use Raster Calculator to calculate the index.

7.2 Take a look at the result, does it match what you would expect? Remember to include some observations in your report.

8. Comparing with Other Data

NDVI is just one measure of vegetation health and any remotely sensed measure should be checked with other sources and/or field work.

8.1 Create a point shapefile with at least 10 points within Humboldt County. For each of the following, extract the values at the points and compare them in a table.

• The values from your NDVI layer
• The land cover type from a National Land Cover Data (NLCD) data set
• The land use value from the Humboldt County Parcels shapefile
• A description of what GoogleMaps shows in satellite view

Finally, visit at least 2 of the sites and record what you see within 15 meters (radius) of the point.

Resources

Below is a link to the original tutorial this material was taken from:

Calculating Vegetation Indices's from Landsat 5 TM and Landsat 7 ETM+ Data