Quantcast

Documentation Center

  • Trial Software
  • Product Updates

Finding Vegetation in a Multispectral Image

This example shows how differences between the visible red and near-infrared (NIR) bands of a LANDSAT image can be used to identify areas containing significant vegetation. It uses a LANDSAT Thematic Mapper image covering part of Paris, France, made available courtesy of Space Imaging, LLC. Seven spectral bands are stored in one file in the Erdas LAN format.

Variations in the reflectivity of surface materials across different spectral bands provide a fundamental mechanism for understanding features in remotely-sensed multispectral imagery.

Step 1: Import CIR Bands from a BIL Image File

The LAN file, paris.lan, contains a 7-band 512-by-512 Landsat image. A 128-byte header is followed by the pixel values, which are band interleaved by line (BIL) in order of increasing band number. They are stored as unsigned 8-bit integers, in little-endian byte order.

The first step is to read bands 4, 3, and 2 from the LAN file using the MATLAB® function multibandread.

Thematic Mapper 4, 3, and 2 bands cover the near infrared (NIR), the visible red, and the visible green parts of the electromagnetic spectrum. When they are mapped to the red, green, and blue planes, respectively, of an RGB image the result is a standard color-infrared (CIR) composite. The final input argument to multibandread specifies which bands to read, and in which order, so that you can construct a composite in a single step.

CIR = multibandread('paris.lan', [512, 512, 7], 'uint8=>uint8',...
                    128, 'bil', 'ieee-le', {'Band','Direct',[4 3 2]});

Variable CIR is a 512-by-512-by-3 array of class uint8. It is an RGB image, but with false colors. When the image is displayed, red signifies the NIR band, green signifies the visible red band, and blue signifies the visible green band. In the CIR image, water features are very dark (the Seine River) and green vegetation appears red (parks and shade trees). Overall, the contrast is low and the colors are subtle.

figure
imshow(CIR)
title('CIR Composite (Un-enhanced)')
text(size(CIR,2), size(CIR,1) + 15,...
  'Image courtesy of Space Imaging, LLC',...
  'FontSize', 7, 'HorizontalAlignment', 'right')

Step 2: Enhance the CIR Composite with a Decorrelation Stretch

It's helpful, before analyzing the CIR composite, to enhance it for more effective visual display. Because of the subtle color differences in the original composite, a decorrelation stretch is suitable. You can use the function decorrstretch, which enhances color separation across correlated channels, and also specify an optional linear contrast stretch to saturate the brightest and darkest one percent of the pixels in each band.

decorrCIR = decorrstretch(CIR, 'Tol', 0.01);
figure
imshow(decorrCIR)
title('CIR Composite with Decorrelation Stretch')

The surface features have become much more obvious and the image is more colorful. This is because the spectral differences across the scene have been exaggerated and the contrast has been increased.

Much of the image appearance is due to the fact that healthy, chlorophyll-rich vegetation has a high reflectance in the near infrared. Because the NIR band is mapped to the red channel in our composite, any area with a high vegetation density appears red in the display. A noticeable example is the area of bright red on the left edge, a large park (the Bois de Boulogne) located west of central Paris within a bend of the Seine River.

By analyzing differences between the NIR and red bands, you can quantify this contrast in spectral content between vegetated areas and other surfaces such as pavement, bare soil, buildings, or water.

Step 3: Construct an NIR-Red Spectral Scatter Plot

A scatter plot is a natural place to start when comparing the NIR band (displayed as red) with the visible red band (displayed as green). It's convenient to extract these bands from the original CIR composite into individual variables. (We return to the original bands, because the decorrelation-stretched image is appropriate for visual display but not for spectral analysis.) It's also helpful to convert from class uint8 to class single so that the same variables can be used in the NDVI computation below, as well as in the scatter plot.

NIR = im2single(CIR(:,:,1));
red = im2single(CIR(:,:,2));

Viewing the two bands together as grayscale images, you can see how different they look.

figure
imshow(red)
title('Visible Red Band')
figure
imshow(NIR)
title('Near Infrared Band')

