‬ TUTORIAL 3 -

IDRISI IMAGE PROCESSING

INTRODUCTORY IMAGE PROCESSING EXERCISES

Image Georegistration Using RESAMPLE

Image Exploration

Image Restoration and Transformation

Image Restoration: Landsat 8

Principal Components Analysis

Supervised Classification

Unsupervised Classification

Change Analysis--Pairwise and Multiple Image Comparison

Data for the exercises in this section are in the \TerrSet Tutorial\Introductory IP folder. The TerrSet Tutorial data can be
downloaded from [Link]

ADVANCED IMAGE PROCESSING EXERCISES

Bayes' Theorem and Maximum Likelihood Classification

Segmentation Classification

Soft Classifiers I: BAYCLASS

Hardeners

Soft Classifiers II: Dempster-Shafer Theory and BELCLASS

Dempster-Shafer and Classification Uncertainty

Vegetation Analysis in Arid Environments

Data for the exercises in this section are in the \TerrSet Tutorial\Advanced IP folder. The TerrSet Tutorial data can be downloaded
from [Link]

TUTORIAL 3 IDRISI IMAGE PROCESSING 191

‬ EXERCISE 3-1
IMAGE GEOREGISTRATION USING
RESAMPLE

Resampling is a procedure for spatially georeferencing an image to its known position on the ground. Often, this procedure is used to
register an image to a universally recognized coordinate reference system such as Lat/Long or Universal Transverse Mercator
(UTM). If the image is already georeferenced but needs to be transformed into another reference system (e.g., from Lat/Long to
UTM), it is advised that you follow the method outlined in Exercise 6-3 on Changing Reference Systems with PROJECT. Resampling
should only be performed when an image is not georeferenced, or when it is not possible to project it. For more information, refer to
the chapter on Georeferencing in the TerrSet Manual.

Even though satellite imagery and other data may often be bought already georeferenced, there are two reasons why you should
consider purchasing non-georeferenced data and doing it yourself. First, you can monitor and reduce the positional error that is
inevitably introduced during any resampling process. A pre-georeferenced image has positional error that is not always documented,
and that may be larger than what you can achieve by doing it yourself.

Second, you can choose the reference system into which the image will be transformed. Resampling is a rubber-sheet transformation
that stretches and warps an image to fit a particular reference system. This process introduces spatial distortion. Some reference
systems, and their associated projections, will introduce more spatial distortion than others for your area. By choosing to do the
resampling yourself, you can choose the reference system that introduces the least amount of spatial distortion. You can also
reference the data to match the reference system of other data you are using.

The resampling procedure may be summarized in three steps as follows:

1. The user identifies the X,Y coordinates of pairs of points that represent the same
place within both the input and output coordinate systems (Figure 1). These are
often referred to as control points or ground control points (GCPs). The
coordinates of the output system may be taken from a map, from another already
georeferenced image, from a vector file, or through surveying either with
traditional instruments or with Global Positioning Systems (GPS).

2. TerrSet derives an equation that describes the relationship between the two coordinate systems.

3. Using this equation, TerrSet converts the input file to the output reference system through what is termed a rubber-sheet
transformation.

In this exercise, we will georeference a raw Landsat Thematic Mapper (TM) image (input reference system) to a previously
resampled Landsat image in a UTM coordinate system (output reference system). TM imagery has a pixel resolution of 30 meters,
and this will be maintained through analysis. The input image, called PAXTON, is within the Paxton quadrangle, just west of Howe Hill

EXERCISE 3-1 IMAGE GEOREGISTRATION USING RESAMPLE 192

in central Massachusetts. We will use a Band 4 TM image from a previous exercise to derive the UTM control points. This image,
called P012R31_5T870916_NN4, is found in the Introductory Image Processing tutorials.

A Open TerrSet Explorer and the Projects tab. Create a new project with the working folder set the Introductory IP tutorial
folder. These folders should be found within the TerrSet Tutorial folder.

B Display the image PAXTON with the Autoscale, Equal Intervals option and the Greyscale palette. This is the infrared band.

Move the cursor across the image and notice that the column positions match the X coordinates (as reported at the
bottom of the screen). From Layer Properties on Composer, choose the Properties tab. Note the values for the minimum
and maximum X and Y and the number of rows and columns. This reference system was entered into the image's
documentation file when it was imported to TerrSet. The reason this particular 'arbitrary' reference system is used will be
explained at the end of this exercise when we consider the positional error introduced during the resampling.

1 When you move the cursor across the image, the row positions and Y coordinates don't match. Why?

C Display the image P012R31_5T870916_NN4 with the Autoscale, Equal Intervals option and the Greyscale palette. This is
also an infrared band.

Move the cursor across the image and notice X and Y coordinates (as reported at the bottom of the screen). From Layer
Properties on Composer, choose the Properties tab. Note the reference system and the values for the minimum and
maximum X and Y coordinates.

The first step in the resampling procedure is to find points that can be easily identified within both the input image and some already
georeferenced map or data layer, i.e., P012R31_5T870916_NN4. The X,Y coordinates of these points in the georeferenced map or data
layer will be the "output" coordinate pairs, while the coordinates from the currently arbitrarily referenced image (PAXTON) will be
the "input" coordinate pairs. Places that make good control points include road and river intersections, dams, airport runways,
prominent buildings, mountain ridges, or any other obvious physical feature.

The input image P012R31_5T870916_NN4 is an entire TM band. Since only a small portion in the upper-left corner corresponds to
the town of Paxton, we will window out the portion we need. This will make displaying and finding good control points easier.

D Run the module WINDOW. Enter as the input filename P012R31_5T870916_NN4 from the Introductory Image Processing
folder. Give an output image name of BAND4UTM. Select Geographical positions as the method to specify the window. For
the coordinates specify:

Minimum X coordinate = 252000

Maximum X coordinate = 264000

Minimum Y coordinate = 4681000

Maximum Y coordinate = 4697000

When WINDOW has completed, display BAND4UTM alongside PAXTON using the Greyscale palette with Autoscale Equal
Intervals. We will use BAND4UTM to determine all the output control points for the rest of this exercise.

Before continuing, close all images and forms on the TerrSet desktop (Ctl-Shift-W).

We are now ready to begin the resample process.

EXERCISE 3-1 IMAGE GEOREGISTRATION USING RESAMPLE 193

E Run the module RESAMPLE. The input file type specifies the type of file to be resampled and can be a raster or a vector file,
or a group of files entered as an RGF. Leave the input file type as raster and specify the input image as PAXTON and the
output image as PAXTONUTM. We will fill in the output reference parameters later.

The input and output reference files to be specified next refers to the set of images to be used to create the GCPs. For the
input reference image select PAXTON and for the output reference image select BAND4UTM. With each selection, the
image will be displayed in a separate window. Although in this case the files specified for the input and reference images
are the same as those specified for the input and output images, the reference images can be any set of images with
corresponding reference systems used in the creation of ground control points.

Before continuing we need to specify the background value, mapping function and the resampling type.

F
Enter 0 as the background value.

A background value is necessary because, after fitting the image to a projection, the actual shape of the data may be angled. In this
case, some value needs to be put in as a background value to fill out the grid. The value 0 is a common choice. This is illustrated in
Figure 2.

The best mapping function to use depends on the amount of warping required to transform the input image into the output
registered image. You should choose the lowest-order function that produces an acceptable result. A minimum number of control
points are required for each of the mapping functions (three for linear, six for quadratic, and 10 for cubic).

G Choose the linear mapping function.

The process of resampling is like laying the output image in its correct orientation on top of the input image. Values are then
estimated for each output cell by looking at the corresponding cells underneath it in the input image. One of two basic logics can be
used for the estimation. In the first, the nearest input cell (based on cell center position) is chosen to determine the value of the
output cell. This is called a nearest neighbor rule. In the second, a distance weighted average of the four nearest input cells is
assigned to the output cell. This technique is called bilinear interpolation. Nearest neighbor resampling should be used when the data
values cannot be changed, for example, with categorical data or qualitative data such as soils types. The bilinear routine is
appropriate for quantitative data such as remotely sensed imagery.

H
Since the data we are resampling is quantitative in character, choose the bilinear resampling type.

EXERCISE 3-1 IMAGE GEOREGISTRATION USING RESAMPLE 194

We are now ready to digitize control points. It is critical to obtain a good distribution of control points. The points should be spread
evenly throughout the image because the equation that describes the overall spatial fit between the two reference systems will be
developed from these points. If the control points are clustered in one area of the image, the equation will only describe the spatial fit
of that small area, and the rest of the image may not be accurately positioned during the transformation to the new reference system.
A rule of thumb is to try to find points around the edge of the image area. If you are ultimately going to use only a portion of an image,
you may want to concentrate all the points in that area and then window out that area during the resampling process.

I
To illustrate how control points are found, zoom into the PAXTON image around the coordinates X 93 and Y 359. This is a
long narrow reservoir in the upper-left portion of the image. Look for a pixel that defines the road intersection going
across the reservoir. This intersection corresponds to the intersection found at X,Y position 253712, 4693988 in
BAND4UTM. Zoom into BAND4UTM at the corresponding location. Notice how difficult it is to determine a precise location
for the intersection in PAXTON because of cell resolution. This is what makes resampling a time-consuming and exacting
task.

We will select a total of 18 well distributed control points throughout the two images. As we develop these points, you can refer to
Figure 3 at the end of this exercise for the approximate location of all the control points. Before you begin to locate and digitize
control points, you may want to adjust the contrast of each image.

J
Zoom back out of both images to the default extent. You can use the Home key when the image is in focus. With PAXTON in
focus, select Layer Properties from Composer. Try adjusting the display maximum down to around 120 and notice that
many features, particularly roads, are more visible. Make a similar adjustment to BAND4UTM. Keep this in mind, that as
you try to discern features in both images, it will be helpful to adjust either of the contrast settings.

K Let’s digitize our first control point. From the RESAMPLE dialog box, notice the Digitize GCP input and output buttons. The
input button refers to the input reference image PAXTON, and the output button refers to the output reference image
BAND4UTM. Click the input button. Notice that a control point is placed in the center of the PAXTON image. We will now
place this point at the location mentioned in Step (i) above, i.e., on the road as it crosses the reservoir at approximately 93
Xand 359 Y. You may want to move the point to the general location and then zoom into the image to place it more
precisely. Notice that as you move the point, the input X and Y values on the RESAMPLE grid change. In addition, you will
notice that as you move the cursor through either the input or output reference images, the area around the cursor will be
magnified and displayed on the right side of the RESAMPLE form.

Once you have placed the GCP on the PAXTON image, click the Output Digitize GCP button. This will place the first GCP in
the BAND4UTM image. Move the first output GCP to the same location on the road as above, in the BAND4UTM image. It
should be placed approximately at X,Y position 253702, 4693981.

We now will place the next three points. We will place one point at a time.

L Zoom back out on both images. We will place GCP 2 at a location below GCP 1 approximately at X, Y position 112.8, 158.3
in PAXTON. Digitize another input GCP, move it to this location. This is at the exit of the reservoir. Select another output
GCP and place it at approximately 252989, 4688317 in the X and Y position in BAND4UTM.

M Next, we will place GCP 3 on the input and output reference images. We will place this GCP on the center of an island in a
reservoir at approximately 213.5 X and 29.5 Y in the input image PAXTON. Locate this reservoir on both images and zoom
into the area. Digitize GCP 3 on both the input and output image and place it at the center of the brightest cell on the island.

2 What was the X,Y coordinate pairs GCP 3 for both the input and reference image?

EXERCISE 3-1 IMAGE GEOREGISTRATION USING RESAMPLE 195

Next, we will place the 4th GCP. As this and subsequent points are digitized you will notice several features, the calculation of the
RMS and residual values and the automatic placement of the corresponding coordinate pair. Each of these features are described
below. First, we will digitize the point.

N On both the input and output image, zoom into the airport in the lower-right corner of each image. We will place GCP 4 at
the intersection of two airstrips at approximately 492.5 X and 65.5 Y on the input image. Place the input coordinate pair at
this location. Notice that once you place the 4th GCP, the 4th output GCP will be interpolated and automatically placed on
the output reference image. Initially, as you add more GCPs, the corresponding interpolated point will be more accurately
placed. Move the output GCP to the correct location at approximately 262981 X and 4683370 Y.

The interpolation is dependent on the mapping function selected. In our case we are using a linear mapping function, so after the
third point, all subsequent points will be interpolated based on the linear polynomial equation, or best fit. As you digitize the
remaining points, the corresponding coordinate pair will be automatically placed, but some adjustment will need to be made
manually.

Also notice that the total root mean square (RMS) and residuals for each control point are now calculated. The residuals express how
far the individual control points deviate from the best fit equation. Again, the best fit equation describes the relationship between the
input image's arbitrary reference system and the output reference system into which it will be resampled. This relationship is
calculated from the control points. A point with a high residual may suggest that the point's coordinates were ill chosen, in either the
input system, the output system, or both.

The total RMS describes the typical positional error of all the control points in relation to the equation. It describes the probability
that a mapped position will vary from the true location. According to US national map accuracy standards, the RMS for images should
be less than 1/2 the resolution of the input image. Recall that TM imagery has a resolution of 30 meters. So in our case, one would
expect that the RMS should be less than 15 meters. The RMS is expressed, however, in input units. Here, we need to understand the
'arbitrary' reference system for PAXTON.

O Open TerrSet Explorer and select the raster file PAXTON, then view its metadata. Notice the properties for this image, in
particular the number of rows and columns and the minimum and maximum X and Y coordinates.

PAXTON's X,Y coordinate system matches the number of rows and columns in the image. This means that one unit in the reference
system is equal to the width of one pixel. In other words, by moving one unit in the X direction, you move one pixel. Therefore, 0.5
units in the reference system is equal to 1/2 the pixel width. The goal, therefore, is to reduce the total RMS error to less than 0.5.

During the resample process and the placement of GCPs, one should constantly be aware of the overall RMS and the residual values.
Notice that some points have higher residuals relative to others. This is not unexpected nor uncommon. As we saw earlier, choosing
control points is not easy. Fortunately, we can choose not to include the bad points and calculate a new equation. Before omitting
points, however, recall a critical issue mentioned earlier: maintaining a good distribution of points. While those points with a very
high residual value tend to be poor, this is not always the case. A few bad points in another part of the image may be "pulling" the
equation and making one good point appear bad. You might choose to remove the most questionable points first. Alternatively, re-
examine the X,Y coordinate positions of your coordinate pairs and reposition them if necessary.

P Let’s digitize the remaining 14 GCPs. Refer to their physical locations in Figure 3 and their precise locations in Table 1. You
can place them by either typing in the coordinates or digitizing each GCP. If you wish to type each GCP, digitize the point
first, then edit the X and Y coordinate.

Remaining GCPs

Point Input X Input Y Output X Output Y

EXERCISE 3-1 IMAGE GEOREGISTRATION USING RESAMPLE 196

 5 417.0 141.5 261354 4685938

6 285.0 200.3 258062 4688399

7 429.0 396.0 263283 4692957

8 216.8 408.2 257460 4694610

9 384.7 501.6 262678 4696160

10 354.4 22.7 258876 4683022

11 186.8 445.8 256850 4695840

12 407.2 309.8 262141 4690689

13 232.9 331.3 257403 4692371

14 152.4 20.6 253231 4684224

15 427.3 458.7 263638 4694714

16 467.8 215 263220 4687664

17 191.6 191.5 255403 4688725

18 297.3 145 258048 4686790

As you enter the GCPs, you should be aware of the total RMS and the residuals for each point. High residual values, for example, over
1.0, are a clue that the coordinate pairs need to be adjusted, or alternate locations found altogether. Remember, our goal is to achieve
an RMS below the input resolution of 0.5.

Q
After completing the placement of all the GCPs, save all the coordinates to a GCP file called PAXTON. Use the Save GCP as
button to save the file. This file can be called up later to add more points or to adjust during the development of your GCPs.

Once you are satisfied with entering and adjusting the GCPs, the final stage is to specify the output reference parameters. These are
the parameters that the resampled image will acquire after the resample process.

R Click on the Output Reference Parameters button. We first need to determine the number of columns and rows for the
output image. These depend upon the extent of the output image, however, so we will first fill in the minimum and
maximum X and Y coordinates.

Enter the following minimum and maximum X and Y's:

Minimum X coordinate = 253000

EXERCISE 3-1 IMAGE GEOREGISTRATION USING RESAMPLE 197

 Maximum X coordinate = 263500

Minimum Y coordinate = 4682000

Maximum Y coordinate = 4695000

This is the bounding rectangle of the output file that will be created. Any bounding rectangle may be requested, and it is quite
common to window out a study area that is smaller than the original image during this process. Note that if the bounding rectangle
extends beyond the limits of the original image, those pixels will be assigned the background value.

We can now calculate the number of columns and rows for the output image. The number of columns for the output file is calculated
from the following equation:

# Columns = (MaxX-MinX)/Resolution

PAXTON is a Landsat Thematic Mapper image which has a resolution of 30 meters. This is the cell resolution that we will want to
retain for the output image. The equation is therefore:

# Columns = (263500-253000)/30 = 350

3 What is the equation for determining the number of rows, and what is the correct number of rows (round the result)?

S Enter 350 columns and the correct number of rows.

T Next select the reference system parameter file UTM-19N from the Georef sub-folder of your TerrSet program folder.
Retain the default meters for reference units and enter 1.0 as the unit distance. Press OK on both dialog boxes.

UTM-19N is the name of the reference system parameter file that corresponds to the Universal Transverse Mercator projection in
Zone 19 (covering Massachusetts). A full discussion of reference system parameter files is found in the chapter on Georeferencing in
the TerrSet Manual.

U
After entering all the required output reference parameters, select OK and again, select OK to run RESAMPLE. The
computer is now performing the last step of the resampling process. The entire image is being transformed into an output
reference system according to the equation calculated from the GCPs.

4 What is the RMS? The overall RMS should be just below the US map accuracy standard.

V
When the resampling is complete, give focus to the output image, PAXTONUTM, then use Layer Properties from Composer
to enable Autoscaling, Equal Intervals and change the palette to Greyscale, if necessary.

W Display the original image, PAXTON, and notice that a clockwise twisting has occurred during the resampling. This spatial
transformation is most evident when looking at the long lakes and the airport runways on the right side of the image.

Georegistering images is an exacting process. Any spatial inaccuracies in the registered images will carry through in all other
analyses derived from the registered data. As with many of the processes we have explored in this Tutorial, the best approach is often
an iterative one, with many rounds of assessment and adjustment.

EXERCISE 3-1 IMAGE GEOREGISTRATION USING RESAMPLE 198

EXERCISE 3-1 IMAGE GEOREGISTRATION USING RESAMPLE 199
‬ EXERCISE 3-2
IMAGE EXPLORATION

With this exercise, we begin an extensive exploration of remotely sensed imagery and image processing techniques. Because
remotely sensed imagery is a common source of data for GIS analysts, and has a raster structure, many raster geographic information
systems provide some image processing capabilities. If you have not already read the Remote Sensing chapter in the TerrSet Manual,
do so now before continuing with this set of exercises.

We will explore different ways to increase the contrast of remotely sensed images to aid visual interpretation, a process known as
image enhancement. We introduced this concept in the display exercises at the beginning of the Tutorial, but we will review and
extend the discussion here because of its importance in image processing and interpretation. We will also learn about the nature of
satellite imagery and the information it carries.

We will use remotely sensed data for the region just west of Worcester, Massachusetts called Howe Hill. Four bands of Landsat
Thematic Mapper (TM) imagery that were acquired by the satellite on September 10, 1987, constitute the data set for this small area.
They are called HOW87TM1, HOW87TM2, HOW87TM3 and HOW87TM4, and correspond to the blue visible, green visible, red visible
and near infrared wavelength bands, respectively.

We begin our investigation of image enhancement by questioning why we need to increase visual contrast in the imagery. In working
with satellite imagery, we will almost always want to use a grey-scale palette for display. This palette choice for auto-display, as well
as other aspects of the display, may be customized in User Preferences.

A Choose User Preferences from the File menu. On the System Settings tab, enable the option to automatically display the
output of analytical modules. Then on the Display Settings tab, set the default quantitative palette to be Greyscale. Choose
to automatically show the title, but not the legend.

B Display the image HOW87TM4 with the Greyscale palette and no autoscaling. Notice that the whole image has a medium
grey color and therefore has very low contrast. The Greyscale palette ranges from black (color 0) to white (color 255), yet
there don't appear to be any white or light grey pixels in the display. To see why this is the case, click Layer Properties on
Composer. Note that the minimum value in HOW87TM4 is 0 and the maximum value is 190. This explains why the image
appears so dark. The brighter colors of the palette (colors 191-255) are not being used.

C
To further explore how the range of data values in the image affects the display, run HISTO from the Display menu. Enter
HOW87TM4 as the input image, choose to produce a graphic output, use a class width of one and the default minimum and
maximum values. When finished, move the histogram to the side to view both the image and the histogram at the same
time.

The horizontal axis of the histogram may be interpreted as if it were the Greyscale palette. A reflectance value of zero is displayed as
black in the image, a reflectance value of 255 is displayed as white, and all values in between are displayed in varying shades of grey.
The vertical axis shows how many pixels in the image have that value and are therefore displayed in that color. Notice also the

EXERCISE 3-2 IMAGE EXPLORATION 200

bimodal structure of the histogram. We will address what causes two peaks in the near infrared band later in the exercise, when we
learn about the information that satellite imagery carries.

As verified by the histogram, none of the pixels in the image have the value of 255. Corresponding to the histogram, there are no
bright white pixels in the image. Notice also that most of the pixels have a value around 90. This value falls in the medium grey range
in the Greyscale palette, which is why the image HOW87TM4 appears predominantly medium grey.

1 If the image HOW87TM4 had a single pixel with reflectance value 0 and one other with the value 255 (all the other data values
remaining as they are) would the contrast of the image display be improved? Why or why not?

Contrast Stretches
To increase the contrast in the image, we will need to stretch the display so that all the colors of the palette, ranging from black to
white, are used. There are several ways to accomplish this in TerrSet, and the most appropriate method will always depend on the
characteristics of the image and the type of visual analysis being performed.

There are two outcomes of stretch operations in TerrSet: changes only to the display (the underlying data values remain unchanged)
and the creation of new image files with altered data values. The former are available through options in the display system, while
the latter are offered through the module STRETCH. There are also two types of contrast stretches available in TerrSet: linear
stretches, with or without saturation, and histogram equalization. All these options will be explored in this section of the exercise.

Simple Linear Stretches

The simplest type of stretch is a linear stretch using the minimum and maximum data values as the stretch endpoints. The term
stretch is quite descriptive of the effect. If the histogram you displayed earlier were printed on a rubber sheet, you could hold the
histogram at the minimum and maximum data values and stretch the histogram to have a wider X axis. With a simple linear stretch,
the endpoints of the data distribution are pulled to the endpoints of the palette and all values in between are re-scaled accordingly.

