Created By: Grant J.

Firl, Lane Carter

Advisors: Paul Evangelista, Jim Graham
Date: January 2011
Software: ArcGIS v9.2, Windows Notepad

Lesson 10: Calculating Vegetation Indices from Landsat 5 TM and
Landsat 7 ETM+ Data

Vegetation indices like the Normalized Difference Vegetation Index (NDVI), Soil-
Adjusted Vegetation Index (SAVI), and Tasseled Cap transformations are useful
measures of vegetation calculated from remotely sensed data. It is relatively
straightforward to calculate these indices using data from the Landsat sensors and
ArcMap. In this tutorial, a method for doing so is described.

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. In addition, if one is interested in calculating the tasseled cap brightness,
greenness, and wetness indices, further transformation is needed. As described by Huang
et al. (2002), the tasseled cap coefficients for the Landsat 5 TM data were developed for
both DN format data (Crist and Cicone, 1984) and for surface reflectance data (Crist,
1985). Reflectances calculated from Landsat data are a so-called top of atmosphere
(TOA) measurement, however, so the tasseled cap coefficients calculated in the two
papers from the 1980s are not applicable to Landsat data unless atmospheric correction is
applied, a process that may introduce considerable error. Tasseled cap coefficients are
calculated for the TOA reflectance data from the Landsat 7 ETM+ sensor by Huang et al
(2002). These coefficients are directly applicable for Landsat 7 ETM+ TOA reflectance
data, and can be used with Landsat 5 TM data using a further transformation described in
Vogelmann et al. (2001).

Process:
1. Reclassify the Landsat data so that all 0 values are mapped to NoData using the
method described in step 3 of the Mosaicking and Clipping Landsat Data. We
do not want to calculate reflectances or vegetation indices on the sections where
data is missing.


Figure 1. Original Landsat 5 TM DN data with state and county boundaries
overlaid


Figure 2. Reclassified Landsat 5 TM data

2. Convert Landsat 5 TM data to the equivalent Landsat 7 ETM+ data
a. If using Landsat 7 ETM+ data, go to step 3.
b. In order to be able to use the tasseled cap coefficients developed by
Huang et al. (2002) for the Landsat 5 TM sensor, one must convert the
Landsat 5 TM DN data into data that is equivalent to data recorded by
the Landsat 7 ETM+ sensor (because the two sensors have slightly
different calibration). This process is described by Vogelmann et al.
(2001) in reverse; that is, they converted from Landsat 7 ETM+ data to
Landsat 5 TM equivalent. To convert from Landsat 5 TM DN data to
Landsat 7 ETM+ DN data, we use the following expression:

DN7 = slope

* DN5 ( )+ intercept

where DN7 is the Landsat 7 ETM+ equivalent DN data, DN5 is the
Landsat 5 TM DN data, and the slope and intercept are band-specific
numbers given by the inverse of those found in Vogelmann et al. (2001).
The needed values are given in the following table (the inverse has already
been performed).

Band Slope Intercept
1 0.943 4.21
2 1.776 2.58
3 1.538 2.50
4 1.427 4.80
5 0.984 6.96
7 1.304 5.76

Use the Raster Calculator tool in the Spatial Analyst toolbar to perform
the calculation.
After the Raster Calculator performs the provided evaluation, it may be
wise to save each new evaluated raster by exporting the data. To do this
right click on the layer within the table of Contents and select Export
Data. When the window below is visible, change the data format to a
GRID and be sure to browse to the correct workspace while naming the
raster something memorable.

While calculating, be sure to use the correct slope and intercept values for
the band youre working with. Repeat this process for all bands needed for
the tasseled cap indices (1-5 and 7).


Figure 3. Raster calculator function to perform the transformation from Landsat 5 TM
DN data to the equivalent Landsat 7 ETM+ DN data


Figure 4. Resulting equivalent Landsat 7 ETM+ DN data

IMPORTANT: Be sure to save the calculations if needed

3. Convert DN data to radiance data
a. 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 * m * ster)],

DN7 is the Landsat 7 ETM+ DN data (or the equivalent calculated in step
2), 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

Use the Raster Calculator tool in the Spatial Analyst toolbar 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
youre working with. Repeat this process for all bands needed for the
tasseled cap indices (1-5 and 7).

The process of this calculation and following calculations may take a
significant amount of time so be sure to save any progress at intervals
along the way.


Figure 5. Raster Calculator function to perform radiance calculation

Figure 6. Radiance data resulting from the Raster Calculator calculation

4. Convert radiance data to reflectance data
a. 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 suns energy that is reflected by the surface.
It can be calculated using the following expression:

R

=
t * L

*d
2
E
sun,
*sin(u
SE
)


where R

is the reflectance (unitless ratio), L

is the radiance calculated in

step 3, d is the earth-sun distance (in astronomical units), E
sun,
is the
band-specific radiance emitted by the sun, and u
SE
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.







Band E
sun,
[Watts / (sq. meter * m)]
1 1997
2 1812
3 1533
4 1039
5 230.8
7 84.9

The second and third pieces of information are d, the earth-sun distance, and u
SE,
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 and the day of the year labeled
DATE_HOUR_CONTACT_PERIOD. The solar elevation angle is given in degrees
and the date is in the format YYDDDHH where the 3 D digits denote the day of the
year. For example, 0624117 means the 241
st
day of 2006 at 17 UTC. Once the day of
the year is acquired, use the table reproduced from Chander et al. (2009) to find the earth-
sun distance. For example, for day 241, the earth-sun distance is 1.00992 astronomical
units.