With one simple call to the plot command in MATLAB, you can create a scatter plot displaying one point per pixel (as a blue cross, in this case), with its X-coordinate determined by its value in the red band and its Y-coordinate by the value its value in the NIR band.

figure
plot(red, NIR, '+b')
set(gca, 'XLim', [0 1], 'XTick', 0:0.2:1,...
         'YLim', [0 1], 'YTick', 0:0.2:1);
axis square
xlabel('red level')
ylabel('NIR level')
title('NIR vs. Red Scatter Plot')

The appearance of the scatter plot of the Paris scene is characteristic of a temperate urban area with trees in summer foliage. There's a set of pixels near the diagonal for which the NIR and red values are nearly equal. This "gray edge" includes features such as road surfaces and many rooftops. Above and to the left is another set of pixels for which the NIR value is often well above the red value. This zone encompasses essentially all of the green vegetation.

Step 4: Compute Vegetation Index via MATLAB® Array Arithmetic

Observe from the scatter plot that taking the ratio of the NIR level to red level would be one way to locate pixels containing dense vegetation. However, the result would be noisy for dark pixels with small values in both bands. Also notice that the difference between the NIR and red should be larger for greater chlorophyll density. The Normalized Difference Vegetation Index (NDVI) is motivated by this second observation. It takes the (NIR - red) difference and normalizes it to help balance out the effects of uneven illumination such as the shadows of clouds or hills. In other words, on a pixel-by-pixel basis subtract the value of the red band from the value of the NIR band and divide by their sum.

ndvi = (NIR - red) ./ (NIR + red);

Notice how the array-arithmetic operators in MATLAB make it possible to compute an entire NDVI image in one simple command. Recall that variables red and NIR have class single. This choice uses less storage than class double but unlike an integer class also allows the resulting ratio to assume a smooth gradation of values.

Variable ndvi is a 2-D array of class single with a theoretical maximum range of [-1 1]. You can specify these theoretical limits when displaying ndvi as a grayscale image.

figure
imshow(ndvi,'DisplayRange',[-1 1])
title('Normalized Difference Vegetation Index')

The Seine River appears very dark in the NDVI image. The large light area near the left edge of the image is the park (Bois de Boulogne) noted earlier.

Step 5: Locate Vegetation -- Threshold the NDVI Image

In order to identify pixels most likely to contain significant vegetation, you can apply a simple threshold to the NDVI image.

threshold = 0.4;
q = (ndvi > threshold);

The percentage of pixels selected is thus

100 * numel(NIR(q(:))) / numel(NIR)
ans =

    5.2204

or about 5 percent.

The park and other smaller areas of vegetation appear white by default when displaying the logical (binary) image q.

figure
imshow(q)
title('NDVI with Threshold Applied')

Step 6: Link Spectral and Spatial Content

To link the spectral and spatial content, you can locate above-threshold pixels on the NIR-red scatter plot, re-drawing the scatter plot with the above-threshold pixels in a contrasting color (green) and then re-displaying the threshold NDVI image using the same blue-green color scheme. As expected, the pixels having an NDVI value above the threshold appear to the upper left of the rest and correspond to the redder pixels in the CIR composite displays.

% Create a figure with a 1-by-2 aspect ratio
h = figure;
p = get(h,'Position');
set(h,'Position',[p(1,1:3),p(3)/2])
subplot(1,2,1)
% Create the scatter plot
plot(red, NIR, '+b')
hold on
plot(red(q(:)), NIR(q(:)), 'g+')
set(gca, 'XLim', [0 1], 'YLim', [0 1])
axis square
xlabel('red level')
ylabel('NIR level')
title('NIR vs. Red Scatter Plot')
% Display the thresholded NDVI
subplot(1,2,2)
imshow(q)
set(h,'Colormap',[0 0 1; 0 1 0])
title('NDVI with Threshold Applied')

See also decorrstretch, im2single, MultispectralImageEnhancementExample, multibandread.

Was this topic helpful?