The easiest way to accomplish a simple linear stretch for display purposes is by autoscaling the image. When autoscaling is used, the
minimum value in the image is displayed with the lowest color in the palette and the maximum is displayed with the highest color in
the palette. 1 All of the values in between are distributed through the remaining palette colors.

D
With the HOW87TM4 display in focus, choose Layer Properties on Composer. For the Autoscaling options, click between
Equal Intervals (on) and None (off) a few times, closely examining the overall change in contrast as well as the effects in
the darkest and lightest areas of the image. Notice that the contrast increases with Equal Intervals on.

2 Draw a rough sketch of the histogram for HOW87TM4 with autoscaling. Label the X axis with palette indices 0-255 rather than
data values. On that axis, note where the minimum and maximum data values lie and also mark where the palette colors black,
white, and medium grey lie.

1 Autoscaling uses the Display Min and Display Max values from the image documentation file and matches those to the autoscaling minimum and maximum
values in the palette file. We will return to this later. For now, assume that the minimum and maximum data values are equal to the minimum and maximum
display values for the image and the autoscaling minimum and maximum values are 0 and 255 in the palette file.

EXERCISE 3-2 IMAGE EXPLORATION 201

Note that autoscaling does not change the data values stored in the file; it only changes the range of colors that are displayed.
Although autoscaling often improves contrast, this is not always the case.

E Display HOW87TM1 with the Greyscale palette. Again, open Layer Properties from Composer and click autoscaling on and
off. Notice how little contrast there is in either case. Then, also in Composer, move to the Properties tab and click the
Histogram button on the Layer Properties dialog box. (The module HISTO is called and uses the data values from the file,
and is therefore not affected by any display contrast enhancements, such as autoscaling, that are in effect in the display.)

3 What are the min and max values in the image? What do you notice about the shape of the histogram? How does this explain
why autoscaling does not improve the contrast very much?

Autoscaling alters the display of an image. If it is desirable to create a new image with stretched data values, then the module
STRETCH is used. To achieve a simple linear stretch with STRETCH, choose the linear option and accept the default to use the
minimum and maximum data values as the endpoints for stretching. The stretched image, when displayed, will be identical to the
autoscaled display. (You may try this with one of the images if you wish.)

Linear Stretches with Saturation

We can achieve better contrast by applying a linear stretch with saturation to the image. When we use saturation with a stretch, we
set new minimum and maximum display values that are within the original data value range (i.e., the minimum display value is
greater than the minimum data value and the maximum display value is less than the maximum data value). When we do this, all the
values that lie above the new display maximum are assigned to the same last palette color (e.g., white) and all those below the new
display minimum are assigned to the same first palette color (e.g., black). We therefore lose the ability to visually differentiate
between those "end" values. However, since most remotely sensed images have distributions with narrow tails on one or both ends,
this loss of information is only for a small number of pixels. Most pixels may then stretch across more palette colors, yielding higher
visual contrast and enhancing our ability to perform visual analysis with the image.

The data values that are assigned the lowest and highest palette colors are called the saturation points. There are two ways to
produce a linear stretch with saturation in TerrSet. You may set the saturation points interactively through Composer/Layer
Properties, or you may use the STRETCH module. The former affects the display only, while the latter produces a new image that
contains stretched values. We will experiment with both methods.

F Bring the HOW87TM1 display window into focus (or re-display it if it is closed). Choose Layer Properties in Composer.
The Contrast Settings area of the dialog box is active only when autoscaling is turned on, so turn it on. The default setting
corresponds to a simple linear stretch, with the minimum and maximum data values as the endpoints (11 and 215). Since
the histogram showed a very long thin tail at the upper end of the distribution, it is likely that lowering the Display Max
value will have the greatest effect on contrast. Slide the Display Max down by clicking to the left of the marker. Each time
you click, note the change in the display and the new saturation point value shown in the box to the right of the slider.

G
Click the Revert button to go back to the original autoscaled settings. Now move the Display Min marker up incrementally.

4 Why does contrast actually become worse as you increase the amount of saturation on the lower end of the distribution?
(Hint: recall the image histogram.)

EXERCISE 3-2 IMAGE EXPLORATION 202

Saturation points for display are stored in the image documentation file's Display Min and Display Max fields. By default, these are
equal to the minimum and maximum data values. These may be changed by choosing Save Changes and OK in the Layer Properties
dialog. They may also be changed through the Metadata utility in TerrSet Explorer. Altering these display values does not affect the
underlying data values, and therefore will not affect any analysis performed on the image. However, the new Display Min and Max
values will be used by Display when autoscaling is in effect.

Now we will turn to the linear stretch with saturation options offered through the module STRETCH. A linear stretch with saturation
endpoints may be created with the linear stretch option, setting the lower and upper bounds for the stretch to be the desired
saturation points. This works in the same way as setting saturation points in Layer Properties. The difference is that with STRETCH, a
new image with altered values is produced.

STRETCH also offers the option to saturate a user-specified percentage (e.g., 5%) of the pixels at each end (tail) of the distribution. To
do so, choose the linear with saturation option and give the percentage to be saturated.

H Run STRETCH with HOW87TM4 to create a new file called TM4SAT5. Choose the linear with saturation option and give 5
as the percentage to be saturated on each end. Do the same with HOW87TM1, calling the output image TM1SAT5.
Compare the stretched images to the originals.

The amount of saturation required to produce an image with "good" contrast varies and may require some trial and error
adjustment. Generally, 2.5-5% works well.

Histogram Equalization
The histogram equalization stretch is only available through the STRETCH module and not through the display system. It attempts to
assign the same number of pixels to each data level in the output image, with the restriction that pixels originally in the same
category may not be divided into more than one category in the output image. Ideally, this type of stretch would produce a flat
histogram and an image with very high contrast.

I Try the histogram equalization option of STRETCH with HOW87TM4. Call the output stretched image TM4HE. Compare
the result with the original, then display a histogram of TM4HE.

The histogram is not exactly flat because of the restriction that pixels with the same original data value cannot be assigned to
different stretch values. Note, though, that the higher the frequency for a stretched value, the more distant the next stretched value is.

J Use HISTO again with TM4HE, but this time give a class width of 20. In this display, the equalization (i.e., flattening) of the
histogram is more apparent.

According to Information Theory, the histogram equalization image should carry more information than any other image we have
produced since it contains the greatest variation for any given number of classes. We will see later in this exercise, however, that
information is not the same as meaning.

Exploring Reflectance Values

We will now move on to explore what these remotely sensed images "mean." To facilitate this exploration, we will first create a raster
group file of the original images and one of the enhanced images created earlier. This will allow us to link the zoom and window
actions as well as use the Identify tool across all the images belonging to the group.

EXERCISE 3-2 IMAGE EXPLORATION 203

K Close any display windows that may be open.

L Create a raster group file in TerrSet Explorer. 2 From the Files pane, select the files HOW87TM1, HOW87TM2, HOW87TM3,
HOW87TM4 and TM4SAT5. Then right-click and select Create Raster Group file. By default, a file named RASTER
[Link] is created. Select this file, right-click and rename it to HOW87TM.

M
Open DISPLAY Launcher and activate the Pick List. Note that the group file, HOW87TM, now appears in the list of raster
files in the Working Folder and that there is a plus sign next to it. This indicates that it is a group file. Clicking on the plus
sign expands the Pick List to show all the members of the group. If you wish to use any of the group display features, group
members must be displayed from within the group file and with their full "dot-logic" names. The easiest way to do this is
to invoke the Pick List, expand the group file, then choose the file from the list of group file members. Choose TM4SAT5
from the list. Note that the name in the DISPLAY Launcher file input box reads HOW87TM.TM4SAT5. This is the full "dot
logic" name that identifies the image and its group. Choose the Grey Scale palette and display the image. (Alternatively, you
can display members of a group with the dot-logic from TerrSet Explorer.)

N Also display the four original images, HOW87TM1 through HOW87TM4, in the same manner with the Grey Scale palette.
Do not apply autoscaling or change the contrast for any of these images. We want to be able to visually compare the actual
data values in these original bands. Arrange the images next to each other on the screen so that you can see all five at once.
If you need to make them smaller so they can all be seen, follow this procedure:

Position the cursor over the lower right edge of each map window until the cursor becomes a double arrow, then drag the
map window to the desired size. If necessary, you can always return to the original display size by pressing the End key.

Because the contrast is low in all of the original images, we will use the stretched image, TM4SAT5, to locate specific areas to query.
However, it is the data values of the original files in which we are interested.

There are three land-cover types that are easily discernible in the image: urban, forest and water. We want to now explore how these
different cover types reflect each of the electromagnetic wavelengths recorded in the four original bands of imagery.

O Draw three graphs as in Figure 1 and label them water, forest and urban.

Figure 1

2 Note that all the files of a group must be stored in the same folder. If you are working in a laboratory situation, with input data in a Resource Folder and your
output data in the Working Folder, you will need to copy the input files HOW87TM1-HOW87TM4 into your Working Folder, where TM4SAT5 is stored, before
continuing with the exercise.

EXERCISE 3-2 IMAGE EXPLORATION 204

To examine reflectance values in all four images we will use the Identify tool feature that allows simultaneous query of all the images
within the same map window.

P This time, from TerrSet Explorer, select all five images, H0W87TM1, H0W87TM2, H0W87TM3, H0W87TM4, and
TM4SAT5. When all five images are selected, right-click and select Add Layer. Then click the Identify icon on the toolbar.
(Note that by default Identify mode is activated.) A small Identify box opens to the right of the map window. Find three to
four representative pixels in each cover type and click on the pixels to check their values. The reflectance values of the
queried pixel in all five images in the map window appear in the table. Determine the reflectance value for water, forest
and urban pixels in each of the four original bands. Fill in the graphs you drew in step o) for each of the cover types by
plotting the pixel values.

5 What is the basic nature of the graph for each cover type? (In other words, for each cover type, which bands tend to have high
values and which bands tend to have low values?)

You have just drawn what are termed spectral response patterns for the three cover types. With these graphs, you can see that
different cover types reflect different amounts of energy in the various wavelengths. In the next exercises, we will classify satellite
imagery into land cover categories since land cover types have unique spectral response patterns. This is the key to developing land
cover maps from remotely sensed imagery.

We will now return to two outstanding issues that were mentioned earlier but not yet resolved. First, let's reconsider the shape of the
histogram of HOW87TM4. Recall its bimodal structure.

6 Now that you have seen how different image bands (or electromagnetic wavelengths) interact with different land cover types,
what do you think is the land cover type that is causing that small peak of pixels with low values in the near infrared band?

"Information" versus "Meaning"

Now, let us return briefly to our stretched images and reconsider how stretching images may increase contrast and therefore
"information," but not actually add any "meaning."

Q Use STRETCH with HOW87TM1, choosing a histogram equalization and 256 levels. Call the output TM1HE. Then also
display TM1SAT5.

Note how different these images are. The histogram equalized version of Band 1 certainly has a lot of variation, but we lose the sense
that most of the cover in this image (forest) absorbs energy in this band heavily (because of moisture within the leaf as well as plant
pigments). It is best to avoid the histogram equalization technique whenever you are trying to get a sense of the
reflectance/absorption characteristics of the land covers. In fact, in most instances, a linear with saturation stretch is best. Remember
also that stretched images are for display only. Because the underlying data values have been altered, they are not reliable for
analysis. Use only raw data for analysis unless you have a clear reason for using stretched data.

EXERCISE 3-2 IMAGE EXPLORATION 205

Creating Color Composites
In the final section of this exercise, we will explore the creation of color composite images as a type of image enhancement. Up to this
point in the exercise, we have been displaying single bands of satellite imagery. Color composite images allow us to view the
reflectance information from three separate bands in a single image.

In TerrSet, the 24-bit color composite image is used for display and visual analysis. It contains millions of colors and the contrast of
each of the three bands can be manipulated interactively and independently in Composer on the display system.

We will now create a 24-bit natural color composite image using the three visible bands of the same imagery for Howe Hill as we
examined above. 3

R
Run COMPOSITE from the Display menu. Specify HOW87TM1 as the blue image band, HOW87TM2 as the green image
band and HOW87TM3 as the red image band. Give COMPOSITE123 as the output filename. Choose a linear with saturation
points stretch. Choose to create a 24-bit composite with the original values. Do not omit zeros and saturate 1%.

The resulting composite image retains the original data values but display saturation points are set such that 1% on each end of the
distribution of each band is saturated. These can be further manipulated from the Layer Properties dialog box. However, for now,
leave these as they are.

S Use the Identify tool to examine some of the values in the composite image. Note that the values of the red, green, and blue
bands are all displayed. Try to interpret the values as spectral response patterns across the three visible bands.

7 Look back at the spectral response patterns you drew above for water, forest and urban cover types. Given the bands we have
used in the composite image, describe why each of these cover types has its particular color in the composite image.

Compositing is a very useful form of image enhancement, as it allows us to see simultaneously the information from three separate
bands of imagery. Any combination of bands may be used, and the choice of bands often depends upon the application. In this
example we have created a natural color composite in which blue reflectance information is displayed with blue light in the computer
display, green information with green light and red information with red light. Our interpretation of the spectral response patterns
underlying the particular colors we see in the composite is therefore quite intuitive—what appears as green in the display is
reflecting relatively high on the green band in reality. However, it is very common to make color composite images from other bands
as well, some of which may not be visible to the human eye. In these cases, it is essential to keep in mind which band of information
has been assigned to which color in the composite image. With practice, the interpretation of composite images becomes much
easier. 4

T Create a new composite image using the same procedure as before, except give HOW87TM2 as the blue band, HOW87TM3
as the green band, HOW87TM4 as the red band and FALSECOLOR as the output image name.

This type of composite image is termed a false color composite, since what we are seeing in blue, green and red light is information
that is not from the blue, green and red visible bands, respectively.

3 See Exercise 3-1 on Composites for creating 24-bit RGB composites on the fly from Composer.
4 For practice in interpreting colors as mixes of red, green and blue light, open Symbol Workshop from the toolbar. Choose one palette color index and vary the
amount of red, green, and blue, observing the resulting colors. Experienced image analysts can estimate the relative reflectance values of the three input
images just by looking at the colors in the composite image.

EXERCISE 3-2 IMAGE EXPLORATION 206

8 Why does vegetation appear in bright red colors in this image?

Satellite imagery is an important input to many analyses. It can provide timely as well as historical information that may be
impossible to obtain in any other way. Because the inherent structure of satellite imagery is the same as that of raster GIS layers, the
combination of the two is quite common. The remainder of the exercises in this section illustrate the use of satellite imagery for land
cover classification.

EXERCISE 3-2 IMAGE EXPLORATION 207

‬ EXERCISE 3-3
IMAGE RESTORATION AND
TRANSFORMATION

In this exercise, we will explore the use of several techniques for image restoration. Restoration techniques are preprocessing
techniques for the removal of noise or flaws in imagery due to either sensor detection errors or natural noise from atmospheric
effects. TerrSet provides a range of techniques to address these issues. The modules DESTRIPE, PCA, and ATMOSC will be used here
to explore radiometric correction and noise removal in imagery. With DESTRIPE and PCA we will explore the removal of noise due
primarily to sensor errors. These errors are common since satellites transfer and receive vast amounts of digital data from many
miles above the earth. We will also explore the removal of noise caused by the scattering of solar radiation, which can result in haze.
Given the components of the atmosphere, reflectances can be affected by the interaction between incoming and outgoing
electromagnetic radiation, which alters the true ground-leaving radiance. The module ATMOSC attempts to account for these effects
by removing or dampening the resulting haze.

Removing Sensor Error using DESTRIPE

In the first part of this exercise, we will attempt to address image noise due to sensor error, often occurring in the form of striping or
banding. This is very typical with older imagery but can occur with any sensor platform. Striping or banding is systematic noise in an
image that results from variation in the response of the individual detectors used for a particular band. This usually happens when a
detector goes out of adjustment and produces readings that are consistently much higher or lower than the other detectors for the
same band.

The procedure that corrects systematically bad scan lines in an image is called destriping. It involves the calculation of the mean (or
median) and standard deviation for the entire image and then for each detector separately. It works on both horizontal and vertical
scan lines. Examples of a horizontal scan line detector include MSS and TM, while SPOT is an example of a vertical scan line detector.

A With the image NJOLO2 displayed, use the Add Layer option in Composer to add the other two raster bands, NJOLO1 and
NJOLO3, to the same map composition. Then in Composer, highlight NJOLO1 and select the blue icon on Composer to
assign it the blue component. Then use the cursor to highlight NJOLO2 and select the green icon to assign it the green
component. Finally, highlight NJOLO3 and select the red icon to assign it the red component. Once the red band is selected
and assigned the red component, you will see the false color display in the map window.

EXERCISE 3-3 IMAGE RESTORATION AND TRANSFORMATION 208

The false color composite highlights the severity of the detector error. Given that the striping is perfectly vertical, we can easily
reduce this error using the module DESTRIPE. This module works on perfectly horizontal or vertical noise by calculating a mean and
standard deviation for the entire image and then for each detector separately. Then the output from each detector is scaled to match
the mean and standard deviation of the entire image. Details of the calculation are given when the module has finished running.

B Open the module DESTRIPE. Enter NJOLO2 as the input image. Call the output image NJOLO2D. Set the number of
detectors equal to the number of columns (509) and select Vertical orientation for the striping. Then hit OK to run the
module.

C
With the result, replace NJOLO2 in Step A above with NJOLO2D and display the new false color composite.

The new composite should show a remarkably less noisy image. Since only band 2 was noisy, all the bands are now ready for
analysis. In the next section, we will look at removing noise due to a combination of factors.

Removing Sensor Error and Haze Removal

with PCA
In this section of the exercise, we will explore using Principal Components Analysis (PCA) for removing noise in imagery that has
already been geocorrected.

D Using DISPLAY Launcher, display the Map Composition VIETNAM. By default, only band 1, VIET1, is displayed. In
Composer, however, you will notice that all the bands are present in the Map Composition. You can display each band,
moving from band 2 to band 7, by selecting the check box just to the left of the filename. Select each of the bands to view
each band’s level of noise.

EXERCISE 3-3 IMAGE RESTORATION AND TRANSFORMATION 209

As each band is displayed, you will see a reduction in the level of noise, although all bands are affected to some degree. The other
striking feature is that, although the noise is striped as in the previous section, it is neither horizontal nor vertical. If satellite data are
acquired from a distributor already fully georeferenced, then radiometric correction via DESTRIPE is no longer possible. This
Landsat TM image from the coast of Vietnam was already geocorrected when it was received. In this case, Principal Components
Analysis can be used on the group of input bands.

Running PCA transforms a group of bands into statistically separate components. The last few components usually represent less
than 1 percent of the total information available and tend to hold information relevant to noise, and in our case, striping. If these
components are removed completely and the rest of the components are reassembled, the improvement can be dramatic. The
striping effect may even disappear.

E Open the module PCA. Specify Forward T-Mode as the analysis type and the covariance matrix unstandardized option.
Insert the layer group file VIETNAM. Specify 7 for the number of components to be extracted. Give an output prefix of PCA
and select the option to output the complete text. Then click OK.

When PCA finishes, it will produce a set of images with the prefix PCA and output a table of statistics from the transformation.

F Display each of the seven component images, either within a single window or independently.

Once the images are displayed, notice how each subsequent image contains more and more noise. Also, according to the table of
statistics produced from PCA, notice that Component 1 explains 93% of the total variability across all the bands (read from the %var.
line under each component)

1 What is the total percent variance explained by the last four bands? By the last three bands?

The Results table from the PCA module shows the statistics from the transformation, including the variance/covariance matrix, the
correlation matrix, and the component eigenvectors and loadings. Analyzing the components section of the table, the rows are
ordered according to the band number and the column eigenvectors, reading from left to right, represent the set of transformation
coefficients required to linearly transform the bands to produce the components. Similarly, each row represents the coefficients of
the reverse transformation from the components back to the original bands. Multiplying each component image by its corresponding
eigenvector element for a particular band and summing the weighted components together reproduces the original band of
information. If the noise components are simply dropped from the equation, it is possible to compute the new bands, free of these
effects. This can be achieved manually using Image Calculator in TerrSet. But an easier method is just to use the inverse PCA option in
the PCA module.

G In PCA, select to perform an inverse t-mode PCA. Specify PCA_T-MODE_COMPS as the RGF component filename and PCA_T-
MODE as the eigen filename. Then enter 1-2 for list of components to be used, NEWBAND as the prefix for the output files,
and 1-7 for the output bands to be created. Then click OK.

Once the operation is completed, display and compare the original band 1, VIET1, to the transformed PCA NEWBAND1. Notice the
significant improvement in the display. Doing the reverse transformation with only the first two components significantly reduced
the noise. These two components also contain 97.3% of the overall variance in the original band. We can add more components to
capture more of the original variability, but we will need to weigh this against increasing noise.

H Run PCA again. This time only use the eigenvector for component 1 to create a new band 1. Call the new output
NEWBAND1_.

EXERCISE 3-3 IMAGE RESTORATION AND TRANSFORMATION 210

2 If you were to use Image Calculator to calculate the new band, what equation would you use to create NEWBAND1_1 above?
What was the total variance explained in the result compared to the original image?

You can experiment with entering any number of components and their respective eigenvectors. You would not want to do a reverse
transformation on all the bands. If you recall, only bands 1, 2, 3, and 6 seem to contain radiometric noise. The other bands should be
left as is for further analysis. Also, once you have the reverse bands the way you would like, they will need to be stretched to a byte
level from 0 to 255 for further use with the other original bands. You can use the modules STRETCH or FUZZY.

Atmospheric Correction to Remove Haze with

ATMOSC
In the previous sections of this exercise, we demonstrated the removal of systematic noise due to bad scan data. The modules
DESTRIPE and PCA substantially reduced banding and striping noise due to sensor errors. In this section, we will explore the removal
of radiometric errors due to haze and demonstrate atmospheric correction using the module ATMOSC.

The images we will use to demonstrate atmospheric correction are taken from a Landsat 5 TM image of Southeastern New England,
USA, including Boston, Worcester and Cape Cod, Massachusetts, and Providence, Rhode Island. The date of the image is September
16, 1987. The goal is to reduce or remove any atmospheric influence by eliminating haze or other interferences. First, we will remove
haze using the Cos(t) model, and then we will verify our results using “pure” spectral libraries.

ATMOSC calls for several inputs, particularly for the full model, which requires the calculation of Optical Thickness. Usually, most of
the data input required for the module can be found or calculated from the accompanying metadata. You can also consult with basic
Image Processing texts for some of the required parameters.

ATMOSC also needs the meteorological conditions for that day. For the vicinity of Boston, Worcester, and Providence, we contacted
the local weather bureau for Worcester, Massachusetts and were provided with the following weather information for that day:

Sept 16, 1987, Worcester Regional Airport (KORH)

10:00 DST Temperature 67 F Dew Point 51 F Visibility 30 mi Station Pressure 28.95

11:00 DST Temperature 70 F Dew Point 53 F Visibility 30 mi Station Pressure 28.92 SLP 30.01

I Display the image P012R31_5T870916_NN3, using the Greyscale palette and autoscaling (Equal Intervals). The last
character in the band filenames is the band number. You are now displaying band 3. You will notice banding especially in
the ocean areas east of the center of the image. This is Boston Harbor.