Figure 7. An example of where the solar elevation angle is found within the header file.

Figure 8. An example of where the day of the year is found within the header file



Figure 9. Table from Chander et al. (2009) showing the earth-sun distance as a function
of day of the year









Use the Raster Calculator tool in the Spatial Analyst toolbar to perform
the calculation, saving the newly created reflectance raster as something
memorable. Be sure to use the correct E
sun,
for the band youre working
with and the correct earth-sun distance and solar elevation angle for the
scene. Keep in mind that the sine function within ArcMap requires the
solar elevation angle to be in radians instead of degrees. Convert from
degrees to radians using

radians = (degrees * t)/180

Below, this conversion has been incorporated into the calculation. Repeat
this process for all bands needed for the tasseled cap indices (1-5 and 7). If
an error occurs while the Raster Calculator attempts to evaluate this
expression, it may be a simple syntax error. Before doing anything
extravagant to fix the problem, try counting parentheses () to make sure
that all are present. Another remedy may include adding parentheses ()
around the numerator and denominator.


Figure 10. Raster Calculator function for calculating the reflectance

Figure 11. Reflectance data resulting from the Raster Calculator function

5. Enforce positive reflectances

a. During the conversion from DN data to reflectance, it is possible to
create small negative reflectances. 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. Use raster calculator to check for
negative values and replace them with zero.

corrected_reflectance = CON( [reflectance] < 0.0, 0.0, [reflectance])

The CON () statement can either be utilized by directly typing it into the
Raster Calculator or by clicking the appropriate function button within
certain versions of ArcMap.

Figure 12. Raster Calculator function for enforcing positive reflectance values

6. Calculate NDVI
a. With TOA reflectances 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 TOA reflectances are used for each band. Use Raster Calculator
to calculate the index.

Figure 13. Raster Calculator function for calculating NDVI from reflectance data

Figure 14. NDVI calculated from the Raster Calculator calculation

7. Calculate SAVI
a. SAVI is similar to NDVI, but attempts to account for the background
soil conditions. SAVI is given by

SAVI = (1+ L)*(band 4 band 3) / (band 4 + band 3 +L)

where the TOA reflectances are used for each band and L is a soil
brightness correction factor. From Huete (1988), L=0.5 is used in most
conditions. Use Raster Calculator to calculate the index.

Figure 15. Raster Calculator function for calculating SAVI from reflectance data


Figure 16. SAVI calculated from the Raster Calculator expression

8. Calculate tasseled cap indices
a. Tasseled cap indices give a measure of the brightness, greenness, or
wetness of a pixel and utilize a linear combination of 6 of Landsats
frequency bands. Each of the brightness, greenness, and wetness
indices can be calculated using the expression

tas. cap
i
= coeff
1
*band
1
( )+ coeff
2
*band
2
( )+ coeff
3
*band
3
( )
+ coeff
4
*band
4
( )+ coeff
5
*band
5
( )+ coeff
7
*band
7
( )

where tas. cap
i
is the calculated tasseled cap index for brightness,
greenness, or wetness depending on the coefficients used, the bands are
the TOA reflectances calculated in this tutorial, and the coefficients are
given by Huang et al. (2002):

Index Band 1 Band 2 Band 3 Band 4 Band 5 Band 7
Brightness 0.3561 0.3972 0.3904 0.6966 0.2286 0.1596
Greenness -0.3344 -0.3544 -0.4556 0.6966 -0.0242 -0.2630
Wetness 0.2626 0.2141 0.0926 0.0656 -0.7629 -0.5388

Use Raster Calculator to calculate the desired index. Outputs will be
displayed with a black-white color ramp. The color ramp can be changed
within the layers properties like Figures 18-20.


Figure 17. Raster Calculator function for calculating the tasseled cap greenness index

Figure 18. Tasseled cap brightness index


Figure 19. Tasseled cap greenness index

Figure 20. Tasseled cap wetness index

References:

Chander, G., B. L. Markham, and D. L. Helder, 2009: Summary of current radiometric
calibration coefficients for Landsat MSS, TM, ETM+, and EO-1 ALI sensors.
Remote Sensing of Environment, 113, 893-903.

Crist, E. P., and R. C. Cicone, 1984: A physically-based transformation of Thematic
Mapper data the TM Tasseled Cap. IEEE Trans. On Geosciences and Remote
Sensing, GE-22, 256-263.

Crist, E. P., 1985: A TM tasseled cap equivalent transformation for reflectance factor
data. Remote Sensing of Environment, 17, 301-306.

Huang, C., B. Wylie, L. Yang, C. Homer, and G. Zylstra, 2002: Derivation of a tasseled
cap transformation based on Landsat 7 at-satellite reflectance. International
Journal of Remote Sensing, 23, 1741-1748.

Huete, A. R., 1988: A soil-adjusted vegetation index (SAVI). Remote Sensing of
Environment, 25, 295-309.

Vogelmann, J. E., S. M. Howard, L. Yang, C. R. Larson, B. K. Wylie, and J. N. Van
Driel,
2001: Completion of the 1990s National Land Cover Data Set for the
conterminous United States. Photogrammetric Engineering and Remote Sensing,
67, 650-662.