J We will next create a false color composite. With band 3 in focus, add two more raster layers to this image. Either hit the
‘R’ key or click Add Layer in Composer. Add bands P012R31_5T870916_NN2 and P012R31_5T870916_NN4 to the layer
P012R31_5T870916_NN3.

Once all three layers are present within the same map window, you can use the features in Composer to assign each band to
represent blue, green and red in a combined window display.

EXERCISE 3-3 IMAGE RESTORATION AND TRANSFORMATION 211

K With the map layer containing all three bands in focus, move the cursor over to Composer. Highlight band 2 and then
select the blue icon on Composer to assign it the blue component. Then use the cursor to highlight band 3 and select the
green icon. Likewise, highlight band 4 and select the red icon to assign it the red component. Once the red band is selected
and assigned the red component, you will see the false color display in the map window.

Once the composite image is displayed, the distortion caused by haze is evident. The composite image, especially when the visible
bands are used, best illustrates the distortions caused by energy scattering in the atmosphere as well as noise due to sensor
problems. Look particularly at water bodies in the interior of the image, as well as along the coast. Explore the image, looking at the
ocean, lakes, urban areas and vegetation areas.

To correct for these errors, the user must collect the required metadata for the imagery used. The data used in this exercise was
originally downloaded from the University of Maryland with its accompanying metadata. Let’s examine this file.

L Open the file [Link] using Edit and examine the data. We will use the time and date information, the sun
elevation, the satellite name, and for each band, the wavelength, gain and bias. You may find that printing this four-page
file will be useful over the course of this exercise.

At this point, we have everything needed to run the ATMOSC’s cos(t) model.

M Open the module ATMOSC and select the cos(t) model. Enter the input image as P012R31_5T870916_NN4. By entering the
input image first, the module will read the minimum and maximum values from its documentation file and enter in default
Dn min and Dn max. These can be edited later.

N Next, we need to enter the year, month, date and GMT (Greenwich Mean Time). If not already open, open the metadata file:
[Link]. Locate the line that begins: “Start_Date_Time.” This line will list the year, month, day, and time (1987, 09,
16, 14:53:59.660, respectively). However, the ATMOSC module requires that the time be in decimals. Round the minutes,
seconds and milliseconds to the nearest minute (54) and divide by 60 to get the minutes to a single decimal place, i.e.,
14.90.

O Next, we need to determine the wavelength of the band center. We will again use the metadata file. Each band has its own
section in the metadata file. Locate the section for band 4. The second line for band 4 should read: “File_Description=Band
4.” Below this line we find the wavelength information in the line: “Wavelengths.” The values here, 0.76 and 0.90, are the
minimum and maximum wavelengths in microns for this band. Average these values to find the wavelength of the band
center (.83) and enter it into ATMOSC.

The next input is the DN haze value, which refers to the Digital Number or value that must be subtracted to account for visible haze.
This can be determined by isolating extremely low reflectance values in the images such as deep lakes or fresh burn scars.

P To assess the DN haze value, once again display the false color composite that you created in Step (j), and find a large deep
lake. These are areas that should have very low reflectance at all wavelengths. The Wachusett Reservoir, just north of
Worcester in the upper-left quadrant (column 2840, row 1375) of the image is a good location for this. Zoom into this
area, and using the Identify tool, find the lowest red (band 4) value in the lake. Make sure you have band 4 selected in
Composer. Enter this value (it should be about 6 for this band) into ATMOSC for the DN haze. Remember that the image is
bordered by background values that are 0.

For the next set of inputs, we need to calibrate the radiance. This is done using the gain and bias values in the band 4 section of the
metadata file. Please note that the metadata uses the term “bias” while ATMOSC uses the term “offset.” This exercise will use the
ATMOSC terminology.

Q Select the offset/gain radiance calibration option. Then, reading from the metadata file for the fourth band from the gains
and biases line, the gain and offset given are 8.14549 and -1.51, respectively. The module requires that these inputs be in

EXERCISE 3-3 IMAGE RESTORATION AND TRANSFORMATION 212

 mWcm-2sr-1um-1. To test the units, multiply the gain by the highest possible image value (255 for byte images) and add
the offset. If the result is between 10 and 30, then the units are correct. If the result is too large by a factor of 10, then the
units are Wm-2sr-1um-1 and the decimal for both offset and gain must be shifted one place to the left. Enter an offset of -
.0151 and a gain of .0814549. For more details, see the Notes section of the ATMOSC Help file.

R
The next input is the satellite viewing angle, which is 0 for all Landsat Satellites. This is the default setting in ATMOSC. For
other satellite platforms, the user must check to determine the viewing angle for the scene, although it is usually zero.

S Finally, the sun elevation must be entered. Near the beginning of the metadata file, look for the line Solar_Elevation. Enter
the sun elevation, 45.18. Give the output image name BAND4COST and click on OK.

T
Repeat this process for bands 2 and 3.

3 What were the values used to correct bands 2 and 3?

Next we need to create a composite of the corrected images.

U
Repeat the steps in (j) and (k) above using the transformed, atmospherically corrected images to create a false color
composite. Compare this composite with the one made from the pre-transformed bands. Explore the differences,
particularly in the shallow coastal regions, urban areas, and Wachusett Reservoir.

You will notice that much of the haze has been removed. This haze is most likely a result of attenuation due to particles, both
moisture and solid materials, in the atmosphere. If you look closely, however, you will notice that other noise is present that is not
due to atmospheric effects, but due to possible errors with the sensors on board the satellite. The image could be corrected further
through PCA.

Evaluation of ATMOSC
Researchers at the USGS Spectroscopy Lab have measured the spectral reflectance of hundreds of materials in a laboratory setting.
The result is a compilation of a spectral library for each material measured. For each material, a “pure” signature of spectral
reflectance is produced. It is pure because, done in a lab setting, the measurement of the spectral response is void of any atmospheric
effects and other attenuation effects. The library can be used as a reference for material identification in remotely sensed imagery,
particularly hyperspectral imagery. After running ATMOSC, the values in the output images are reflectances, the same scale of values
found in the spectral library. Although used to calibrate remote sensors, they can also be used to validate our results. One of the
materials measured by the USGS is that of “lawn grass.” As a pure spectral signature, lawn grass elicits the following spectral
response pattern across six bands of TM:

Band Spectral Reflectance

1 4.043227E-02

EXERCISE 3-3 IMAGE RESTORATION AND TRANSFORMATION 213

 2 7.830066E-02

3 4.706150E-02

4 6.998996E-01

5 3.204015E-01

7 1.464245E-01

Table 1: Spectral reflectance values for lawn grass as reported from the USGS spectral library for TM. (More detail on spectral
libraries can be found at [Link]

We have digitized a test area for the TM images used in this exercise from a golf course. This will approximate a large contiguous
“lawn grass” area needed to validate the atmospheric correction on each of the bands. The golf course is at approximately column
2405 and row 1818. A raster file with the name LAWN GRASS exists that can be used to overlay on the images to verify its location.

Using the image LAWN GRASS, we will extract out the average values from the three bands created from ATMOSC.

V
Run the module EXTRACT. Specify the feature definition image as LAWN GRASS and the image to be processed as
BAND2COST. Select average as the summary type and tabular output type. Run EXTRACT again on BAND3COST and
BAND4COST.

4 What were the reflectance values extracted for each of the three corrected bands? The raster file LAWN GRASS is a Boolean
image with values of 1 for the areas of interest (lawn grass) and zero for the background.

The results should indicate very similar reflectance values for our three corrected bands, TM bands 2, 3, and 4.

EXERCISE 3-3 IMAGE RESTORATION AND TRANSFORMATION 214

‬ EXERCISE 3-4
IMAGE RESTORATION: LANDSAT 8-9

This tutorial demonstrates the Landsat module used for the import and restoration of any of the Landsat satellite archives, including
MSS, TM, ETM+ and Landsat 8-9 OLI and TIRS. The Landsat module will not only import raw DN, but also undertake atmospheric
correction on the multispectral data to produce reflectance or radiance imagery, as well as calibrate the thermal bands to
temperature data.

The LANDSAT module utilizes the MTL metadata text files distributed with each Landsat scene in the USGS archive. This revised
archive provides data in GeoTIFF format with an accompanying MTL metadata text file. One site that provides the simplest method
for searching and downloading these data is the USGS Earth Explorer portal, at [Link] The USGS Glovis site
is another portal that provides the Landsat data in the required format.

The Landsat satellite archive is one of the most recognizable remote sensing platforms with imagery spanning four decades
providing the longest time series of imagery data available to earth scientists. In this tutorial, we will import an example of Landsat 8
OLI/TIRS imagery from coastal British Columbia, Canada, and perform some preprocessing steps to improve the quality and contrast
of the imagery. If you already have Landsat 8 imagery in the required format, you can certainly use this in place of the British
Columbia data.

A The first step is to download the required scene. Open your browser to the USGS Earth Explorer portal at
[Link]. You will need to register and login before you can download data. If you have not already done so,
please register and then log in to the Earth Explorer site before continuing.

B Once you have logged in to the Earth Explorer site, search for the scene with path and row of 50/24 and an acquisition
date of August 3, 2014. On the Earth Explorer site, select the data set to use, the Landsat archive L8-9 OLI/TIRS. Having
selected the dataset, you can then search using the scene ID above.

C Once you have located this scene, locate the option to view the scene in Earth Explorer.

The interactive map in Earth Explorer shows a pre-generated composite for this Landsat scene. We can see that this image is mostly
cloud free, except for some fog in the lower portion of the image. Much of the mountainous areas of this image are covered by
glaciers, while the lower elevation areas are covered by temperate rainforest. Landsat 8-9 imagery can be used effectively to map
both surface types. In Earth Explorer, your interactive window may look like the image below.

EXERCISE 3-4 IMAGE RESTORATION: LANDSAT 8-9 215

D
Once you have successfully located the scene, proceed to download. You will be presented with several download options.
Select the Level 1 GeoTIFF data product. Save it to your computer and store it in a folder you will set as your working
folder in TerrSet.

E Once downloaded you must unzip the archive file so that you have an individual GeoTIFF file for each band and the
accompanying MTL text file in the same folder.

F After the data is unzipped and ready for import, open TerrSet and set your working folder to the folder containing your
unzipped Landsat 8 scene.

G Open the Landsat data archive module from the Import menu. Then load the Landsat metadata file:
LC08_L1TP_050024_20140803_20200911_02_T1.[Link]. Then click OK.

Opening this MTL will automatically enter all band information for this scene. Landsat 8 imagery contains 4 visible, 5 NIR/SWIR, and
2 thermal bands, along with a quality check band. The top of the form also displays basic data about the scene, including the sensor,
image date, and path and row.

EXERCISE 3-4 IMAGE RESTORATION: LANDSAT 8-9 216

To investigate the imagery, we will first import 3 bands to create a false color composite.

H Under the include column, click on all bands so that only 4, 5, and 6 are set to “yes”. In the output image name column, add
the suffix “_raw” to these 3 bands (e.g., LC08_L1TP_050024_20140803_20200911_02_T1_B4_RAW). Leave all other
defaults and click OK.

Once the module completes running, the first image will autodisplay. Landsat 8 data improves upon previous versions by having a
16-bit radiometric resolution with raw DN values ranging from 0 and 65,535. Next, we will create a quick composite to better view
the data.

I With the raw band 4 image displayed, use the add layer button on Composer to add bands 5 and 6 to the map composition.
Again, from Composer, use the composite red, green, and blue icons for the three bands in Composer. In Composer, select
band 4 and then click the blue icon. Select band 5 and click the green icon. Then select band 6 and click the red icon. You
may have to extend the size of Composer to adequately view the band names.

EXERCISE 3-4 IMAGE RESTORATION: LANDSAT 8-9 217

The result will look very similar to the quick browse image in Earth Explorer. This is a false color composite using the red, NIR and
SWIR bands. Forest appears as green, and ice and snow as blue, though brighter than reality.

Next, we will import the imagery and convert it to reflectances.

J In the LANDSAT module, select the “Convert to reflectance” option under Multispectral bands and set the reflectance
correction to None. Change the output image names for bands 4, 5, and 6 by adding the suffix “_reflecNone”. Click OK.

Once the module runs, investigate the outputs. Note that while the output appears visually similar, the values have changed, now
roughly between 0 and 1. Unlike raw DN, reflectance values are a physical property of the surface, where values near 0 represent
surfaces that are very absorptive at a particular wavelength and those near 1 very reflective. Note that the highest values in Band 4
(red visible) are the snow and ice surfaces, while everything else is mostly absorptive.

Next, we will correct imagery for atmospheric haze. Atmospheric correction is typically performed by the ATMOSC module in
TerrSet. For Landsat versions prior to Landsat 8, the Landsat module will make a call to the ATMOSC module. For Landsat 8, the
Landsat module will use the metadata contained in the MTL text file. See the tutorial on atmospheric correction using ATMOSC for
more detail on this restoration technique.

K In the Landsat module, select the Dark-object subtraction radio button under Reflectance correction. Change the suffixes
for bands 4, 5, and 6 to “_darkObj”. Click OK.

The Dark-object subtraction algorithm uses the input radiances and sun angles, as specified in the MTL text file, along with the lowest
non-zero DN in an image to remove atmospheric haze. While this image is mostly devoid of haze, contrast will be improved.

EXERCISE 3-4 IMAGE RESTORATION: LANDSAT 8-9 218

L Open the reflectance with no correction output for band 4 with the grayscale and compare it with the output using Dark-
object subtraction. Choose the “stretch current view” option from Composer for each image output.

Note that while the images look similar, the range of values is greater for the atmospherically corrected image, indicating better
contrast. Even when haze is minimal in an image, atmospheric correction is suggested to improve contrast and standardize
reflectance.

Next, we will investigate the thermal imagery captured by Landsat 8.

M In the LANDSAT module, select not to include Bands 4, 5 and 6 from the include column, but set bands 10 and 11 to “yes”.
The Thermal band options will now become available. Leave the default, Raw DN, and change the file name suffixes for
bands 10 and 11 to “_raw”. Click OK. When the module finishes running, click the Instant stretch button in Composer to
improve the contrast.

Note that this image looks the opposite than the visual imagery we imported above, with the lowest values (lowest emission) from
cold surfaces (ice and snow) and highest from warmer surfaces (areas in the rain shadow). As with the other raw imagery, values are
in the native 16-bit format. However, they have no real physical meaning. To remedy this, we can import the thermal imagery and
convert to at-satellite brightness temperature.

N From the LANDSAT module, select the “Convert to at-satellite brightness temperature” option under the thermal bands
options. Change the output suffixes to “_atSat”. Click OK. Once the LANDSAT module is done running, click on the image to
make it active, then in Composer, select the stretch icon “stretch current view”. Then from Layer Properties specify a new
palette RADAR.

Raw DN values have now been converted to at-satellite brightness temperatures, in degrees Kelvin. To get a better idea of what these
values mean, we will zoom in to a subsection of the image.

O Using the zoom window option, select a box around the right-hand side of the image, about 1/3 from the bottom. Maximize
your map window. Use the identify tool to select and investigate pixel values.

EXERCISE 3-4 IMAGE RESTORATION: LANDSAT 8-9 219

Identifying the darkest pixels in the upper-left hand corner will return values near 273 Kelvin, or approximately 0 Celsius, a value we
would expect given that this is the surface of an icefield in midsummer. Identifying the greenish pixels near the bottom of the image
returns values between 281 and 282, or 7 or 8 Celsius. The difference in water temperature between the Klinaklini River (draining
from the north from the Ha-Iltzuk Icefield in the upper left) and the ocean water in Knight Inlet is apparent. Note also the high values
on some south facing slopes, reflecting the relative warmth of the exposed rock of the mountainsides in the area. Converting to at-
satellite brightness temperatures accurately reflects expected temperatures of the surfaces in this image.

EXERCISE 3-4 IMAGE RESTORATION: LANDSAT 8-9 220

‬ EXERCISE 3-5
PRINCIPAL COMPONENTS ANALYSIS
FOR MULTI-SPECTRAL IMAGERY

In a priro exercise, we explored the use of principal components analysis for the removal of noise in multi-spectral imagery. Although
PCA is commonly used for this purpose, in this exercise, we will explore its additional wide use as a method for data compaction. In
satellite imagery, it is not uncommon to find that a strong degree of correlation exists between the multispectral bands. Such
correlation indicates that if reflectances are high at a particular location on one band, they are also likely to be high on the other
band. In the extreme case, if two bands were perfectly correlated, they would essentially describe the same information. It is not
unusual to find that an image with 7 bands, such as Landsat Thematic Mapper, contains far fewer than 7 bands of true information.

The question then arises as to whether just a few of the bands provide an adequate characterization of earth surface reflectances. To
answer this, let's explore the information-carrying characteristics of the Landsat imagery we used in a prior exercise through
Principal Components Analysis..

Principal Components Analysis (PCA) is related to Factor Analysis and can be used to transform a set of image bands such that the
new bands (called components) are uncorrelated with one another and are ordered in terms of the amount of image variation they
can explain. The components are thus a statistical abstraction of the variability inherent in the original band set.

Since each of the components produced by this transformation is uncorrelated with the other, each carries new information. Also,
because they are ordered in terms of the amount of information they carry, the first few components will tend to carry most of the
real information in the original band set while the later components will describe only minor variations. One application of Principal
Components Analysis then is data compaction—by retaining only the first few components, one can keep most of the information
while discarding a large proportion of the data.

With high processor speeds and disk capacity, data compaction is less of an issue now than it was in the past. Most classifiers will
allow the input of many bands, and it is common to use all bands in a classification, whether they are highly correlated or not.
However, reducing the number of bands may increase efficiency since noise may possibly be eliminated and the classifiers have less
information to discriminate. In this exercise, through PCA, we will learn about the information-carrying characteristics of our Landsat
data. 1

A Set the working folder back to the Introductory IP tutorial folder. Display H87TM4 (the near infrared band) with the
Greyscale palette and Equal Interval autoscaling. Use the Instant Stretch tool (found on Composer) as well. Now display
each of the remaining bands in the same way.

1 For more detail on the PCA methodology, please review the Help for the PCA module.

EXERCISE 3-5 PRINCIPAL COMPONENTS ANALYSIS FOR MULTI-SPECTRAL IMAGERY

221
1 Do any other bands look like band 4 (H87TM4)? Which one(s)?

B Now run PCA (Principal Components Analysis) from the Image Processing/Transformation menu. Choose Forward t-mode
and the covariance matrix unstandardized option. Indicate that seven bands will be used. Click into the Image Band Name
list, click the Pick List button, and choose H87TM1. Do the same for each of the seven bands. Indicate that seven
components should be extracted. Alternatively, you can use the Insert layer group option and select the H87 raster group
file. Enter H87 as the new prefix for the output files and the complete text output option.

PCA will then proceed to calculate the transformation equations and write out the new component files with names that
range from H87_T-MODE_CMP1 through H87_T-MODE_CMP7.

The results will appear on the screen as summary tables when the PCA module has finished working. You may print this if
you wish.

2 Look at the correlation matrix. Is there much correlation between bands? Which band correlates most with band 1? Do any
bands correlate with band 4? How does this compare to your answer for question 1?

C Now scroll down the screen to look at the component summary table where the eigenvalues and eigenvectors for each
component (listed as columns) are displayed. The eigenvalues express the amount of variance explained by each
component and the eigenvectors are the transformation equations. Notice that this has been summarized as a percent
variance explained (% var.) measure at the top of each column.

3 How much variance is explained by components 1, 2 and 3 separately? How much is explained by components 1 and 2
together (add the amount explained by each)? How much is explained by components 1, 2 and 3 together?

D Now scroll down the screen and look at the table of loadings. The loadings refer to the degree of correlation between these
new components (the columns) and the original bands (the rows).

4 Which band has the highest correlation with component 1? Is it a high correlation?

5 Which band has the highest correlation with component 2?

If you did not print the tables, do not close this window since you will need to refer to this information later. Merely minimize it to
make room for the display of other images.

E Now display the following four images, all with Equal Interval autoscaling and the Greyscale palette: H87_T-
MODE_CMP1(component 1), H87_T-MODE_TM2 (band 2), H87_T-MODE_CMP3 (component 3), H87_T-MODE_TM4 (band
4).

EXERCISE 3-5 PRINCIPAL COMPONENTS ANALYSIS FOR MULTI-SPECTRAL IMAGERY

222
 Try arranging all these images on the screen at the same time so all are visible. Remember, you can reduce the size of the
layer frame by double-clicking in the image, dragging one of the sizing handles, clicking outside the image, then clicking
the Fit Map Window to Layer Frame toolbar icon.

6 How similar does component 1 look to the infrared image? How similar does component 2 look to the red image?

F Now look at component 7 (H87_T-MODE_CMP7) with the autoscale option.

7 How well does this correlate with the original seven bands (use the loadings chart to determine this)? Judging by what you
see, what do you think is contained in component 7? How much information will be lost if you discard this component?

The relationships we see in this example will not be the same in every landscape. However, this is not an uncommon experience. If
you had to choose only one band to work with, it is often the case that the near infrared band (TM band 4) carries the greatest
amount of information. After this, it is commonly the red visible band that carries the next greatest degree of information. After this it
will vary. However, the green visible (TM band 2) and middle infrared (TM band 5) bands are two good candidates for a third band
to consider.

Going back to our original question, it is clear that three bands can carry an enormous amount of information. In addition, we can
also see that the bands that are used in the traditional false color composite (green, red and infrared) are also very well chosen—
they clearly carry the bulk of the information in the full data set. Thus, for the purpose of unsupervised classification, which we will
explore in the next exercise, it makes sense that we could use just three bands of imagery to carry out the image classification.

You may delete the seven component images (H87_T-MODE_CMP1-7).

EXERCISE 3-5 PRINCIPAL COMPONENTS ANALYSIS FOR MULTI-SPECTRAL IMAGERY

223
‬ EXERCISE 3-6
SUPERVISED CLASSIFICATION
In an earlier exercise of this section, we drew the spectral response patterns for three kinds of land covers: urban, forest and water.
We saw that the spectral response patterns of each of these cover types were unique. Land covers, then, may be identified and
differentiated from each other by their unique spectral response patterns. This is the logic behind image classification. Many kinds of
maps, including land cover, soils, and bathymetric maps, may be developed from the classification of remotely sensed imagery.

There are two methods of image classification: supervised and unsupervised. With supervised classification, the user develops the
spectral signatures of known categories, such as urban and forest, and then the software assigns each pixel in the image to the cover
type to which its signature is most similar. With unsupervised classification, the software groups pixels into categories of like
signatures, and then the user identifies what cover types those categories represent.

The steps for supervised classification may be summarized as follows:

1. Locate representative examples of each cover type that can be identified in the image (called training sites).

2. Digitize polygons around each training site, assigning a unique identifier to each cover type.

3. Analyze the pixels within the training sites and create spectral signatures for each of the cover types.

4. Classify the entire image by considering each pixel, one by one, comparing its signature with each of the known signatures.
So-called hard classifications result from assigning each pixel to the cover type that has the most similar signature. Soft
classifications, on the other hand, evaluate the degree of membership of the pixel in all classes under consideration,
including unknown and unspecified classes. Decisions about how similar signatures are to each other are made through
statistical analyses. There are several different statistical techniques that may be used. These are often called classifiers.

This exercise illustrates the hard-supervised classification techniques available in TerrSet. Soft classifiers are explored in the
Advanced Image Processing Exercises of the Tutorial. A more detailed discussion of both types of classification may be found in the
section Classification of Remotely Sensed Imagery in the TerrSet Manual.

Training Site Development

We will begin by creating the training sites. The area we will classify is a small, windowed area around Howe Hill, immediately
northwest of the airport, that we saw in the HOW87TM1-4 images in an earlier exercise. Figure 1 shows the results of a field visit to
this area. The training sites created in this exercise will be based on the knowledge of land cover types identified during this visit.

EXERCISE 3-6 SUPERVISED CLASSIFICATION 224

 Urban (Streets) Agriculture

Conifers

Deep
Water

Shallow
Water

Deciduous

Conifers

Urban
(Abandoned
Airport)

Conifers

Deep Water
Deciduous
Figure 1

Each known land cover type will be assigned a unique integer identifier, and one or more training sites will be identified for each.

A Write down a list of all the land cover types identified in Figure 1, along with a unique identifier that will signify each cover
type. While the training sites can be digitized in any order, they may not skip any number in the series, so if you have ten
different land-cover classes, for example, your identifiers must be 1 to 10.

The suggested order (to create a logical legend category order) is:

1-Shallow water
2-Deep water
3-Agriculture

EXERCISE 3-6 SUPERVISED CLASSIFICATION 225

 4-Urban
5-Deciduous Forest
6-Coniferous Forest

B Display the image called H87TM4 using the Greyscale palette, with autoscaling set to Equal Intervals. Use the on-screen
digitizing feature of TerrSet to digitize polygons around your training sites. On-screen digitizing in TerrSet is made
available through the following three toolbar icons:

Digitize Delete Feature Save Digitized Data

C Zoom into the image to focus closely around the deep water lake at the left side of the image. Then select the Digitize icon
from the toolbar.

Enter TRAININGSITES as the name of the layer to be created. Use the Default Qualitative palette and choose to create
polygons. Enter the feature identifier you chose for deep water (e.g., 2). Press OK.

The vector polygon layer TRAININGSITES is automatically added to the composition and is listed on Composer. Your
cursor will now appear as the digitize icon when in the image. Move the cursor to a starting point for the boundary of your
training site and press the left mouse button. Then move the cursor to the next point along the boundary and press the left
mouse button again (you will see the boundary line start to form). The training site polygon should enclose a
homogeneous area of the cover type, so avoid including the shoreline in this deep water polygon. Continue digitizing until
just before you have finished the boundary, and then press the right mouse button. This will finish the digitizing for that
training site and ensure that the boundary closes perfectly. The finished polygon is displayed with the symbol that
matches its identifier.

You can save your digitized work at any time by pressing the Save Digitized Data icon on the toolbar. Answer yes when
asked if you wish to save changes.

If you make a mistake and wish to delete a polygon, select the Delete Feature icon (next to Digitize). Select the polygon
you wish to delete, then press the delete key on the keyboard. Answer yes when asked whether to delete the feature. You
may delete features either before or after they have been saved.

Use the navigation tools to zoom back out, then focus in on your next training site area, referring to Figure 1. Select the
Digitize icon again. Indicate that you wish to add features to the currently active vector layer. Enter an identifier for the
new site. Keep the same identifier if you want to digitize another polygon around the same cover type. Otherwise, enter a
new identifier.

Any number of training sites, or polygons with the same ID, may be created for each cover type. In total, however, there should be an
adequate sample of pixels for each cover type for statistical characterization. A general rule of thumb is that the number of pixels in
each training set (i.e., all the training sites for a single land cover class) should not be less than ten times the number of bands. Thus,
in this exercise, where we will use seven bands in classification, we should aim to have no less than 70 pixels per training set.

D Continue until you have training sites digitized for each different land cover. Then save the file using the Save Digitized
Data icon from the toolbar.

Signature Development

EXERCISE 3-6 SUPERVISED CLASSIFICATION 226

After you have a training site vector file you are ready for the third step in the process, which is to create the signature files. Signature
files contain statistical information about the reflectance values of the pixels within the training sites for each class.

E Run MAKESIG from the Image Processing/Signature Development menu. Choose Vector as the training site file type and
enter TRAININGSITES as the file defining the training sites. Click the Enter Signature Filenames button. A separate
signature file will be created for each identifier in the training site vector file. Enter a signature filename for each identifier
shown (e.g., if your shallow water training sites were assigned ID 1, then you might enter SHALLOW WATER as the
signature filename for ID 1). When you have entered all the filenames, press OK.

Indicate that seven bands of imagery will be processed by pressing the up arrow on the spin button until the number 7 is
shown. This will cause seven input name boxes to appear in the grid. Click the Pick List button in the first box and choose
H87TM1 (the blue band). Click OK, then click the mouse into the second input box. The pick button will now appear on
that box. Select it and choose H87TM2 (the green band). Enter the names of the other bands in the same way: H87TM3
(red band), H87TM4 (near infrared band), H87TM5 (middle infrared band), H87TM6 (thermal infrared band) and
H87TM7 (middle infrared band). Click OK.

F
When MAKESIG has finished, open the TerrSet Explorer from the File menu. Select the filter to display the signature file
type (sig) and check that a signature exists for each of the six land cover types. If you forgot any, repeat the process
described above to create a new training site vector file (for the forgotten cover type only) and run MAKESIG again.

To facilitate the use of several subsequent modules with this set of signatures, we may wish to create a signature group file. Using
group files (instead of specifying each signature individually) quickens the process of filling in the input information into module
dialog boxes. Similar to a raster image group file, a signature group file is an ASCII file that may be created or modified with TerrSet
Explorer. MAKESIG automatically creates a signature group file that contains all our signature filenames. This file has the same name
as the training site file, TRAININGSITES.

G Open TerrSet Explorer from the File menu. From the Filters pane select to also display signature and signature group files.
Then from the Files tab choose TRAININGSITES. In the Metadata pane, verify that all the signatures are listed in the group
file.

To compare these signatures, we can graph them, just as we did by hand in the previous exercise.

H Run SIGCOMP from the Image Processing/Signature Development menu. Choose to use a signature group file and choose
TRAININGSITES. Display their means.

1 Of the seven bands of imagery, which bands differentiate vegetative covers the best?

I Close the SIGCOMP graph, then run SIGCOMP again. This time choose to view only 2 signatures and enter the urban and
the conifer signature files. Indicate that you want to view their minimum, maximum, and mean values. Notice that the
reflectance values of these signatures overlap to varying degrees across the bands. This is a source of spectral confusion
between cover types.

2 Which of the two signatures has the most variation in reflectance values (widest range of values) in all of the bands? Why?

Another way to evaluate signatures is by overlaying them on a two-band scatterplot or scattergram. The scattergram plots the
positions of all pixels on two bands, where reflectance of one band forms the X axis and reflectance of the other band forms the Y axis.

EXERCISE 3-6 SUPERVISED CLASSIFICATION 227

The frequency of pixels at each X,Y position is signified by a quantitative palette color. Signature characteristics can be overlayed on
the scattergram to give the analyst a sense of how well they distinguish between the cover types in the two bands that are plotted.

To create such a display in TerrSet, use the module SCATTER. It uses two image bands as X and Y axes to graph relative pixel
positions according to their values in these two bands. In addition, it creates a vector file of the rectangular boundary around the
signature mean in each band that is equal to two standard deviations from this mean. Typically, one would create and examine
several scattergrams using different pairs of bands. Here we will create one scattergram using the red and near infrared bands.

J Run SCATTER from the Image Processing/Signature Development menu. Indicate H87TM3 (the red band) as the Y axis
and H87TM4 (the near infrared band) as the X axis. Give the output the name SCATTER and retain the default logarithm
count. Choose to create a signature plot file and enter the name of the signature group file TRAININGSITES. Press OK.

K
Move the cursor around in the scatterplot. Note that the X and Y coordinates shown in the status bar are the X and Y
coordinates in the scatterplot. The X and Y axes for the plot are always set to the range 0-255. Since the range of values in
H87TM3 is 12-66 and that for H87TM4 is 5-136, all the pixels are plotting in the lower-left quadrant of the scatterplot.
Zoom in on the lower-left corner to see the plot and signature boundaries better. You may also wish to click on the Full
Extend Maximized icon on the toolbar (or press the End key) to enlarge the display.

The values in the SCATTER image represent densities (log of frequency) of pixels, i.e., the higher palette colors indicate
many pixels with the same combination of reflectance values on the two bands and the lower palette colors indicate few
pixels with the same reflectance combination. Overlapping signature boxes show areas where different signatures have
similar values. SCATTER is useful for evaluating the quality of one's signatures. Some signatures overlap because of the
inadequate nature of the definition of land cover classes. Overlap can also indicate mistakes in the definition of the
training sites. Finally, overlap is also likely to occur because certain objects truly share common reflectance patterns in
some bands (e.g., hardwoods and forested wetlands).

It is not uncommon to go through several iterations of training site adjustment, signature development, and signature evaluation
before achieving satisfactory signatures. For this exercise, we will assume our signatures are adequate and will continue on with the
classification.

Classification
Now that we have satisfactory signature files for all of our land cover classes, we are ready for the last step in the classification
process—to classify the images based on these signature files. Each pixel in the study area has a value in each of the seven bands of
imagery (H87TM1-7). As mentioned above, these are respectively: blue, green, red, near-infrared, middle infrared, thermal infrared
and another middle infrared band. These values form a unique signature which can be compared to each of the signature files we just
created. The pixel is then assigned to the cover type that has the most similar signature. There are several different statistical
techniques that can be used to evaluate how similar signatures are to each other. These statistical techniques are called classifiers.
We will create classified images with three of the hard classifiers that are available in TerrSet. Exercises illustrating the use of soft
classifiers and hardeners may be found in the Advanced Image Processing section of the Tutorial.

EXERCISE 3-6 SUPERVISED CLASSIFICATION 228

L We will be producing a variety of classified images. To make the automatic display of these images more informative, open
User Preferences from the File menu and on the Display Settings tab check on the option to automatically show the legend
(in addition to the title).

The first classifier we will use is a minimum distance to means classifier. This
classifier calculates the distance of a pixel's reflectance values to the spectral
mean of each signature file, and then assigns the pixel to the category with the
closest mean. There are two choices on how to calculate distance with this
classifier. The first calculates the Euclidean, or raw, distance from the pixel's
reflectance values to each category's spectral mean. This concept is illustrated in
two dimensions (as if the spectral signature were made from only two bands) in
Figure 2. 1 In this heuristic diagram, the signature reflectance values are
indicated with lower case letters, the pixels that are being compared to the
signatures are indicated with numbers, and the spectral means are indicated
with dots. Pixel 1 is closest to the corn (c's) signature's mean, and is therefore
assigned to the corn category. The drawback for this classifier is illustrated by
pixel 2, which is closest to the mean for sand (s's) even though it appears to fall
within the range of reflectances more likely to be urban (u's). In other words,
the raw minimum distance to mean does not consider the spread of reflectance
values about the mean.

M All of the classifiers we will explore in this exercise may be found in the Image Processing/Hard Classifiers menu. Run
MINDIST (the minimum distance to means classifier) and indicate that you will use the raw distances and an infinite
maximum search distance. Click on the Insert Signature Group button and choose the TRAININGSITES signature group file.
The signature names will appear in the corresponding input boxes in the order specified in the group file. Call the output
file MINDISTRAW. Click OK to start the classification. Examine the resulting land cover image. (Change the palette to Qual
if necessary.)

We will try the minimum distance to means classifier again, but this time
with the second kind of distance calculation—normalized distances. In this
case, the classifier will evaluate the standard deviations of reflectance
values about the mean—creating contours of standard deviations. It then
assigns a given pixel to the closest category in terms of standard deviations.
We can see in Figure 3 that pixel 2 would be correctly assigned to the urban
category because it is two standard deviations from the urban mean, while
it is at least three standard deviations from the mean for sand.

N To illustrate this method, run MINDIST again. Fill out the dialog
box in the same way as before, except choose the normalized
option, and call the result MINDISTNORMAL.

3 Compare the two results. How would you describe the effect of
standardizing the distances with the minimum distance to means classifier?

1 Figures 2-5 are adapted from Lillesand and Kiefer, 1979. Remote Sensing and Image Interpretation. First edition. New York, Chichester, Brisbane and Toronto:
John Wiley & Sons.

EXERCISE 3-6 SUPERVISED CLASSIFICATION 229

The next classifier we will use is the maximum likelihood classifier. Here, the
distribution of reflectance values in a training site is described by a
probability density function, developed based on Bayesian statistics (Figure
4). This classifier evaluates the probability that a given pixel will belong to a
category and classifies the pixel to the category with the highest probability
of membership.

O
Run MAXLIKE. Choose to use equal prior probabilities for each
signature. Click the Insert Signature Group button, then choose the
signature group file TRAININGSITES. The input grid will then
automatically fill. Leave the minimum likelihood at 0.0 and call the
output image MAXLIKE, then click OK. Maximum likelihood is the
slowest of the techniques, but if the training sites are good, it tends
to be the most accurate.

Finally, we will look at the parallelepiped classifier. This classifier creates 'boxes' (i.e., parallelepipeds) using minimum and
maximum reflectance values or standard deviation units (z-scores) within the training sites. If a given pixel falls within a signature's
'box,' it is assigned to that category. This is the simplest and fastest of classifiers and the option using Min/Max values was used as a
quick-look classifier years ago when computer speed was quite slow. It is prone, however, to incorrect classifications. Due to the
correlation of information in the spectral bands, pixels tend to cluster into cigar- or zeppelin-shaped clouds. As illustrated in Figure 5,
the 'boxes' become too encompassing and capture pixels that probably should be
assigned to other categories. In this case, pixel 1 will be classified as deciduous
(d's) while it should probably be classified as corn. Also, the 'boxes' often overlap.
Pixels with values that fall at this overlap are assigned to the last signature,
according to the order in which they were entered.

Run PIPED and choose the Min/Max option. Click the Insert Signature Group
button and choose TRAININGSITES. Call the output image PIPEDMINMAX. Then
click OK. Note the zero-value pixels in the output image. These pixels did not fit
within the Min/Max range of any training set and were thus assigned a category of
zero.

The parallelepiped classifier, when used with minimum and maximum values, is
extremely sensitive to outlying values in the signatures. To mediate this, a second
option is offered for this classifier that uses z-scores rather than raw values to
construct the parallelepipeds.

P Run PIPED exactly as before, only this time choose the z-score option,
and retain the default 1.96 units. This will construct boxes that include 95% of the signature pixels. Call this new image
PIPEDZ.

4 How much did using standard deviations instead of minimum and maximum values affect the parallelepiped classification?

The final supervised classification module we will explore in this exercise is FISHER. This classifier is based on linear discriminant
analysis.

Q Run FISHER. Insert the signature group file TRAININGSITES. Call the output image FISHER. Press OK.

EXERCISE 3-6 SUPERVISED CLASSIFICATION 230

R Compare each of the classifications you created: MINDISTRAW, MINDISTNORMAL, MAXLIKE, PIPEDMINMAX and PIPEDZ.
To do this, display all of them with the Default Qualitative palette. You may need to make the window frames smaller to fit
all of them on the screen.

5 Which classification is best?

As a final note, consider the following. If your training sites are very good, the Maximum Likelihood or FISHER classifiers should
produce the best results. However, when training sites are not well defined, the Minimum Distance classifier with the standardized
distances option often performs much better. The Parallelepiped classifier with the standard deviation option also performs rather
well and is the fastest of the considered classifiers.

Keep MINDISTNORMAL and MAXLIKE for use in Exercise 3-7. You may delete the other images created in this exercise.

EXERCISE 3-6 SUPERVISED CLASSIFICATION 231

‬ EXERCISE 3-7
UNSUPERVISED CLASSIFICATION

Unsupervised classification is another technique for image classification. In the unsupervised approach, the dominant spectral
response patterns that occur within an image are extracted, and the desired information classes are identified by means of ground
truthing. In TerrSet, unsupervised classification is provided by way of two modules named CLUSTER and ISOCLUST. This exercise
focuses on CLUSTER.

CLUSTER uses a histogram peak selection technique. This is equivalent to searching for the peaks in a one-dimensional histogram,
where a peak is defined as a value with a greater frequency than its neighbors on either side. Once the peaks have been identified, all
possible values are assigned to the nearest peak. Thus, the divisions between classes tend to fall at the midpoints between peaks.
Because this technique has specific criteria for what constitutes a peak, you do not need to make a prior estimate (as some
techniques require) of the number of clusters an image contains—it will determine this for you.

CLUSTER evaluates a multi-dimensional histogram based on the number of input bands. We will use all seven bands used in the
previous exercise of TM for the Howe Hill area to illustrate this technique.

A
To facilitate data entry, we will create a raster group file for all these bands. Open TerrSet Explorer and from the File pane
select six of the seven H87TM bands. Do not include band six, H87TM6. With the six bands highlighted, right-click and
select Create Raster Group. By default, a file named RASTER [Link] is created. Select this file, right-click and rename
it to HOWEHILL.

B
Now run CLUSTER from the Image Processing/Hard Classifiers menu. Choose to insert the layer group, HOWEHILL. This
will insert all six bands into the input filename grid. Call the output image BROAD. Then choose the broad generalization
level and elect to drop the least significant clusters with 10%. Leave the option Grey levels at its default of 6. The result,
BROAD, will be displayed with the Qualitative color palette.

C To facilitate visual analysis of this image, you may wish to use the category "flash" option. Place the cursor over a legend
color box and press and hold down the left mouse button. This will cause that category to be displayed in red, while every
other category is displayed in black. When you release the mouse button, the display will return to normal.

You can also display three of the H87TM bands behind the broad classification as a composite. With the broad
classification image in focus, add the three raster layers to the composition: H87TM2, H87TM3, H7TM4. Then select each
band and assign an RGB component, using the icons in Composer. Select H87TM2 and assign it the Blue component. Select
H87TM3 and assign it the Green component. And select H87TM4 and assign it the Red component. Once the Red
component is assigned, the false color composite will obscure the broad classification. Move the file BROAD so it is the top
layer (bottom of the file list in Composer). Then, by clicking on and off the file BROAD, you can investigate your
assumptions of the classes.

The result from CLUSTER is an image of the very broad spectral classes in the study area.

EXERCISE 3-7 UNSUPERVISED CLASSIFICATION 232

1 How many broad clusters were produced? Given your knowledge of the area from the supervised classification exercise, what
land cover do you think is represented by each of the clusters?

The broad and fine generalization levels use different decision rules when evaluating the frequency histogram for peaks. In broad
clustering, a peak must contain a frequency higher than all of its non-diagonal neighbors. Fine classification allows a peak to have one
non-diagonal neighbor with a higher frequency. This accommodates true peaks which are otherwise missed because nearby peaks of
greater magnitude obscure the usual dip between the peaks. This concept is illustrated in one-dimensional space in Figure 1. Broad
clusters are divided only at the valleys. Fine clusters are divided at both the valleys and the shoulders of the histogram.

D Use CLUSTER again, with the same six H87TM images to create an image called FINE. This time, use the fine generalization
level, and again, elect to drop the 10% least significant clusters. As you can see, the fine generalization produces many
more clusters. Scroll down the legend or increase the size of the legend box to see how many clusters there are.

2 How many clusters are produced? Which cluster is most easily identified? Why do you think this is the case?

E
Image histograms allow us to see the difference in the distribution of pixels among classes, depending on the
generalization level. Run HISTO from the Display menu to create a histogram of FINE, keep the rest of the defaults. In the
output of the CLUSTER module, Cluster 1 is always the one with the highest frequency of pixels. It corresponds to the
largest land cover type detected during classification. The second cluster has a smaller number of pixels and so on.

Note that many of the higher numbered clusters have relatively few pixels. One approach that is often employed is to look for a
natural break in the histogram of fine clusters to estimate the number of significant cover types in the study area. Once determined,
you can run the CLUSTER module again, this time specifying the number of clusters to identify. All remaining pixels are assigned to
the cluster to which they are most similar. (Note that this would not be a good approach if you were specifically looking for a land
cover type that covers a little area.)

F Look at the histogram of FINE. Note that the study area is dominated by two clusters. Several small breaks in the
histogram might be chosen as the cutoff point. One might choose to set the number of clusters to 6, 10 or 15 based on
those breaks in the histogram. For ease of interpretation in the absence of ground truth information, we will choose to
keep the first 6 clusters as our significant land cover types.

G Run CLUSTER with our six bands again. This time give FINE10 as the output filename, choose the fine generalization level
and choose to set the maximum number of clusters to 10. Keep the remaining defaults.

EXERCISE 3-7 UNSUPERVISED CLASSIFICATION 233

The problem we now face is how to interpret these clusters. If you know a region, the broad clusters are often easy to interpret. The
fine clusters, however, can take a considerable amount of care in interpreting. Usually existing maps, aerial photographs and ground
visits are required to adequately identify these finer clusters. In addition, we will often find that we need to merge certain clusters to
produce our final map. For example, we might find that one cluster represents pine forest on shaded slopes while another is the pine
forest on bright slopes. These are two distinct spectral classes. In the final map, we want both to be part of a single pine forest
information class. To group and reassign clusters like this, we can use ASSIGN.

H
Try to interpret the 10 clusters of FINE10. To do so, compare FINE10 with the supervised classification outputs you
created in the previous exercise (MINDISTNORMAL and MAXLIKE). You may also find it useful to look at the original
bands or composite images (create 24-bit composites for a better visual effect) to determine what cover type is
represented by a cluster. When you have determined to which category each cluster should be assigned, use Edit to enter
this information into an attribute values file called LANDCOVER. The cluster numbers should be listed in the first column
and the numeric land cover categories in the second column of the values file. Accept the default integer data type when
asked.

3 What were your class assignments?

I Use ASSIGN to create the new land cover image. The feature definition file is FINE10, the values file is LANDCOVER and call
the output image LANDCOVER. Display it with the Qualitative palette. Use the Metadata utility in TerrSet Explorer to add
meaningful legend captions to LANDCOVER and save. Then redisplay LANDCOVER to cause the new legend information to
appear in the display.

The unsupervised cluster classification is a very quick way to gain knowledge of the study area. Classification is most often an
iterative process where each step yields new information that the analyst can use to improve the classification. Oftentimes,
supervised and unsupervised classifications are used together in hybrid approaches. For example, in FINE10, cluster number 3 is
quite difficult to interpret, yet it is the third most prevalent spectral class in the study area. This might alert us to a land cover
category (e.g. wetlands) that was left out of the original set of cover classes we developed signatures for in the supervised
classification. We could then go back and create a training site and signature for that class and re-classify the image using the
supervised classifiers. The clusters of an unsupervised analysis might also be used as training sites for signature development in a
subsequent supervised classification. The important thing to note is that classification is hardly ever a single-step process.

Finally, no classification is complete without an accuracy assessment. Accuracy assessment provides the means to assess the
confidence with which one might use the classified land cover map. It can also provide information to help improve the classified
map. The Classification of Remotely Sensed Imagery section in the TerrSet Manual describes this important process.

In this set of exercises, we have concentrated on the hard classifiers. The soft classifiers, which delay the assignment of each pixel to a
class, are described in the Advanced Image Processing set of exercises in this Tutorial.

EXERCISE 3-7 UNSUPERVISED CLASSIFICATION 234

‬ EXERCISE 3-8
CHANGE ANALYSIS - PAIRWISE AND
MULTIPLE IMAGE COMPARISON

This exercise will explore some of the ways in which environmental change can be analyzed through image comparison. Explanations
of the techniques that are used can be found in the Change Analysis chapter of the TerrSet Manual and in Lillesand et. al. 2004. 1 The
techniques in this exercise relate to quantitative pairwise and multiple image data only and include simple differencing, thresholding,
image regression, image ratioing, and change vector analysis. Subsequent exercises deal with qualitative image comparison using
Land Change Modeler and the analysis of long time series of quantitative images using Earth Trends Modeler. Bear in mind that while
tools are available to analyze change, there are no standard procedures for applying those tools. As a result, this exercise should be
regarded as an exploration, not as a definitive approach or set of steps.

Simple Differencing
The first technique explores differences in the quantitative distribution of vegetation over the continent of Africa for the same month
over two different years. The first image is a normalized difference vegetation index (NDVI) image derived from NOAA (United States
National Oceanic and Atmospheric Administration) AVHRR (Advanced Very High Resolution Radiometer) satellite imagery for the
month of December 1987 (called AFDEC87). The second is a corresponding image for December 1988 (called AFDEC88). Has any
significant change in vegetation occurred between the two years? If so, what areas are affected?

NDVI is a quantitative measure that correlates highly with the quantity of living vegetative matter in any region. The index is derived
quite simply using the red and near infrared wavelength bands of AVHRR (or any other source) data. In green vegetation, the
presence of chlorophyll causes strong absorption of red wavelengths while leaf structure will tend to cause high reflectance of near
infrared wavelengths. As a result, areas with a dense and vigorous vegetative canopy will tend to show a strong contrast in the
reflectances in these two regions. The index is calculated as follows:

NDVI = (Infrared - Red) / (Infrared + Red)

This operation is available directly in the OVERLAY and VEGINDEX modules of TerrSet and can be duplicated on most systems using
the simple math operators provided. In our case, however, these images were processed directly by NOAA and rescaled to a byte
integer range (i.e., the images measure NDVI directly with a range of values between 0-255).

1 Lillesand, T. M., R.W. Kiefer, and J.W. Chipman. 2004. Remote Sensing and Image Interpretation. John Wiley & Sons.

EXERCISE 3-8 CHANGE ANALYSIS – PAIRWISE AND MULTIPLE IMAGE COMPARISON 235
A First, open TerrSet and from TerrSet Explorer, verify that the Introductory IP folder is listed as either the main Working
Folder or a Resource Folder.

B Next, display the AFDEC87 and AFDEC88 images with DISPLAY Launcher using the NDVI palette and the Equal Intervals
autoscale option.

In these images, low NDVI values are shown in brown colors while high NDVI values are shown in dark green.

1 What are the main differences you can identify through visual comparison?

We will now create a simple difference image to compare the two dates by subtracting the 1987 image from the 1988 image. There
are several ways to accomplish this in TerrSet. We will use the Image Calculator. This facility allows us to use entire images as
arguments in mathematical equations. Operations are performed between corresponding pixels of the input images to produce an
output image. (Image Calculator makes calls to other modules such as OVERLAY, SCALAR and TRANSFORM. You will see these in the
status bar as Image Calculator evaluates expressions.)

C Open Image Calculator. Enter the output filename DIFF88-87 in the first input box. Place your cursor in the Expression to
process input box and click the Insert Image button. Select the image, AFDEC88, from the Pick List. Click on the subtraction
button, then use Insert Image again and select AFDEC87 (see figure below). Then click the Process Expression button.
When the calculation is completed, the output image will automatically display. If necessary, change the display palette to
NDVI from within Composer.

AFDEC88 - AFDEC87 = DIFF88-87

While the image is displayed, press the Add Layer button on the Composer dialog and enter the vector filename COUNTRY
using the Outline white symbol file to overlay the country boundaries.

The legend provides information about the correspondence between image colors and data values. We can also query the
data values at particular points. To do so, make sure that the DIF88-87 image is selected in Composer, then click anywhere
in the image. The value at the cursor location is displayed. Feel free to do this with any input or output images we work
with.

2 The positive value areas on this image are those in which we have a stronger NDVI in 1988 than in 1987 while the negative
value areas are those in which the NDVI is lower. For the sake of discussion, we will call these areas positive and negative
change. What areas have strong positive change? What areas have strong negative change?

D You may find it helpful in the visual analysis to isolate the positive and negative change areas. To do this, make sure
DIFF88-87 is in focus (by clicking anywhere in the image), then click Layer Properties on Composer. The contrast settings
allow you to interactively control the saturation points of the display. To highlight the negative change areas in the image,
set the display maximum endpoint to 0. This causes all pixels that have the value 0 or higher to be displayed in the highest
palette color -- green in this case. After examining the image, choose the Revert button, then set the display minimum to 0.
This causes all the pixels with values less than or equal to 0 to be displayed with the lowest palette color -- black in this

EXERCISE 3-8 CHANGE ANALYSIS – PAIRWISE AND MULTIPLE IMAGE COMPARISON 236
 case. The actual data values have not been altered at all. Feel free at any point in these exercises to use the saturation
values settings to further explore any image. Alternatively, you can use a bipolar palette to achieve the same display. With
DIFF88-87 in focus, select Layer Properties from Composer, then the Advanced Palette selection. Choose the bipolar color
logic low-high-low with the inflection point value at 0. Select the third palette choice that resembles the NDVI palette, dark
green to red. Then hit OK.

What we have created here is a simple difference image. However, whenever we work with a difference image, there is the problem
of distinguishing true change from random variation. This is usually done through a process called thresholding and is explored in
the next part of the exercise.

Thresholding
With thresholding, we try to establish upper and lower limits to normal variation beyond which we consider true change to have
occurred. To establish the threshold limits to normal variation, a histogram is usually required.

E
To display a histogram of the difference image, choose the HISTO module under the Display menu. Enter DIFF88-87 as the
input image, a class width of 1, new min and max values of -120 and 120 and the graphic histogram output type option.

3 What are the mean and standard deviation values?

If you are unfamiliar with the concept of a standard deviation, it would probably be wise to consult an introductory statistics text.
Briefly, the standard deviation is a measure of the degree of variation in a data set that can be used whenever the histogram follows a
normal distribution. A normal distribution has a bell-shaped curve with a single central peak and symmetrical tails that fall off in a
convex fashion to either side.

If the data truly are normal, then the standard deviation (often abbreviated with the Greek letter sigma -- V) measures the
characteristic dispersion of values away from the mean and can be used to evaluate the probability that certain differences from the
mean would be expected. For example, approximately 95% of all values would be expected to fall within plus or minus 2 V from the
mean while almost 99% would be expected to fall within plus or minus 3 V. Data values that are more than 3 Vfrom the mean are
very unusual (see figure).

The mean and standard deviation can thus be used to isolate unusual changes. However, in our case, the distribution is only
somewhat normal in character. Despite this, we will go ahead with this procedure.

EXERCISE 3-8 CHANGE ANALYSIS – PAIRWISE AND MULTIPLE IMAGE COMPARISON 237
To create our thresholds, we will take the mean and subtract three times the standard deviation to get the lower threshold. We will
then add three times the standard deviation to the mean to get the upper threshold. This should isolate the most unusual values that
we can call significant change.

Lower Threshold = Mean - 3 V = -44.5650

Upper Threshold = Mean + 3 V = 41.3257

F
To create the thresholded image, we will use the module RECLASS. The resulting image will have three classes -- class 1
covering all values less than -44.5650 (Mean - 3 V), class 0 covering all values from -44.5650 to 41.3257, and class 2
covering all values greater than 41.3257 (Mean + 3 V).

Open RECLASS and enter DIFF88-87 as the input file and CHG88-87 as the output file (Figure 3). Then assign a new value
of 1 to all values ranging from -75 to those just less than -44.5650. Assign a new value of 0 to all values ranging from -
44.5650 to those just less than 41.3257. Finally, assign a new value of 2 to all values ranging from 41.3257 to those just
less than 112. When finished, press OK to execute the reclassification.

Criteria

1 = <-3 V

DIFF88-87 RECLASS 0 = -3 V to +3 V = CHG88-87

2 = >+3V

G The image will automatically display with a qualitative palette in which value 0 (no significant change in our image) is
represented with the color black. If you find it difficult to distinguish the change areas on this black background, you may
find it useful to make a special palette to use with images like this. Choose Symbol Workshop under the Display menu or
from its toolbar icon. Choose File/New and select the Palette option. Enter the new filename CHANGE and press OK. Adjust
the color mixes of Red, Green, and Blue for palette colors 0, 1 and 2 such that they will be meaningful to you when used to
display change images. (You might consider, for example, light grey for no change (0), bright red for negative change (1)
and bright green for positive change (2)). Under File, choose the Save option and then exit Symbol Workshop.

To apply the new palette to the image, click Layer Properties on Composer. Choose the palette file CHANGE, then press OK
on Layer Properties. Do not autoscale the display since you want values in the image to correspond directly to the color
numbers in the CHANGE palette.

4 Where are the areas of “negative” and “positive” change on this image? (You may want to use Composer to add the layer
COUNTRY.) Does your list of significant change areas differ from your list for question 1? If so, describe this difference.

You probably noticed that the mean value of the difference image is not 0. This suggests that there is an overall change between the
two dates. One possibility is that on average, December 1988 was simply not as wet (NDVI correlates very highly with rainfall) as
December 1987. The other possibility is that the sensor on the satellite was not working identically during the two time periods. In
fact, it is not only differences in the mean that we need to be concerned about, but also differences in variability.

EXERCISE 3-8 CHANGE ANALYSIS – PAIRWISE AND MULTIPLE IMAGE COMPARISON 238
Differences in the mean and variation may be the result of such effects as sensor drift, atmospheric conditions or differences in
illumination, in which case they will lead to non-comparability of data values. The next step in this exercise will review a technique to
try to compensate for these conditions.

Image Regression
To correct for changes in the mean and variation, a technique known as image regression can be used. Regression is used to
determine the relationship between variables. If you are unfamiliar with the technique, you should probably consult an introductory
statistics text. In TerrSet, a module named REGRESS provides a simple linear regression facility for determining the relationship
between the data either in two image or two values files. In this case, we will look at the relationship between two images.

With image regression, we assume that the image at time two is a function of that at time one (i.e., that it is the way it is largely
because of the way it was in the past). The time-one image is thus the independent variable, and the time-two image is the dependent
variable. REGRESS calculates the linear relationship between the two images and plots a graph of individual pixel values using the
two dates as the X and Y axes. The regression equation can then be used with Image Calculator to create a predicted image for time
two based on the data for time one. This predicted image is really the time-one image but adjusted for overall differences in the mean
and for differences in variation about the mean. Thus, we could equally refer to the predicted time-two image as an adjusted time-one
image.

Once an adjusted time-one image has been created, it is then subtracted from the actual time-two image to yield a difference image
that can then be thresholded in the normal way. Let's try this with our data.

H We first must choose between computing the regression between the full images or on samples taken from them. If we use
the full images, AFDEC87 would be entered as the independent variable and AFDEC88 as the dependent variable.
However, with any geographic or finely spaced data, one should consider the presence of spatial autocorrelation because it
produces a false indication of the degrees of freedom in the data (a measure of the effective number of independent
sample points). For our purposes, we will use the regression only to estimate the regression coefficients that will be used
to adjust for atmospheric and instrument calibration effects. Sample spacing will not bias these estimates -- only our
consideration of their significance. One should explore more fully this concept of spatial autocorrelation before utilizing
this technique on actual data.

I
Now run the module REGRESS. Indicate that you will be computing a regression between images and specify AFDEC87 as
the independent variable and AFDEC88 as the dependent variable. We will not use the mask image option. Click OK to run.

In the REGRESS display, the frequency of pixels is indicated in the scatterplot and the best-fit line is shown. The equation
of that line is also provided and should read as follows: Y = -11.822667 + 1.222612 X with a correlation coefficient (r) of
0.945929 and a t statistic of 1337.11.

EXERCISE 3-8 CHANGE ANALYSIS – PAIRWISE AND MULTIPLE IMAGE COMPARISON 239
The equation states that the value in December 1988 is equal to -11.822667 plus 1.222612 times the value in December 1987. The
correlation coefficient (r) of 0.95 is squared to produce the coefficient of determination. This indicates that just over 89% of the
variability in December 1988 can be explained by the variability in 1987! The slope of the equation is 1.222612. You will notice that
REGRESS also provides a t statistic to test whether this slope is significantly different from 1. In our case, the value of t is very high,
suggesting (along with the value of the slope itself) that this can probably be considered a significant difference. However, it should
be noted that a definitive statistical test does require that we have confidence in the stated degrees of freedom, and that an analysis
of spatial autocorrelation would be required for a strongly defensible judgment. In our case, however, there is a significant change in
the variability (as evidenced by the slope) from one date to the next. Let's use this equation then to adjust the 1987 data.

J First, close all the open windows and displays. Then use Image Calculator to evaluate the following mathematical
expression. Remember to use the Insert Image button to add existing images to the Expression to process input box.

ADJUST87 = ([AFDEC87] * 1.222612) - 11.822667

K You may wish to change the palette for the display of ADJUST87. To do so, select Layer Properties from Composer and
enter NDVI for the palette file. Click OK.

L
Now that you have the adjusted 1987 image (or the predicted 1988, depending on how you wish to consider it), let's use it
to create a new difference image. Use Image Calculator to create an image called DIFFADJ that is the difference of AFDEC88
and ADJUST87 (see figure).

AFDEC88 – ADJUST87 = DIFFADJ

M
Now utilize the HISTO module and display a histogram of DIFFADJ. Since this is a real number image, change the minimum
and maximum to new values of -97 and 114 and choose a graphic output with a class width of 1.0.

5 How does this distribution differ from that of DIFF88-87?

EXERCISE 3-8 CHANGE ANALYSIS – PAIRWISE AND MULTIPLE IMAGE COMPARISON 240
6 What are the mean and standard deviation values for this image? How does this compare to the previous difference image?

N Now use the same thresholding procedure as before (using RECLASS) to create an image called CHGADJ (see figure) that
illustrates areas of significant change based on three standard deviations away from the mean.

Criteria

1 = <-3V

DIFFADJ RECLASS 0 = -3V to 3V = CHGADJ

2 = >3V

O At this point, we would like to compare this image with the previous change image you created. Change the palette for the
display of CHGADJ to be the CHANGE palette you created earlier (or use the default qualitative palette). Then display
CHG88-87, also with the CHANGE palette. Place the images side-by-side so you can compare them. Be sure autoscaling is
off for both. You may wish to add the vector COUNTRY layer and use the Identify tool to explore the image values.

7 How does CHGADJ compare with CHG88-87? What are the major differences?

Image regression is a very effective technique for circumventing what are known as "offset and gain" effects between images. These
effects are due to differences in the satellite sensor between the two dates. Offset refers to a shift in the mean while gain refers to a
slope that is significantly different from 1, causing values that should be identical to be different.

However, both differencing and regression differencing techniques consider differences of a given quantity to be equivalent no
matter where they occur on the measurement scale. Sometimes this is not desired. The next technique, image ratioing, provides for a
relative scaling of differences.

We will use the images CHG88-87 and CHGADJ again later when examining qualitative data comparison techniques, so do not delete
them.

Image Ratioing
In some instances, a researcher may wish to give more emphasis to differences at the low end of the scale, not unlike emphasizing the
difference between a pin dropping in a quiet room as compared to one dropping beside a running jet engine. Imagine, for example,
that a researcher is more concerned about change in arid areas than in those areas with a strong vegetative cover. In such instances a
relative scaling of differences is required and may be achieved by image ratioing. Image ratioing can be accomplished in TerrSet
using Image Calculator.

In the result, areas where the data value is identical on both images receive a value of 1.0. Those where the value is higher at time
two will have a value greater than 1.0. For instance, an area with a value two and a half times as large at time two as in time one

EXERCISE 3-8 CHANGE ANALYSIS – PAIRWISE AND MULTIPLE IMAGE COMPARISON 241
would receive a value of 2.5. Those at time two with a lower value will receive values less than 1.0, thus, for example, areas with
values half as large at time two as at time one would receive a value of 0.5. The resulting image often looks quite different from one
produced by image differencing, with change areas at the low end of the original measurement scale given substantially greater
emphasis.

There are, however, several problems with the image ratioing technique. First, the presence of zeros in the images being compared
presents a variety of problems. When the denominator is zero, the value cannot normally be evaluated because division by zero is
undefined. One solution to division by zero is to add a small increment to each image. Be aware, however, that this does affect the
scaling of the ratio. Another solution is simply to mask out from the final result all the cells that contained zeros in the denominator
image. This is an option only if the system being used allows division by zero to be performed.

TerrSet provides some mechanisms for allowing a division-by-zero operation to be completed. Zero divided by zero is evaluated as
1.0, or no change. A positive number divided by zero is evaluated as positive infinity, which is represented by a very large number (1
times 10 to the power of 18). Similarly, a negative number divided by zero is evaluated as negative infinity. Since TerrSet will allow
division by zero to occur, the kind of postprocessing discussed above can be done.

The second problem with image ratioing is that the resulting data scale is not linear. For example, while 1.0 indicates no change and
2.0 indicates twice as much in time two, 0.5 represents twice as much at time one -- a distance only half that when the sequence is
reversed! To correct this problem, Image Calculator can again be used to convert the ratio scale to a log ratio scale. The result will
then be linear and symmetrical about zero. For example, ratios of 0.5, 1 and 2 will produce log ratios of -0.69, 0 and +0.69.

To explore image ratioing, we will use a different data set because we want to illustrate one of the pairwise comparison techniques
using data at a larger scale. We have NDVI images from 1977 and 1979 derived from the Landsat MSS satellite sensor for an area of
Mauritania along the Senegal River.

The part of Mauritania in which we're interested is the Rosso area. This is located in the southwestern corner of Mauritania. Much of
Mauritania's land area is marked by plateau and desert. This Saharan zone gradually merges south into the Sahel. Further south,
along the Senegal River, there is a narrow zone of agriculture. This area is flooded seasonally and produces crops of millet, maize and
sorghum. Rainfall in this region of Mauritania in 1977 was 123.3 mm and in 1979 was 325.9 mm. (The 47-year average is 264.4 mm.)

The first two images we will work with are the NDVI images named MAUR77 and MAUR79. A normalized ratio was used, with the
infrared and red bands from Landsat MSS (Multispectral Scanner) satellite imagery. (Your exercise data includes the original four
bands of MSS imagery for both dates.)

P
Display each of the NDVI images with the NDVI palette. Click on each image in turn, select Layer Properties in Composer
and note the ranges of values.

The range of values includes negative numbers, and this is worth noting. Negative NDVI values may show up in areas in which there
is little or no vegetation. Non-vegetated areas do not display the specific spectral response of vegetation (absorption in the red band
and reflectance in the infrared band) and their NDVI ratios decrease in magnitude. (Note that whenever the red reflectance value is
higher than the infrared reflectance value, the NDVI will be negative.) Areas of snow, sand, bare soil and dead vegetation are
examples of such areas. Given that the Mauritania image covers areas with shifting dunes, the appearance of negative values is not
unusual.

Q To proceed with image ratioing, we will use this information about NDVI values to assume that pixels in our image with
negative or zero values have very little or no vegetation as measured by the satellite sensor. Therefore, we will assign all
non-positive numbers a value of .01. This will alleviate the problems of division by zero as well as that of the
interpretation of negative values. We believe that .01 represents such a low NDVI value (bare soil has been shown to have
an NDVI value of .25), that it will have little effect on our attempt to identify areas of significant "negative" and "positive"
change between the two dates. Essentially, .01 still represents the absence of vegetation. We will then ratio the images
after the zero and non-positive values have been changed.

EXERCISE 3-8 CHANGE ANALYSIS – PAIRWISE AND MULTIPLE IMAGE COMPARISON 242
R Remembering that we want the lowest value in our image to be .01, we first want to reclassify MAUR77 (using RECLASS)
so that all values just less than .01 are changed to .01. All other values will remain the same. Call this image MAUR77P (see
figure).

We now have an image that consists entirely of positive NDVI values with a minimum value of .01.

MAUR77 (Reclass .01=>.01) = MAUR77P

S Repeat the above step for MAUR79 and call the result MAUR79P.

T Now we can begin the steps for the image ratioing technique. Open the OVERLAY module and select the ratio option
(First/Second) to divide MAUR79P by MAUR77P. Call the result IMGRATIO. Make sure that the Change analysis option is
selected for handling division by zero then click OK. As discussed above, the direct result of the ratioing operation is not
linear nor is it symmetrical about zero. To correct this, open the module TRANSFORM, and select the natural logarithm
(ln(x)) transformation to transform IMGRATIO into a new image called LOGRATIO (see figure).

MAUR79P (OVERLAY /) MAUR77P = IMGRATIO

IMGRATIO (TRANSFORM natural logarithm) = LOGRATIO

8 Using Layer Properties, look at the characteristics of LOGRATIO. What are the minimum and maximum values? Why do we
have negative values? What do those negative values indicate about the change in vegetation between 1977 and 1979?

U Now display a histogram to examine the characteristics of LOGRATIO. Change the minimum for the display to -3.0 and the
maximum to 4.5, and use a class width of 0.05.

9 Why does the histogram have a spike at 0? Does the histogram look reasonably symmetrical? From your examination of the
histogram, within what range do most non-zero values occur? What are the mean and standard deviation?

V
Now reclassify LOGRATIO to create an image as before with class 1 for values less than 3 standard deviations below the
mean, class 0 for those values between -3 and 3 standard deviations, and class 2 for those with values greater than 3
standard deviations above the mean. Call this new image CHGRATIO (see figure). Through Layer Properties, change the
palette to CHANGE and turn off autoscaling.

Criteria

1 = <-3V

EXERCISE 3-8 CHANGE ANALYSIS – PAIRWISE AND MULTIPLE IMAGE COMPARISON 243
 LOGRATIO RECLASS 0 = -3V to +3V = CHGRATIO

2 = >+3V

W
It would appear that very little significant change occurred between 1977 and 1979. Use RECLASS again to create another
image using thresholds based on 2 standard deviations away from the mean. Call this image CHGRAT2 (see figure).
Examine this result with the CHANGE palette.

Criteria

1 = <-2V

LOGRATIO RECLASS 0 = -2V to +2V = CHGRAT2

2 = >+2V

10 Describe the differences you see between the two thresholded images.

When we use thresholds with 3 standard deviations, we can say that 99.73% of the values in the image are due to normal variation
and 0.135% in each "tail" represents significant change (the pixels we see). When we use thresholds with 2 standard deviations, we
can say that 95.45% of the values in the image are due to normal variation while 2.275% in each "tail" represents significant change.

How do we decide what is a significant change in vegetation from one year to the next? From a statistical point of view, this is difficult
to answer with certainty. We would need to investigate other records from those years and perhaps ground-truth the area.

This completes our exploration of pairwise comparison techniques for quantitative data. Try to summarize the differences between
them, then think of how they might be used in your own work.

Change Vector Analysis

Change Vector Analysis can be applied to either pairs of multi-band data or whole time series of single band data. Thus, it is a
technique that bridges both pairwise and multiple comparisons. In this exercise, we will use the red and infrared bands from images

EXERCISE 3-8 CHANGE ANALYSIS – PAIRWISE AND MULTIPLE IMAGE COMPARISON 244
of different dates to examine two components of change detection that are important in change vector analysis -- the magnitude of
change and the direction of change.

This exercise uses two SPOT multi-spectral (XS) images for the Gharb Plain area of Morocco. The Gharb Plain is in the northwestern
corner of Morocco and is crossed by the Sebou River. It is a coastal lowland with deep alluvial deposits and is suitable for intensive
agriculture. During the winter of 1985-86, Morocco received good winter rains after seven years of drought.

We have imagery for two dates in 1986 -- May 10 and June 13. Three bands of imagery are provided for each date, band 1 (green), 2
(red), and 3 (infrared).

X Display the MAY3 and the JUNE3 (infrared) images using equal interval autoscaling and the greyscale color palette.

In the May image, many crops have not reached maturity and show up as dark grey. In contrast, the crops that have reached maturity
show up as light grey (the leaf structure causes high reflectance of infrared energy). In the image JUNE3 (infrared), you will notice
some distinctive changes. Many of the crops that were close to maturity are now mature, while many of the fields that had previously
shown up as mature have now been harvested.

Spatial Registration
Before proceeding with change vector analysis with the Gharb Plain data, we need to introduce the important process of spatial
registration for the purpose of change analysis. Whenever you are comparing two or more images that were collected at different
times or from different sources, spatial registration is a crucial step in the process. Typically, we look at changes over time by
examining the differences in the values of corresponding cells in multiple images. This process only makes sense, however, if the
corresponding pixels of each image describe the same location on the ground. In earlier exercises, the step of registering the images
was already done for us. The two image sets for this exercise have not yet been registered. Since they were taken on separate dates
and thus differ slightly in position and orientation, our first task in this exercise will be to register these images using a process
known as rubber sheet resampling. This technique is covered more thoroughly in the Image Georegistration exercise in the Image
Processing section of the Tutorial. If you are not familiar with the technique, you may wish to complete it before proceeding.

To aid in the process of registration, we will create a new image for each date combining information from all three spectral bands.
These images, called color composites, will allow us to more easily complete the registration task.

Y Use the module COMPOSITE with MAY1, MAY2 and MAY3, assigned to the blue, green and red bands respectively to create
a color composite called MAYCOMP. Choose the linear with saturation endpoints stretch type and the 24-bit composite
with original values and stretched saturation points output type.

Do not omit zeros and enter 3 as the percent to saturate. Do the same with your June images (JUNE1, etc.) to create
JUNECOMP.

This procedure produces what is known as a false color composite. When displayed, the green band is assigned to the blue
component in the resulting image, the red band to the green component, and the infrared to the red component. The result is
therefore not what we would see with our eyes.

Z Arrange the two composite images so they are side by side. Note the differences between the two images, especially in the
pink and red areas.

AA
At the beginning of this exercise, you looked at the infrared bands for these two dates. You were given some hints about
which colors indicated immature, mature and harvested crops. Take a moment to review that information, then compare
the single infrared bands to the color composite images.

EXERCISE 3-8 CHANGE ANALYSIS – PAIRWISE AND MULTIPLE IMAGE COMPARISON 245
11 In the color composite images, what colors seem to represent immature crops, mature crops, and harvested areas?

We will now proceed with the registration. We will leave the June image as it is and will register the May image to it. In order to do
this, we need to precisely (within a single cell if possible) identify several locations on both images for which we can record the
geographic coordinates. Road intersections and other easily visible features are often used. These locations are called control points
and will be used to create a mapping function with which the entire May image will be resampled. The accurate collection of control
point information often requires a fair amount of precision and time (and patience!). The remainder of our change analysis depends
upon a good registration between the two images, so the extra time spent doing this step well is certainly worth the effort.

The spatial registration procedure is somewhat lengthy, but it is a procedure that you will undoubtedly need to undertake if you do
change analysis with your own data. Because of this, we recommend that you take time to complete this section. However, if you do
not wish to complete this part of the exercise, first read through the following steps, then use the Rename option in TerrSet File
Explorer to rename the correspondence (.cor) file GHARBTMP, which was included in your data set, to the new filename of GHARB.
Then rejoin the exercise later below when the GCPs are created.

This correspondence file contains the following data:

1168.557481 9497.598907 1351.812384 9567.990134

239.208719 2362.946385 368.662817 2441.580867

8775.445072 2259.587386 8932.436718 2245.876513

9871.579797 7290.177530 10049.083791 7268.100631

4662.593278 5804.216520 4821.436200 5833.687406

5415.710057 9476.473832 5606.497044 9500.786715

5104.974257 663.387245 5231.711292 688.884075

1352.233630 4291.895646 1503.302825 4365.379702

The first line contains a single whole number indicating the number of control points in the file. Each succeeding line contains two
sets of X and Y coordinates for each control point, the first set from the original referencing system, and the second set from the new
referencing system. Complete details for this format can be found in the TerrSet Help System.

BB
Run the module RESAMPLE. The input file type specifies the type of file to be resampled and can be a raster or a vector file,
or a group of files entered as an RGF. Leave the input file type as raster and specify the input image as MAY2 and the
output image as MAY2RES. We will fill in the output reference parameters later.

The input and output reference files to be specified next refers to the set of images to be used to create the GCPs. For the
input reference image, enter MAYCOMP and for the output reference image, enter JUNECOMP. The images will display in
separate windows.

Before continuing, we need to specify the background value, mapping function and the resampling type.

CC
Enter 0 as the background value.

EXERCISE 3-8 CHANGE ANALYSIS – PAIRWISE AND MULTIPLE IMAGE COMPARISON 246
 A background value is necessary because after fitting the image to a projection, the actual shape of the data may be angled.
In this case, some value needs to be put in as a background value to fill out the grid. The value 0 is a common choice.

The best mapping function to use depends on the amount of warping required to transform the input image into the
output registered image. You should choose the lowest-order function that produces an acceptable result. A minimum
number of control points are required for each of the mapping functions (three for linear, six for quadratic, and 10 for
cubic). Choose the linear mapping function.

The process of resampling is like laying the output image in its correct orientation on top of the input image. Values are
then estimated for each output cell by looking at the corresponding cells underneath it in the input image. One of two
basic logics can be used for the estimation. In the first, the nearest input cell (based on cell center position) is chosen to
determine the value of the output cell. This is called a nearest neighbor rule. In the second, a distance weighted average of
the four nearest input cells is assigned to the output cell. This technique is called bilinear interpolation. Nearest neighbor
resampling should be used when the data values cannot be changed, for example, with categorical data or qualitative data
such as soil types. The bilinear routine is appropriate for quantitative data such as remotely sensed imagery. Since the
data we are resampling is quantitative in character, choose the bilinear resampling type.

We are now ready to digitize control points. It is critical to obtain a good distribution of control points. The points should be spread
evenly throughout the image because the equation that describes the overall spatial fit between the two reference systems will be
developed from these points. If the control points are clustered in one area of the image, the equation will only describe the spatial fit
of that small area, and the rest of the image may not be accurately positioned during the transformation to the new reference system.
A rule of thumb is to try to find points around the edge of the image area. If you are ultimately going to use only a portion of an image,
you may want to concentrate all the points in that area and then window out that area during the resampling process.

As you identify control points, note the total RMS and the individual RMS for each point. The RMS provides an indication of how well
the coordinates listed in the correspondence file fit the mapping function and polynomial equation that were specified in the
RESAMPLE dialog. You should strive to have an RMS less than half the size of a cell in the output image. In this case, an overall RMS
less than 10 meters is acceptable.

You may wish to omit control points with high residuals in order to lower the overall RMS error. RESAMPLE will recalculate the
function based upon the remaining control points. You should try to keep as many points as possible and still attain the acceptable
RMS. Also, ensure that the remaining points are well distributed in the image.

DD
Once you are satisfied with the control points, click on the Output Reference Parameters button. Enter the following
reference system parameters to match the June images. Alternatively, you can select to copy the parameters from any of
the June images.

Number of Columns = 512

Number of Rows = 512

Min. X Coordinate = 0

Max. X Coordinate = 10240

Min. Y Coordinate = 0

Max. Y Coordinate = 10240

Reference System = plane

EXERCISE 3-8 CHANGE ANALYSIS – PAIRWISE AND MULTIPLE IMAGE COMPARISON 247
 Reference Units = m

Unit Distance = 1

After you enter the above information, you are now ready to run RESAMPLE on MAY2. When RESAMPLE finishes, it will
automatically display the resampled image. Note the black areas on the left side and bottom of the image. This will be
explained below.

EE Use RESAMPLE in the same way for the May infrared image (MAY3) using the exact same parameters to create MAY3RES.
You can simply change the input image name and run RESAMPLE again.

In order to fit the May images to the June images, a rubber sheet transformation was applied to the May images. In this case, the May
images had to be rotated slightly in a clockwise direction and shifted slightly to the right to achieve registration with the June images.
This can be confirmed by examining the original and resampled May images. Note the zero-value areas at the left and bottom of the
resampled images. When the images were rotated to match the June orientation, some pixels had no corresponding data in the input
image and were therefore filled with the background value of 0, that was specified in the RESAMPLE dialog. This is illustrated in the
figure below.

We now have areas in the May images that are filled with non-data values. We don’t want to identify change between these
background areas in the May images and the corresponding data values in the June images. There are two ways we might approach
this problem. At this point, we could window out the common area from both the May and June image sets. This makes subsequent
processing easier, but also requires that we exclude some pixels for which we do have data values in both May and June (because a
raster image must be rectangular). The other approach is to continue with the data as they are, but mask out the background areas
whenever necessary. This has the advantage of keeping all relevant data values and only discarding the irregularly shaped mask area.
We will take the first approach.

FF Window into one of the resampled May images and determine the corner row/column numbers for the largest rectangular
area that can be extracted such that it contains no background values.

GG Open WINDOW from the Reformat menu. Indicate that 4 files will be windowed. Then click in each grid line and enter
MAY2RES, MAY3RES, JUNE2 and JUNE3. Enter the output prefix WIN and choose the option to add the prefix to the
filename. Select to specify window coordinates based on row/column positions and enter the following:

Upper Left Column: 9

Upper Left Row: 0

EXERCISE 3-8 CHANGE ANALYSIS – PAIRWISE AND MULTIPLE IMAGE COMPARISON 248
 Lower Right Column: 511

Lower Right Row: 507

Though it is time consuming and often tedious work, spatial registration is an extremely important step in change analyses of all
types. Now that this has been accomplished, we are ready to proceed with the change vector analysis.

Change Vector Extraction

We are now ready to explore the change vector techniques that 1) measure the magnitude of change and 2) the direction of that
change. We will relate the latter to the type of change that occurred, i.e., growth or harvesting. Because some agricultural fields have
experienced growth while others have undergone harvest between these two dates, this data set provides a good illustration of
different types of change that occurred in one location.

To measure the magnitude of change between the two dates, we must use an approach that accommodates the multi-band imagery
we have available. Taking the red and infrared bands for each date, we can imagine that each pixel has a "location" in each of the two
bands (see figure above). The difference between the pixels can then be expressed as the Euclidean distance between them in space.
The formula is:

‫ = ܦ‬ඥ(‫ܦܧܴܣܴܨܰܫ‬ଶ െ ‫ܦܧܴܣܴܨܰܫ‬ଵ )ଶ + (ܴ‫ ܦܧ‬ଶ െ ܴ‫ܦܧ‬ଵ )ଶ

With our images, the distance formula becomes:

‫ = ܦ‬ඥ(ܹ‫ ܧܷܰܬܰܫ‬3 െ ܹ‫ ܻܣܯܰܫ‬3ܴ‫)ܵܧ‬ଶ + (ܹ‫ ܧܷܰܬܰܫ‬2 െ ܹ‫ ܻܣܯܰܫ‬2ܴ‫)ܵܧ‬ଶ

This distance formula could easily be evaluated using Image Calculator. However, the module DECOMP can be used to calculate both
the distance and the direction images, so we will use it.

HH First calculate the simple difference images that will be the X and Y component images submitted to DECOMP using
OVERLAY. Call the Band 3 difference image DIF3 and the Band 2 difference image DIF2. (see figure)

WINJUNE3 (OVERLAY -) WINMAY3RES= DIF3

WINJUNE2 (OVERLAY -) WINMAY2RES = DIF2

II Before running DECOMP, we must first convert the DIF images to real data format. Run CONVERT from the Reformat
Menu. Give DIF3 as the input file, DIF3 as the output file, and choose to create a Real Binary file. Click OK when asked
whether to overwrite the file. Do the same with DIF2.

EXERCISE 3-8 CHANGE ANALYSIS – PAIRWISE AND MULTIPLE IMAGE COMPARISON 249
JJ Open the module DECOMP and choose the option to compose X and Y component images into a force pair. Enter DIF3 as
the input X component image and DIF2 as the input Y component image. Enter DISTANCE as the output magnitude
filename and DIRECT as the output direction filename (see figure). When DECOMP finishes, display DISTANCE with the
quantitative palette.

DIF3 (DECOMP) DIF2 = DISTANCE and DIRECT

12 Where are the areas where the magnitude of change is relatively high? Where are the areas where the magnitude of change is
relatively low?

Now we can focus on examining the direction or type of change that has occurred: where crops have been harvested and where crops
have reached maturity. For each cell, DECOMP has calculated the direction from the location of the May pixel to the location of the
corresponding June pixel. These values are measured as azimuths in degrees clockwise from the positive Y-axis. This is most easily
visualized by thinking of plotting each May pixel in a grid system where the X-axis represents the infrared band, and the Y-axis
represents the red band. The location of the May pixel is the origin. Then the June location is plotted. The angle formed between the
positive Y-axis and a line connecting the May and June locations is the change angle recorded by DECOMP (see figure).

KK Display DIRECT with the quantitative palette.

13 What ranges of angles are most common in the image? (Note: you may want to run HISTO twice, once with the graphic option
and once with the numeric option. You may find it useful to set the min-max values for the histogram to 0 and 360 and the
class width to 1.)

14 What percent of the change angles are found in each 90-degree quadrant (upper right, lower right, lower left, upper left)?
(Hint: use RECLASS to divide the direction image into 90-degree quadrants then use HISTO with the numeric option with the
reclassified image.)

EXERCISE 3-8 CHANGE ANALYSIS – PAIRWISE AND MULTIPLE IMAGE COMPARISON 250
Interpreting the Results
Now we are ready to explore the ways this final image can be interpreted. Angles in the lower right quadrant would seem to indicate
areas that have experienced growth between the two dates (Figure 16). This would generally indicate that values in the infrared
increased and values in the red decreased between May and June. The increase in infrared may be due to a fuller canopy cover, while
the decrease in red indicates that more red light is being absorbed for photosynthesis.

We intended to identify areas of harvest as well as growth in this study area. We would expect harvested areas to have a marked
increase in red reflectance, since the cut vegetation would no longer be absorbing red light for photosynthesis. We would also expect
an increase in infrared reflectance since more of the underlying soil would be exposed. Therefore, we would expect harvested areas
to have change angles in the upper right quadrant.

Most change angles, however, fall in the lower left quadrant, where there was a decrease in infrared as well as a decrease in red
reflectance. The interpretation of this change direction is difficult. The effect of soil moisture on reflectance values for vegetation and
soil has not been addressed in the analysis so far yet may provide an explanation for the absence of change angles in the upper right
quadrant and the prevalence of those in the lower left.

Since the reflectance properties of soil are different from those of vegetation, when the vegetation canopy is not very full, the
reflectance recorded by the sensor is mixed. Dry soil has high reflectance in both the red and infrared, while wet soil absorbs both the
red and infrared wavelengths, resulting in low reflectance values. Harvested areas would allow more of the soil signature to reflect,
so it is possible that the lower left quadrant areas really are harvested, but high soil moisture is depressing both the red and infrared
reflectance. Ground truth data would be necessary to verify this hypothesis, however.

We logically would expect that there are areas where no significant change occurred between the two dates. The question becomes
once again one of thresholding.

EXERCISE 3-8 CHANGE ANALYSIS – PAIRWISE AND MULTIPLE IMAGE COMPARISON 251
15 In this case, if we wanted to establish threshold values in order to differentiate significant from non-significant change, would
we work with the distance image, the direction image, or both? Does a change angle of zero indicate no change?

LL
Display a histogram of DISTANCE, giving new min and max values of 0 and 100 and specifying a class width of 1.

Recall that in earlier exercises, we used the mean plus or minus three standard deviations to define our upper and lower threshold
values. We assumed that values falling outside those thresholds represented significant change. In this case, however, that approach
does not make sense, since the lower part of the distribution is the smallest magnitude of change. For this exercise, it is only the
upper tail of the distribution that represents the largest, and perhaps the most significant, changes.

MM Choose a threshold value beyond which you believe significant change has occurred. (Note that you would normally have
ground truth information available to guide you in setting the threshold value.) Make an image of change/no change areas
using this threshold value with RECLASS and the DISTANCE image. Give the change areas the value 1, and no change areas
the value 0. Use this resulting image to find which change angles are most represented by the larger change distances. Use
OVERLAY to multiply the change/no change image by your DIRECT image.

16 Do the largest change distances correspond to a narrow range of change directions, or are they fairly equally distributed
among all the change directions present?

As you can see, there are several factors that may affect our interpretation and conclusions with respect to the change vector analysis
technique. The development of this exercise has been a part of continuing research in vegetative change detection. The importance of
ground truth information in change analysis must be stressed. By knowing with certainty, the amount and type of change that has
occurred in a few places, we are better able to interpret the changes we see in the images we create for the entire study area.

EXERCISE 3-8 CHANGE ANALYSIS – PAIRWISE AND MULTIPLE IMAGE COMPARISON 252
‬ EXERCISE 3-9
BAYES’ THEOREM AND MAXIMUM
LIKELIHOOD CLASSIFICATION

The next six exercises expand the discussion of classification techniques presented in the Introductory Image Processing exercises.
These exercises will focus on information that can be gleaned from an iterative classification process that provides several layers of
information using soft classifiers. The analyst then reduces that information to a single classified image. If you have not already done
so, read the Classification of Remotely Sensed Imagery section in the TerrSet Manual before continuing with these exercises. The data
for these next set of exercises can be found in the Advanced IP tutorial folder.

We will be working with the same dataset for all six exercises, and results from one exercise may be used for comparison with results
from another. Therefore, if possible, keep all the resulting images from each exercise until the entire set has been completed.

The Maximum Likelihood procedure is unquestionably the most used procedure for classification in remote sensing. The foundation
for this approach is Bayes' Theorem which expresses the relationship between evidence, prior knowledge, and the likelihood that a
specific hypothesis is true. Unfortunately, surprisingly little use is made of the ability to incorporate prior knowledge into the
procedure. Most commonly, analysts make no assumptions about the relative likelihood of finding the land cover classes of interest
before considering the evidence, and thus assume that each class is equally likely. In cases of strong evidence, this will usually do
little harm. However, it is in the context of weak evidence that prior knowledge can make a very important contribution. TerrSet is
unusual in that it offers an especially rich set of options for the inclusion of prior knowledge into the classification process. It offers
the special ability to incorporate prior knowledge in the form of probability images, such that the prior probability of any class is
allowed to vary from one location to the next. As demonstrated in this exercise, this offers a significant improvement in the
classification procedure.

A Display the three images named SPWEST1, SPWEST2 and SPWEST3, each in its own display window using the Greyscale
palette. These are the green, red and near infrared bands from the SPOT-HRV multispectral (XS) sensor for the area of
Westborough, Massachusetts. Form a false color composite of these bands using the COMPOSITE module (from the Display
menu). Enter the bands in the order listed above as the blue, green and red input bands. Call the resulting image
SPWESTFC. Choose a linear stretch with saturation points and create a 24-bit composite that retains the original values.
Give 1% as the amount to saturate on each end. Then display the result.

Westborough is a small rural town that has undergone substantial development in recent years because of its strategic location in
one of the major high-tech development regions in the United States. It is also an area with significant wetland coverage—a land
cover of environmental concern.

B
From Composer, add the vector layer named SPTRAIN using the Qualitative symbol file. Select Map Properties and add a
legend for this layer by choosing the Legend tab. Make the first legend visible and choose SPTRAIN as the layer. You may

EXERCISE 3-9 BAYES’ THEOREM AND MAXIMUM LIKELIHOOD CLASSIFICATION 253

 need to enlarge the map window (by dragging its edge) to view the entire legend. This layer contains a set of training sites
for the following land cover types:

1 Older Residential OLDRES

2 Newer Residential NEWRES
3 Industrial / Commercial IND-COM
4 Roads ROADS
5 Water WATER
6 Agriculture / Pasture AG-PAS
7 Deciduous Forest DECIDUOUS
8 Wetland WETLAND
9 Golf Courses / Grass GOLF-GRASS
10 Coniferous Forest CONIFER
11 Shallow Water SHALLOW

The last column in this list is a set of signature names that will be used in this and the following exercises of this set.

C Use MAKESIG (Image Processing/Signature Development) to create a set of signatures for the training sites in the
SPTRAIN vector file. Indicate that the 3 SPOT bands named SPWEST1, SPWEST2 and SPWEST3 should be used. Choose the
Enter Signature Filenames button and give the signature names in the order listed above.

D MAKESIG automatically creates a signature group file with the same name as the training site file. Signature group files
facilitate use of the classifier dialog boxes. Using TerrSet Explorer, select the signature filter to display files with a “.sgf”
extension. Then in the Files pane verify that the signatures are listed in the signature group file. Then save the file with the
new name called SPOTSIGS. Right-click on the signature group file to rename.

E Run MAXLIKE (Image Processing/Hard Classifiers). In this first classification, we will assume that we have no prior
information on the relative frequency with which different classes will appear. Therefore, choose the option for equal
prior probabilities. Then press the Insert Signature Group button and choose SPOTSIGS. This will fill in the names of all 11
signatures. Leave as 0.0 the minimum likelihood for classification and specify SPMAXLIKE-EQUAL as the output filename.
Press the OK button to run.

F When the classification is completed, display the resulting map using a palette named SPMAXLIKE. Opt to also display the
legend and title. Then compare the result to the false color composite named SPWESTFC.

1 Which classes do you feel the classifier performed best on? Which ones appear to be the worst?

The State of Massachusetts conducts regular land use inventories using aerial photography. The date of the SPOT image used here is
1992. Prior to this, land use assessments had been undertaken in 1978 and 1985. Based on these inventories for the town of
Westborough, CROSSTAB was used to determine the relative frequency with which each land cover class changed to each of the other
classes during the 1978-85 period. These relative frequencies are known as transition probabilities and are the underlying basis for a
Markov Chain prediction of future transitions. If we assume that the underlying driving forces and trajectories of change from 1978
to 1985 have remained stable through 1992, it is possible to estimate the probability with which each land cover class in 1985 might
change to any other class in 1992. These transition probabilities were then applied to the 1985 land cover classes as a base, to yield a
set of probability maps expressing our prior belief that each of the land cover classes will occur in 1992. These images have the
following names:

PRIOR-OLDRES
PRIOR-NEWRES

EXERCISE 3-9 BAYES’ THEOREM AND MAXIMUM LIKELIHOOD CLASSIFICATION 254

 PRIOR-IND-COM
PRIOR-ROADS
PRIOR-WATER
PRIOR-AG-PAS
PRIOR-DECIDUOUS
PRIOR-WETLAND
PRIOR-GOLF-GRASS
PRIOR-CONIFEROUS
PRIOR-SHALLOW

G Display a selection of these prior probability maps using the Default Quantitative palette. Notice that these spatial
definitions of prior probability only extend to the Westborough town boundary. Outside the town boundary, the prior
probability has been expressed as a non-spatial transition probability, much as one would traditionally specify in the use
of the Bayesian Maximum Likelihood Procedure. For example, in the PRIOR-NEWRES image, the area outside the town
boundary has a prior probability of 0.18, which simply represents the likelihood that any area might be expected to be a
newer residential one in 1992. 1 However, the spatially specific prior probabilities range from anywhere up to 0.70,
depending on the existing land cover in 1985.

H Now run MAXLIKE again. Repeat the same steps as were undertaken previously, but this time indicate that you wish to
specify a prior probability image for each signature. Insert the group file SPOTSIGS. Click into the Probability Definition
column of the grid for the first signature. A Pick List button will appear. Click it, then choose the corresponding prior
probability image. For example, the first signature listed should be OLDRES. The probability definition for that line should
be PRIOR-OLDRES. Click into each line in turn and select the prior probability image for that signature. Call the resulting
image SPMAXLIKE-PRIOR. Then click OK to run.

I Display SPMAXLIKE-PRIOR with the SPMAXLIKE palette and indicate that you wish to have a legend. Then add the vector
layer WESTBOUND with the Outline Black symbol file. This layer shows the boundary of the town.

2 Describe those classes in which the most obvious differences have occurred as a result of including the prior probabilities.

J
Use the CROSSTAB module to create a crossclassification image of the differences between SPMAXLIKE-EQUAL and
SPMAXLIKE-PRIOR. Call the crossclassification map EQUAL-PRIOR. Then display EQUAL-PRIOR using the Qualitative
palette, a title and legend. (You may find it useful to create a palette in which the colors for those classes that are the same
between the two images are all white or black. 2 The legend highlight may also be helpful. To highlight a particular
category, hold down the left mouse button on a legend color box.) Add the WESTBOUND vector layer onto your map to
facilitate examination of the effect of the prior probability scheme.

3 Do you notice any other significant differences that were not obvious in question 2 above?

1 This figure is simply the area of the image divided by the area of the newer residential class in 1992.
2 To do so, first open the documentation file for EQUAL-PRIOR with TerrSet Explorer and view its metadata. View its legend categories. Write down the category
numbers of those representing no change (e.g., 1|1, 2|2). There will be 11 of these. Now, open Symbol Workshop. Open the palette file Qual from the TerrSet
program folder's Symbols folder. Choose File/Save As and save it to your Working Folder with a new name, e.g., Equal-Prior. Click on the color boxes for each
of the 11 no-change categories, each time changing that color to be white or black. If there are other palette colors that are white or black that are not on your
list, change their colors to something else. Save the file and use Layer Properties to apply it to the image.

EXERCISE 3-9 BAYES’ THEOREM AND MAXIMUM LIKELIHOOD CLASSIFICATION 255

4 How would you describe the pattern of differences in areas outside the town boundary versus those differences inside?

EXERCISE 3-9 BAYES’ THEOREM AND MAXIMUM LIKELIHOOD CLASSIFICATION 256

‬ EXERCISE 3-10
SEGMENTATION CLASSIFICATION

This tutorial introduces the concept of image segmentation for classification. It builds on the previous exercise using the data for
Westborough, Massachusetts. Set the data path of your Working Folder to Advanced IP in your TerrSet Tutorial data folder.

Classification from segments is a three-step process. The first step is the segmentation of the imagery to the correct level of
generalization. The second step is the development of training sites from the segmentation result. The third and final step is the
classification, based on the the training sites developed in step two as well as a previously classified image.

Segmentation is a process by which pixels are grouped that have homogeneous spectral similarity. The module SEGMENTATION is
used to create an image comprised of segments that have spectral similarity. Across space and over all input bands, a moving window
assesses this similarity, and segments are defined according to a stated similarity threshold. The smaller the threshold, the more
homogeneous the segments. A larger threshold will result in a more heterogeneous and generalized segmentation result. These
segments are then assigned to specific land cover types as we develop training site data and refine the classification process.

A Display the composite image SPWESTFC from the previous exercise.

This is the false color composite image derived from green, red, and near-infrared SPOT imagery, SPWEST1, SPWEST2, and
SPWEST3, respectively. It is from these three bands we will segment and find spectral similarity.

B Open the module SEGMENTATION. Insert the three band files SPWEST1, SPWEST2, and SPWEST3. Specify “0,30,50” in the
Similarity tolerance input box (without quotation marks). Enter the output prefix SPSEG, leave the other defaults and click
OK.

When the SEGMENTATION module has finished,, it will have created three vector files SPSEG_0, SPSEG_30, and SPSEG_50.
We will add these three files, one at a time, to the composite image.

C Display the composite SPWESTFC. Next, in Composer, select Add Layer. (You alternately can hit the V key with the map
window selected.) Add the first vector file SPSEG_50 using the “outline white” palette. Add the other two vector files
SPSEG_30 and SPSEG_0 (in that order) with the same palette. Once all three segment files have been added, you can click
their display on and off from Composer to view the different levels.

Notice that SPSEG_50 contains fewer segments than the other two files, i.e., it is more generalized. The similarity tolerance controls
the level of homogeneity within the segments. Zero is the smallest number that can be used and represents the base watershed, i.e.,
the most homogeneous segments. Numbers greater than zero will result in a more generalized segmentation. We will use SPSEG_30
for the classification process.

EXERCISE 3-10 SEGMENTATION CLASSIFICATION 257

D Close all your map windows and launch the module SEGTRAIN. Select the option to Create a new segment training file.
Enter SPSEG_30 as the segmentation file for sampling and SPWESTFC as the composite background image file. Enter
SEGTRAIN as the output segment training filename. Once the files have been entered into the SEGTRAIN dialog, the display
icon on the center-right of the dialog will become enabled. Click on this icon to display both the segmentation file and
composite image in one map window.

Each segmentation file created contains in its documentation file (.rdc) the names of the bands from which it was created. In
SEGTRAIN, we will interactively select segments that pertain to our classes of interest. When we are finished selecting training
classes and run the module, SEGTRAIN will isolate selected segments as the training classes and then feed these segments to the
module MAKESIG. MAKESIG will create the signatures from the bands from which the segments were derived using the class names
defined in SEGTRAIN. Let’s begin the class selection process.

We are going to select training sites for the following seven classes and select segments that are as homogeneous as possible:

a. deciduous

b. coniferous

c. grass or pasture

d. wetland

e. water

f. residential

g. urban or built

E Make sure that the SPSEG_30 vector file is the highlighted file in Composer. To select segments for training, click the Pick
new sample button on the SEGTRAIN dialog. Then move the cursor to a water body in the map window at approximately
column 230 and row 65. Click once on the segment containing the water body. Notice that it will display the segment ID.
Now double-click the segment to select it.

F
You will notice that the segment ID populates the segment training samples grid in the SEGTRAIN dialog. Enter a class ID
of 5 for this newly selected segment and a class name of Water. Click on the color icon for this selection and choose a basic
blue from the color ramp. When we have finished creating all our training samples, a symbol file with the output filename
will be generated with the colors selected.

G Next, let’s select a segment for deciduous forest. Click the Pick new sample button on the SEGTRAIN dialog and select a
segment in the map window at approximately column 500 and row 180. Double-click to select and enter a Class ID of 1
and a class name of Deciduous.

H We will now select one segment each for the remaining classes. Use the following table as a guide and give the classes
appropriate colors.

Column Row Class ID Class Name

420 505 2 Coniferous
380 100 3 Grass or Pasture
335 275 4 Wetland

EXERCISE 3-10 SEGMENTATION CLASSIFICATION 258

 125 175 6 Residential
407 91 7 Urban or Built

I
J
Select a few more segments per class. Refer to the ones already digitized to find similar segments per class. See the table
below for additional predefined segments.

Column Row Class ID Class Name

105 135 2 Coniferous
115 215 7 Urban or Built
410 75 3 Grass or Pasture
395 110 7 Urban or Built
360 94 7 Urban or Built
355 235 4 Wetland
390 410 1 Deciduous
335 510 1 Deciduous
230 390 6 Residential
110 350 4 Wetland
120 110 6 Residential
360 25 2 Coniferous
185 55 3 Grass or Pasture
480 450 5 Water
470 357 7 Urban or Built
255 140 5 Water

K Once you have selected your segments, click the Create button on the SEGTRAIN dialog.

Now that the training sites are defined, we can begin the classification stage. The module SEGCLASS is used to classify segments
based on an existing reference image. The reference image is a classification image obtained either through a supervised or
unsupervised method. In our case, we will input the training segments we just created to run a maximum likelihood classifier and use
that result as our reference image for the segmentation-based classification.

L
Open the module MAXLIKE, the maximum likelihood classifier. Leave the default to use equal probabilities for each
signature. Click the Insert signature group button and select SEGTRAIN. This is the signature group file created when you
ran the SEGTRAIN module. Call the output image MAX and click OK.

The result is a classification map for our seven classes, based on the training sites developed earlier. We will now refine this result
with the module SEGCLASS running the majority rule classifier.

M Open the module SEGCLASS. Enter the input segmentation file SPSEG_30 and MAX as the pixel classification reference
image. Call the output SEGCLASSMAX and click OK.

N Place the two images SEGCLASSMAX and MAX side by side, set their palettes to SEGTRAIN, and compare the results.

EXERCISE 3-10 SEGMENTATION CLASSIFICATION 259

Notice that the segmentation classification result shows a more generalized map-like result. It may or may not be any more accurate
than the maximum likelihood result, however. Only through ground truthing validation can we make this assessment. You may want
to review your steps and experiment with different segmentation levels for the classification. You also may add additional bands such
as a texture band during the training process.

EXERCISE 3-10 SEGMENTATION CLASSIFICATION 260

‬ EXERCISE 3-11
SOFT CLASSIFIERS I: BAYCLASS

In this exercise, we introduce the concept of a soft classifier. A soft classifier is one that evaluates the degree to which each pixel
belongs to each of a set of land cover classes. Thus, instead of making a definitive (i.e., hard) decision about the class membership of
each pixel, a soft classifier outputs a separate real-number image for each class that expresses set membership on a 0-1 scale. TerrSet
offers a group of soft classifiers, of which the BAYCLASS module is the most approachable.

We will again be using the Westborough SPOT data and the signatures developed in the previous exercise.

A
Run the module BAYCLASS from the Image Processing/Soft Classifiers menu. You will notice that the interface for this
module is almost identical to that of MAXLIKE. Select equal prior probabilities. Indicate that you wish to use the signature
group file named SPOTSIGS. Enter the prefix BAY for the output images. Click OK.

The output from BAYCLASS is in the form of a series of posterior probability maps (BAYOLDRES, BAYNEWRES, BAYIND-COM, etc.).
The values in each represent the evaluated probability that each pixel belongs to that class. BAYCLASS automatically creates two
additional outputs, a raster group file and a classification uncertainty image. The raster group file's name is the same as the prefix
you specified for the output files (i.e., “[Link]” in this case). We will use it to facilitate identifying pixel values across the entire set
of output images. The classification uncertainty image, discussed below, is named BAYCLU.

Because BAYCLASS produces multiple output images, only the classification uncertainty image automatically displays.

B
Display BAYCLU and give it focus. Then from TerrSet Explorer, use the Add Layer command to add BAYDECIDUOUS and
BAYCONIFER to the BAYCLU map window. Then click on the Identify tool icon from the toolbar. Click in the map window
and explore the images (the values of all the images are shown in the Identify box).

1 Compare the BAYDECIDUOUS and BAYCONIFER images. How would you characterize the ability of the classifier to ascertain
whether a pixel belongs to the deciduous class versus the conifer class?

C Notice the distinct forest stand near the top of the BAYDECIDUOUS image that includes the cell at column 324 and row 59.
Use the Zoom Window icon option to window in on this stand. Notice that there is a comparatively greater amount of
uncertainty about many of these pixels compared to other deciduous stands. Use the Identify tool to query several of these
pixels. Activate the Graph option in the Identify box to facilitate your examination.

2 Many of the pixels in this stand have a degree of membership in the deciduous class that is less than 1 (i.e., there is some
uncertainty that the pixel belongs to the deciduous class). For what other class(es) has the classifier indicated some
probability of membership for these pixels?

EXERCISE 3-11 SOFT CLASSIFIERS I: BAYCLASS 261

3 What are the posterior probabilities of all non-zero classes at the cell located at column 326 and row 43. How do you interpret
these data (consider all of these classes in your answer)?

4 Since this is a 20 meter resolution image, each pixel represents 0.04 hectares. For the cell at column 326 and row 43, how
many hectares of deciduous species do you think might exist in this pixel?

5 Examine the cells at column 325, row 43 and column 326, row 43 in the BAYCLU image. What are the uncertainty values at
these locations? What accounts for the difference between them?

6 Examine the cell at column 333, row 37. Notice that the probabilities are fairly evenly spread between three classes. How
many classes were they spread between column 325, row 43? What has been the effect on the uncertainty value? Why?

7 Looking at the BAYCLU uncertainty image as a whole, what classes have the least uncertainty associated with them? Given that
the deciduous category is such a heterogeneous group of species, why do you think the classifier was able to be so conclusive
about this category? (Don't worry too much about your answer here—this is simply a chance to speculate on the reason why.
The reason will be covered in more depth in the next exercise).

D Use EXTRACT (from the IDRISI GIS Analysis/Database Query menu) to extract the average uncertainty associated with
each of the land cover classes in SPMAXLIKE-EQUAL (the Maximum Likelihood classified result created in the first exercise
of this section). Since this image was also created using equal prior probabilities and the non-fuzzy signatures, it
corresponds exactly to the images produced by BAYCLASS. Specify SPMAXLIKE-EQUAL as the feature definition image and
BAYCLU as the image to be analyzed. Then ask for the average summary type and tabular output.

8 What classes have the highest average uncertainties? Can you give a reason why this might be so?

9 Examine the cells in the vicinity of column 408, row 287 on the BAY-CONIFER image. These cells show similar probabilities of
belonging to the wetland and conifer classes. How might you interpret this area? Would you have been able to uncover this if
you had used the MAXLIKE module (compare to the output of SPMAXLIKE-EQUAL)?

EXERCISE 3-11 SOFT CLASSIFIERS I: BAYCLASS 262

‬ EXERCISE 3-12
HARDENERS

In the previous exercise, we produced a series of images expressing the posterior probability of belonging to a set of land cover
classes in the Westborough region. This is a characteristic of all of the soft classifiers. They all defer the issue of making an actual
decision about the land cover class of a pixel. Rather, they simply output the state of one's knowledge about those pixels. We can
force a decision, however, by using a hardener—a module that implements a simple decision logic. The result of using a hardener is a
qualitative land cover image in which each pixel is assigned a single class.

A
Run the HARDEN module. You will find it in the IDRISI Image Processing/Soft Classifiers menu. Select to harden using
posterior probabilities from BAYCLASS. This is the appropriate hardener for use with the output from BAYCLASS. 1 Press
the Insert Layer Group button and choose the group file named BAY (created in the previous exercise). (Do not include
BAYCLU). The number of files indicated should be 11.

Indicate that 4 output levels should be produced. Note that 0 has been entered as the minimum probability value for each
class (by default). 2 Specify BAYMAX as the output file name.

B Display each of the images BAYMAX_1, BAYMAX_2, BAYMAX_3 and BAYMAX_4 from beneath the BAYMAX group file. Use
the SPMAXLIKE palette in each case and specify that a legend should be used. BAYMAX_1 indicates the result of assigning
the class with the maximum probability from the BAYCLASS results. Thus, it will be essentially the same result as that
produced from MAXLIKE (SPMAXLIKE-EQUAL, in this case). 3 BAYMAX_2 indicates the class of the second highest
probability while BAYMAX_3 and BAYMAX_4 indicate the third and fourth highest probabilities respectively.

1 Examine the large stand of deciduous forest in the vicinity of column 583, row 307. Compare the results in BAYMAX1 and
BAYMAX2. How do you interpret those areas where the second highest probability has come out as conifer, wetland or
golf/grass? Examine the probabilities associated with these classes (from the previous exercise) in developing your answer.

1 All the hardener options make calls to MDCHOICE to undertake the analysis. The reason there are separate options for HARDEN is that they have been tailored
to the specific needs of these forms of output.
2 Pixels will be given a value of 0 if they are less than or equal to the value specified for the minimum probability.
3 The result is in fact identical except for the way they may have treated the minimum probability issue. Since HARDEN will assign the value 0 to any pixel with
a probability of belonging to all classes equal to 0, while MAXLIKE will assign an arbitrary choice, the default options may yield a few small differences related
to areas that clearly don't have representation in the classification.

EXERCISE 3-12 HARDENERS 263

2 Notice the striping that is apparent in the third and fourth level images (BAYMAX_3 and BAYMAX_4). Why do you think this
exists? Note also the distinct change that occurs in the vicinity of column 73. This is also related to the same problem as the
striping.

EXERCISE 3-12 HARDENERS 264

‬ EXERCISE 3-13
SOFT CLASSIFIERS II: DEMPSTER-
SHAFER THEORY AND BELCLASS

BELCLASS is the third classifier in the soft classification group and an important counterpart to BAYCLASS. While BAYCLASS is based
on Bayesian probability theory, BELCLASS is based on the variant of Bayesian probability theory known as Dempster-Shafer theory.
If you have not already done so, read the section on BELCLASS in the chapter Classification of Remotely Sensed Imagery in the TerrSet
Manual. You may also wish to read the section on Dempster-Shafer in the Decision Support: Uncertainty Management chapter.

A Run the module named BELCLASS (Image Processing/Soft Classifiers). You will notice that the interface for this module is
quite like that of BAYCLASS. Indicate that you wish to use equal prior probabilities. Then choose the Insert Signature
Group button and select the signature group file named SPOTSIGS that you created in an earlier exercise. Choose the Belief
output option and enter the prefix BEL for the output images. (A raster group file named BEL will also automatically be
created.) Click OK.

B The output from BELCLASS is in the form of a series of Dempster-Shafer belief images (BELOLDRES, BELNEWRES,
BELIND-COM, etc.) and a classification uncertainty image (BELCLU). The latter is autodisplayed. Display the classification
uncertainty image created with BAYCLASS, from the previous exercise, with the Default Quantitative palette and arrange
the two so you can see them both.

1 Describe the difference between BELCLU and the BAYCLU image created in the previous exercise. Given what you have read in
the chapter Classification of Remotely Sensed Imagery, what do you think can account for the fundamental difference between
these images?

C Display several belief images in one map window. Use the Identify tool to examine the values across the images. The values
in each represent the evaluated belief (a form of probability) that each pixel belongs to that class.

D
If not already added, add the image BELDECIDUOUS to the map window along with BAYDECIDUOUS (created with
BAYCLASS in a previous exercise). Look at the large stand of deciduous forest that surrounds the cell at column 215, row
457, and query using the Identify tool.

2 Use the Identify tool with BELDECIDUOUS to examine the beliefs associated with the cells in this stand. What are typical beliefs
for the deciduous class? What are the typical posterior probabilities found in BAYDECIDUOUS for this same area?

EXERCISE 3-13 SOFT CLASSIFIERS II: DEMPSTER-SHAFER THEORY AND BELCLASS 265
3 Notice that the beliefs or probabilities associated with other classes are typically zero or near zero in both cases. How then
does BAYCLASS produce such large probabilities and BELCLASS produce much lower beliefs (remember that they both share
the same underlying mathematical basis)?

4 What do you think might cause the variation in belief in this stand on the BELCLASS image (Hint: consider the issue of the
representativeness of training sites)?

E
Run the module HARDEN to harden these results. Select to calculate beliefs from BELCLASS. Then choose to insert the
layer group BEL. Remove the uncertainty image BELCLU from the set of images to be processed if it is present. Name the
output image BELMAX. Note that you are not asked how many levels to produce. This is because each pixel has a non-zero
belief in only one class. Belief in all other classes is 0. When the result is displayed, change the palette to be SPMAXLIKE.
Then also display the first level image produced the previous exercise with HARDEN, called BAYMAX_1 with that same
palette.

5 How similar are these images? (You may wish to use CROSSTAB with the two images to help you answer this question.)

6 What are the belief and posterior probability values at column 229, row 481? Clearly BAYCLASS (and thus MAXLIKE) has
concluded overwhelmingly that this is an example of deciduous forest. However, given the belief you have determined, is this
reasonable? Is there perhaps another reason other than that given in the answer to question 4 that might account for the
strong difference between these two classifiers? (Hint: BELCLASS implicitly incorporates the concept of an OTHER class in its
calculations—i.e., something other than the classes given in the training sites.)

F Run BELCLASS again and now specify only two signatures: DECIDUOUS and IND-COM. Use the prefix BEL2 for the output.
Then run BAYCLASS and do the same thing using the prefix BAY2.

7 Compare BELDECIDUOUS with BAYDECIDUOUS and BEL2IND-COM with BAY2IND-COM. Given everything you have learned so
far about the difference between these modules, how do you account for the differences/similarities between these two
classifiers in handling this problem? In formulating your answer, compare your results with BAYDECIDUOUS, BELDECIDUOUS,
BAYIND-COM and BELIND-COM.

EXERCISE 3-13 SOFT CLASSIFIERS II: DEMPSTER-SHAFER THEORY AND BELCLASS 266
‬ EXERCISE 3-14
DEMPSTER-SHAFER AND
CLASSIFICATION UNCERTAINTY

In the previous exercise, we saw that BELCLASS provides information on the degree of support for each of a set of land cover classes
independent of the support which is (or is not) provided for the other classes. Dempster-Shafer provides a very rich description of
uncertainty in the classification process, as will be illustrated in this exercise.

A Run BELCLASS with equal prior probabilities. Choose to insert the signature group file named SPOTSIGS that you created
in an earlier exercise. However, this time indicate that you wish to output plausibilities rather than beliefs. Enter the prefix
PLAUS for the output images.

B The output from BELCLASS with this option is in the form of a series of Dempster-Shafer plausibility images
(PLAUSOLDRES, PLAUSNEWRES, PLAUSINDCOM, etc.). The values in each represent the evaluated plausibility, a form of
probability that expresses the highest potential probability that each pixel belongs to that class. Examine these plausibility
images with the Default Quantitative palette. Also examine the PLAUSCLU classification uncertainty image (note that the
PLAUSCLU image is the same as the BELCLU image).

While belief indicates the degree of hard support for a hypothesis, plausibility expresses the degree to which that hypothesis cannot
be disbelieved—i.e., it expresses the degree to which there is a lack of evidence against the hypothesis.

1 Examine PLAUSDECIDUOUS and compare it to BELDECIDUOUS. Overall, how would you describe the plausibility of deciduous
compared to the belief in deciduous? What is the nature of that plausibility in areas in which BELDECIDUOUS is high? Compare
PLAUSDECIDUOUS also to BAYDECIDUOUS from an earlier exercise. How does PLAUSDECIDUOUS compare to BAYDECIDUOUS
in areas where BAYDECIDUOUS is high?

C
Use OVERLAY to subtract BELDECIDUOUS from PLAUSDECIDUOUS (i.e., PLAUSDECIDUOUS - BELDECIDUOUS). Call the
result BELINTDECID. Examine this result using the Default Quantitative palette. This image displays what is called a belief
interval. A belief interval is the difference between the plausibility and the belief for a particular class and expresses a
measure of uncertainty about the state of knowledge about that class.

2 Create similar belief interval images for conifers and wetland. Call the results BELINTCONIF and BELINTWETLAND. How
similar are these images to BELINTDECID?

EXERCISE 3-14 DEMPTER-SHAFER AND CLASSIFICATION UNCERTAINTY 267

D Display the image named PLAUSCLU using the Default Quantitative palette. This is the same image that BELCLASS created
while calculating beliefs, called BELCLU. It is included as an output for use in cases where beliefs have not been output.

3 How similar is PLAUSCLU to the individual uncertainty images BELINTDECID, BELINTCONIF and BELINTWETLAND?

The BELCLU and PLAUSCLU images created by BELCLASS express a very specific form of uncertainty known in Dempster-Shafer
theory as ignorance. Ignorance is different from a belief interval in that a belief interval is category-specific while ignorance applies
to the whole state of knowledge. Ignorance expresses the degree to which the state of knowledge is such that it is unable to
distinguish between the classes. In BELCLASS, we have modified Dempster-Shafer theory to implicitly include an additional class
which we call OTHER, in recognition of the possibility that a pixel belongs to a class for which we have not given a training site. Thus,
ignorance expresses the degree to which we are unable to tell what class the pixel belongs, including the possibility that it is not one
of the classes we are examining.

In the TerrSet implementation of BELCLASS, we also recognize a further aspect of uncertainty that we call ambiguity. Given that
belief expresses the extent of evidence that specifically supports a particular class, ambiguity expresses the degree to which support
is ambiguous because it also supports other classes.

Ambiguity can be calculated as the difference between the belief interval for a specific class and overall ignorance.

E Create an ambiguity image for deciduous by running OVERLAY and subtracting BELCLU (or PLAUSCLU) from
BELINTDECID. Call the result AMBDECID. Notice the degree of ambiguity in the forest stand in the vicinity of the cell at
column 324 and row 59. In the previous exercise on BAYCLASS, we identified this as an area with a significant mixture of
coniferous and deciduous species. The presence of ambiguity gives direct support for the presence of mixtures involving
the class being examined.

F Create a similar ambiguity image for conifers and call it AMBCONIF.

4 How extensive is ambiguity involving either conifers or deciduous?

5 Considering that the total uncertainty of a class (e.g., BELINTDECID) is composed of both ignorance (BELCLU) and ambiguity
(AMBDECID), what is the larger component of uncertainty, ignorance or ambiguity?

As a final note, it is worth considering the issue of sub-pixel classification. The concept of sub-pixel classification assumes that all
uncertainty in the classification of a pixel arises because of the presence of indistinguishable mixtures. However, as has been evident
from this exploration based on Dempster-Shafer theory, ambiguity is not always a major component of uncertainty. Clearly,
ignorance can be a major element. With the range of uncertainty exploration tools provided in TerrSet, however, it is possible to
distinguish between these concepts and focus quite specifically on that aspect which is of greatest concern.

EXERCISE 3-14 DEMPTER-SHAFER AND CLASSIFICATION UNCERTAINTY 268

‬ EXERCISE 3-15
VEGETATION ANALYSIS IN ARID
ENVIRONMENTS

In this exercise, we will explore the use of different vegetation index calculation models available in the VEGINDEX, TASSCAP and
PCA modules to analyze vegetation cover. Before continuing, you may find it useful to read or review the Vegetation Indices section in
the TerrSet Manual. That chapter provides an extensive overview of many vegetation indices, only some of which will be used in this
exercise.

Introduction to Vegetation Indices

Vegetation cover was an early focus of research in natural resources management using space-born satellite images, especially with
the release of the Earth Resources Technology Satellites known as Landsat in 1972. Landsat, SPOT and NOAA data offer time series
images that are widely used to monitor and assess the status of vegetation at the global, regional, national and local levels. Vegetation
indices use various combinations of multi-spectral satellite data to produce a single image representing the amount of vegetation
present, or vegetative vigor. Low index values usually indicate less healthy vegetation while high values indicate more healthy
vegetation. 1 Different indices have been developed to better model the actual amount of vegetation on the ground. The index that is
most appropriate for use in a particular environment can best be determined through calibration with sample measurements of
biomass. In the absence of biomass measurements, these index images can be useful indicators of the relative amount of vegetation
present.

Vegetation has a characteristic spectral response pattern 2 in which visible blue and red energy is absorbed strongly, visible green
light is reflected weakly (hence giving vegetation its green color) and near infrared energy is very strongly reflected. Because of this
characteristic spectral response pattern, many of the vegetation index models use only red and near-infrared imagery bands.

Introduction to the Data and the Study Area

In this exercise, we will assess vegetation cover and its changes in an area of southern Mauritania.

1 Of the 19 vegetation indices produced in the VEGINDEX module, only the RVI and NDVI produce images with high values indicating little vegetation and low
values indicating more vegetation. If you are using a vegetation index model not provided in VEGINDEX, you must determine whether the index values are
proportional or inversely proportional to the amount of vegetation present before you can properly interpret the image.
2 See the Introduction to Remote Sensing and Image Processing chapter in the TerrSet Manual for a discussion of spectral response patterns.

EXERCISE 3-15 VEGETATION ANALYSIS IN ARID ENVIRONMENTS 269

The area covered by the images in this exercise is near the Senegal/Mauritania border and contains part of the Senegal River flood
plain as well as the lower section of the Gorgol River flood plain (partially visible at the upper-left corner of the image). This is a
tributary of the Senegal River. These sections of the two rivers are covered by riverine vegetation dominated by the Acacia nilotica
species, the preferred species for fuelwood and charcoal. Other woody species such as Borassius flabelifer and Iphaene thebaica are
used as building material. Rainfed and flood recessional agriculture and grazing are also practiced in this region.

Once a relatively humid area, persistent rainfall deficits since the late 1960s have left the study area, as well as more and more of the
Sahel, semi-arid. Much vegetation has shifted from savanna to steppe. Relics of the savanna vegetation are only found along river
valleys on clay, clay sand and sandy clay soils, since these retain moisture better than other soils in the area. Increasing pressure
from populations trying to adapt to the continuous drought conditions has been the main cause of vegetation cover degradation in
this environment.

Quantifying the low-density vegetation cover that characterizes arid and semi-arid lands is especially challenging because vegetation
cover is not complete - most pixels contain an average reflectance of vegetation and bare soil. Some of the vegetation index models
we will use have been developed specifically to help account for the effects of background soil reflectance.

The data we will use are Landsat Multi-spectral Scanner (MSS) images. These images were taken on October 10, 1980 and October
12, 1990 by Landsat 4. There are eight images provided in the dataset, four from each year: MAUR80-BAND1, MAUR80-BAND2,
MAUR80-BAND3 and MAUR80-BAND4 for 1980; MAUR90-BAND1, MAUR90-BAND2, MAUR90-BAND3 and MAUR90-BAND4 for
1990. These correspond to MSS bands visible green, visible red, near-infrared and a slightly longer-wavelength near-infrared,
respectively. Since the two scenes were taken at two different dates, they must be registered to one another if we are to do analysis
between them. This task has already been performed using a methodology similar to that described in the exercise on Resample. We
will begin the exercise by producing and comparing several vegetation indices for the 1990 scene, then we will analyze changes
between the two scenes.

Creating Vegetation Index Images

There are three major families of vegetation indices that we will explore: Slope-Based, Distance-Based and Orthogonal
Transformation vegetation indices.

The Slope-Based VI's

The slope-based VI's use the ratio of the reflectance of one band to that of another, usually the red and the near-infrared. The term
slope-based is used because in comparing resulting VI values, we are essentially comparing the slopes of lines passing through the
origin and the pixels as plotted on a graph with the reflectance of one band as the X-axis and the reflectance of the other as the Y-axis.

A Before beginning our exploration of vegetation indices, select User Preferences from the File menu and set the
"Automatically display the output of analytical modules" feature on. We will always display the VI images with a user-
defined palette named NDVI. Go to the Display tab of the User Preferences dialog box and enter NDVI as the Quantitative
Palette. Also, choose to show titles, but do not show legends (this will maximize display space). Click OK to save the
settings and exit User Preferences.

B Use the module VEGINDEX (IDRISI Image Processing/Transformation menu) twice to produce images for two of the slope-
based models: Ratio and NDVI. Use MAUR90-BAND2 as the red band and MAUR90-BAND3 as the near infrared band. Call
the resulting images 90RATIO and 90NDVI. Examine each of the output images. Consult the on-line Help System for details
about the equation used for each index.

EXERCISE 3-15 VEGETATION ANALYSIS IN ARID ENVIRONMENTS 270

1 What similarities and differences do you notice between the two output images? (In answering this question, it may be useful
to look at the pair of images with other quantitative palettes as well, such as Greyscale or Quant.) What is the purpose of
normalizing the Ratio to create NDVI? (You may wish to consult the TerrSet Manual section on Vegetation Indices for help in
answering this question.)

The slope-based VIs are simple linear combinations that use only the reflectance information from the red and infrared bands. In
contrast, the second family of Vegetation Indices that we will explore, the distance-based VIs, uses information about the reflectance
characteristics of the background soil in addition to the red and infrared bands.

The Distance-Based VI's

The reflectance values recorded by the sensor for each pixel constitute an average reflectance of all the cover types in the
instantaneous field of view (i.e., the pixel). When vegetation cover is not complete, which is particularly the case in arid and semi-arid
regions, the average reflectance values are greatly influenced by the background soil type. The distance-based VI's address this
problem of separating information about vegetation from information about soils in remotely sensed data.

The distance-based indices are based on the concept of a soil line and distances from that soil line. A soil line is a linear equation that
describes the relationship between reflectance values in the red and infrared bands for bare soil pixels. This line is produced by
running a simple linear regression between the red and infrared bands on a sample of bare soil pixels. Once that relationship is
known, all unknown pixels in an image that have that same relationship in red and infrared reflectance values are assumed to be bare
soils. Unknown pixels that fall far from the soil line because they have higher reflectance values in the infrared band are assumed to
be vegetation (based on the characteristic spectral response pattern for vegetation where the infrared band reflectance values are
relatively higher than those of the red band). Those that fall far from the soil line because their red reflectances are high are often
assumed to be water (based on the characteristic spectral response pattern for water where the red band reflectance values are
relatively higher than those of the infrared band).

Inputs to the calculation of the distance-based VIs are the red band, the infrared band, the slope of the soil line and intercept of the
soil line. (In addition, some of these VI's also require a scaling factor.)

The first step in calculating the soil line is to identify a sample of bare soil pixels in the image. We will use the 90NDVI image created
earlier to develop a mask image for bare soil. (If better knowledge of the area were available, we could on-screen digitize known bare
soil areas.)

2 If you assume that any pixel having a higher infrared than red reflectance is vegetation and everything else is bare soil, what
threshold value could you use with the 90NDVI image to separate vegetation from bare soils? (Hint: Use the NDVI equation
with some example values to help you answer this question.)

Run RECLASS with 90NDVI to create the image SOILMASK. Assign the new value 1 to bare soil areas and the new value 0
to vegetated areas.

Once the bare soil areas have been identified, the values for those areas in the infrared and red bands are submitted to linear
regression to calculate the soil line. The soil line calculation is not the same, however, for all the distance-based VIs. Some are based
on a regression where the red band is evaluated as the independent variable, and some are based on a regression where the infrared
band is evaluated as the independent variable. Since we will be creating both types of distance-based VI's, you will need to run the
regression twice to determine two soil lines.

EXERCISE 3-15 VEGETATION ANALYSIS IN ARID ENVIRONMENTS 271

C Run REGRESS (from the IDRISI GIS Analysis/Statistics menu) twice, between the MAUR90-BAND2 and MAUR90-BAND3
images, using SOILMASK as the mask image. Write down the slope (a) and intercept (b) values for the case in which the
red band is treated as the independent variable and for the case in which the infrared band is the independent variable. 3

3 What are the slope and intercept when the red band is the independent variable? When the infrared is the independent
variable? What is the coefficient of determination (r2)?

The coefficient of determination is quite high, indicating that the relationship between red and infrared reflectance for these bare soil
pixels is described well by a linear equation.

D Run VEGINDEX three times to produce the distance-based VI's PVI, PVI3, and WDVI. For each VI, refer to the Help System
VEGINDEX Notes section to determine which soil line parameters to use for each VI. Also refer to the Vegetation Indices
section in the TerrSet Manual for details about the equation used for each index.

4 What are the major differences you see in the displays of the three distance-based vegetation index images produced?

5 Is there a noticeable difference between these three images (on average) and the two slope-based images (on average)
produced earlier? In other words, would you be able to separate the five output images into two families based solely on the
resulting images?

The Orthogonal Transformation VIs

The final group of vegetation indices we will explore are the Orthogonal Transformation VI's. With these VI's, four or more bands of
imagery are transformed into a set of new images, one of which describes vegetation. We will explore the use of the Tasseled Cap and
Principal Components transformations for producing vegetation images.

The Tasseled Cap transformation uses a set of four MSS multi-spectral images to produce four new images. 4 The Green Stuff or Green
Vegetation Index (GVI) image represents vegetation. Other images produced represent Soil Brightness Index (SBI), Yellow
Vegetation Index (YVI) and Non-Such Index (NSI). The name of the transformation describes the shape of a plot of pixels in GVI-SBI
space for an image having vegetation in many stages of development. The Tasseled Cap was developed to represent the most
important information from a multi-band agricultural scene in only two images - GVI and SBI.

E Run TASSCAP from the IDRISI Image Processing/Transformation menu. Indicate that you will be using MSS data and enter
the four bands for the 1990 scene. Give 90 as the prefix for the output files. This will produce four images called 90GREEN,
90BRIGHT, 90YELLOW and 90NOSUCH. Display the four images.

3 The equation written at the top of the REGRESS display is in the form y=b+ax, where y=independent variable, b=intercept, a=slope, and x=dependent
variable.
4 The transformation can also be used with six TM images. In this case, three output images are produced, representing greenness, brightness and moistness.

EXERCISE 3-15 VEGETATION ANALYSIS IN ARID ENVIRONMENTS 272

6 Why do you think the areas indicated as having high amounts of vegetation in the green vegetation image show low values in
the soil brightness image?

The Tasseled Cap transformation uses global constants (i.e., the values don't change from scene to scene) to weight the bands being
transformed. Because of this, it may not be appropriate to use it in all environments. Principal components analysis, on the other
hand, is a scene-specific transformation of a set of multi-spectral images into a new set of component images. The component images
are uncorrelated and are ordered according to the amount of variation they explain from the original band set. The first of these
component images typically describes albedo, or brightness, (which includes the background soil) and the second typically describes
variation in vegetative cover.

F Run PCA from the IDRISI Image Processing/Transformation menu. Choose Forward T-mode as the analysis type and the
covariance matrix unstandardized option. Enter 4 as the number of input bands and enter the four 1990 MSS images as
input bands. Enter 4 as the number of components to be extracted. Give 90 as the output file prefix. When the processing is
finished, display the resulting four images, 90_T-MODE_CMP1 through 90_T-MODE_CMP4.

The tabular information produced by PCA indicates that the first component describes nearly 93% of the variance in the original set
of four bands. All the input bands have high and positive loadings for component one. We might then interpret this component as
describing the overall image "brightness." The second component has positive loadings for both infrared bands and negative loadings
for the visible green and red bands. It can be interpreted as an image describing vegetation, independent of the overall scene
brightness. Components three and four describe little of the original variance and appear to represent atmospheric and other noise in
the images.

The equation used for the GVI image of the Tasseled Cap transformation 5 also weights the infrared bands positively and the visible
bands negatively, though the weighting values are somewhat different. It is therefore not surprising to see great similarity between
the second component image and the GVI image produced earlier.

Comparing Vegetation Indices

It is possible to visually compare all the vegetation index images we have produced. Some obviously have better contrast than others.
Some seem to show more variation within the low-value areas. However, without ground-truth information about the status of
vegetation in the area in 1990, we cannot determine which indices are most useful. What we will do is analyze the set of images to
see what different characteristics are illustrated by the various indices.

To do this, we will submit all the VI images we have created in this exercise to a principal components analysis (excluding 90NOSUCH
and 90YELLOW).

G
Run the PCA module. Choose forward t-mode as the analysis type and the correlation matrix standardized option. Indicate
7 as the number of files and enter the names of the seven VI images. Choose to extract 4 components. Give VI as the output
image prefix. The output images will be called VI_T-MODE_CMP1, VI_T-MODE_CMP2, VI_T-MODE_CMP3 and VI_T-
MODE_CMP4. Display these images.

The component images describe the most important "patterns" present in the 7 input vegetation index images. The first component
image shows the pattern which is most common to all the input images. The second component image shows the next most

5 GVI = [(-0.386MSS4)+(-0.562MSS5)+(0.600MSS6)+(0.491MSS7)] In the naming of the image files for this exercise, MAUR90-BAND1 corresponds to MSS4
in the equation, MAUR90-BAND2 to MSS5 and so forth.

EXERCISE 3-15 VEGETATION ANALYSIS IN ARID ENVIRONMENTS 273

important pattern remaining after the first has been removed, and so forth. The statistics produced by PCA include information about
the percent variance explained by each component and the weightings (loadings) of each input image on each component.

7 Compare VI_T-MODE_CMP1 with the input VI images. Which resembles it most? Are the loadings of those input images high
compared to the others for that component?

Research 6 has indicated that in a similar study comparing 25 VI images, the first component described a general vegetation index,
including elements of greenness and soil background. The second component represented those VI's that corrected for soil
background, and the third described soil moisture.

Change Analysis using Vegetation Index Images

We will now undertake an analysis between the two dates of imagery. We will be concerned with identifying areas that have
undergone significant change between 1980 and 1990.

H Display MAUR80-BAND3, the near infrared band of the 1980 image, using the Greyscale palette and autoscaling with Equal
Intervals.

Unfortunately, the data we have for 1980 has significant horizontal "striping" effects due to sensor miscalibration. It is, however, the
best available data for that time and study area, so we will use it. 7

I Choose any one of the vegetation indices you used with the 1990 scene and produce a corresponding image for the 1980
data. If you choose a distance-based VI, you will need to find new soil line parameters for the 1980 data since soil moisture
conditions may be quite different between the two dates and areas of bare soil may have changed.

The most elementary of change analysis techniques is visual comparison.

J Look at the VI image pairs for the two dates and try to determine areas where changes in vegetation are evident. The
striping that is apparent in the 1980 scene is an artifact of the sensor system. Use HISTO with the two vegetation images
and note the average value for the entire image.

8 Does it appear that there is generally more or less vegetation in 1990 than in 1980?

The closest rain-gauge station to this area is the town of Mbout, located outside the image to the East. The station recorded
approximately 200 mm of rain in 1980 and 240 mm of rain in 1990. Since rainfall and vegetation cover are highly correlated, we can
expect to see generally higher vegetation index values in the area for 1990 than for 1980.

There are many quantitative methods we can use to analyze change between images. Here we will explore only one, simple
differencing. For a more complete treatment of change analysis techniques, see the Time Series/Change Analysis section in the
TerrSet Manual. You may use the data from this exercise to explore on your own many of the techniques presented in that chapter.

6 Thiam, Amadou, 1997. Geographic Information and Remote Sensing Systems Methods for Assessing and Monitoring Land Degradation in the Sahel Region:
The Case of Southern Mauritania. PhD Dissertation, Clark University, Worcester, Massachusetts.
7 You may wish to try to mitigate the striping by using Fourier analysis with these 1980 images. Use the forward transform, filter out the horizontal elements,
then use the backward transformation. See the section on Fourier Analysis in the TerrSet Manual for more information.

EXERCISE 3-15 VEGETATION ANALYSIS IN ARID ENVIRONMENTS 274

With simple differencing, we merely subtract one image from the other, then analyze the result. The critical issue then becomes one
of setting an appropriate threshold for the difference image beyond which we consider real change, as opposed to ephemeral
variation, to have occurred. Ground truth information would normally be used to identify these thresholds.

K Use OVERLAY to subtract your 1990 image from your 1980 image. Call the resulting image 1980-1990. Use HISTO with
1980-1990 and change the class width to be small in relation to the range of values in 1980-1990. (The class width will
differ depending on the VI you chose to use. Make sure there are at least 100 "bins" or divisions in the histogram.) Note the
distribution of values, as well as the mean and standard deviation.

In the absence of ground truth information to guide our selection of a suitable change/no-change threshold, we will use the standard
deviation. We will consider that only those pixels lying beyond two standard deviations from the mean in either the positive or
negative direction constitute real change and those lying within two standard deviations represent normal variation. In a normal
distribution, 90% of the values fall within two standard deviations. By setting this as our threshold, therefore, we are identifying the
outlying 10% of pixels as our significant change areas.

L Use RECLASS with 1980-1990 and the mean and standard deviation values you found above to create a new image,
CHANGE, in which areas showing a significantly negative change in vegetation from 1980 to 1990 have the value 1, areas
with normal variation have the value 2 and areas with significantly positive change from 1980 to 1990 have the value 3.

9 What is the distribution of positive and negative change areas in the study area? (Try to disregard change that is due to the
sensor miscalibration in the 1980 imagery.)

10 Optional: Repeat steps A through L for several other vegetation indices and compare the results. How much does the choice of
vegetation index influence the final assessment of change?

EXERCISE 3-15 VEGETATION ANALYSIS IN ARID ENVIRONMENTS 275