MATLAB Examples

Using ToR_LandSat8.m

ToR_LandSat8 has the following general form:

output=ToR_LandSat8(Data,operationList,bandList);

where:

  • Data is the LandSat 8 Data that must be loaded using loadLandSat8, which can be downloaded from here,
  • operationList tells the function what to do, more on this later,
  • bandList is the optional field, that can be ommitted and tells the function what bands need to be processed.

Note that operationList is not case sensitive. Currently the following operations are supported:

  • 'TOARad': Top Of Atmosphere Radiance
  • 'TOARef': Top Of Atmosphere Reflectance
  • 'SatBT': at Satellite's Brightness Temperature
  • 'DOS1': At surface reflectance calculated using DOS1 algorithm.
  • 'DOS2': At surface reflectance calculated using DOS2 algorithm.
  • 'COST': At surface reflectance calculated using COST algorithm.

Notes:

  • DOS2 and COST are the same algorithm and they produce the same results.
  • If single operation is needed just pass it as character. But if multiple operation is requested pass it as cell array, eg. use {'TOARad,'SatBT'} to perform both operations.

This example, shows how to load a sample LandSat8 images to calculate brightness temperature, TOA reflectance and surface reflectance using both DOS1 and DOS2/COST.

It also shows Normalized Digital Vegetation Index (NDVI); and how it changes if it is calculated using Digital Numbers (DN), TOA reflectance, and surface reflectance (both DOS1 and DOS2/COST).

Note that NDVI should never be calculated using DN. The equation for all vegetation indices are for surface reflectance (not even TOA Reflectances).

Also note that DOS1 and DOST2/COST algorithm are very simple atmospheric correction algorithm. In many cases they differ considerably from the reflectance that is obtained using a proper atmospheric correction algorithm, such as the one obtained by using MODTRAN. So, make sure to examine the outputs to see if they make sense. Also note that these algorithms were developed for previous versions of LandSat imageries. I have not yet seen any paper evaluating these algorithms on LandSat8. This should serve as a warning for using these two methods.

Contents

Loading LandSat8 data

The data file must be loaded using loadLandSat8.m (download it from MATHWORKS FileExchange). loadLandSat8 requires the path to the metadata file. Note that individual landSat image files (TIF files) must be in the same location as the metadata.

metaFilename='./LC80400372014241LGN00_MTL.txt';

We want to calculate NDVI; hence, band 4 and 5 are needed. Also we want to calculate at satellite brightness temperature; so band 10 and band 11 are also needed. We want to show the true color composite of this scene, so band 2,3, and 4 are also needed. By setting bandList we can load only those LandSat8 bands that we need. If bandList is omitted all bands are loaded.

bandList=[2,3,4,5,10,11];

% Loading the LandSat8 data
Data=loadLandSat8(metaFilename,bandList);

Converting to "At Satellite Brightness Temperature".

To calculate at Satellite Brightness Temperature set operationList to 'SatBT' when calling ToR_LandSat8.m. Note that operationList is not case sensitive. so you can set it to 'satbt' and same results is obtained.

SatBT=ToR_LandSat8(Data,'SatBT')
Calculating At Satellite Brightness Temperature ...

SatBT = 

    SatBT_bandList: [10 11]
             SatBT: {2x1 cell}

SatBT is a structure with two fields.

  • SatBT_bandList: that shows the band number that was processed.
  • SatBT: which is a cell array and each element is a processed band. To check what band number it was, you can check the SatBT_bandList field. The first number in SatBT_bandList field tells the original band number for the first element in the cell array. In general the i-th number in the SatBT_bandList field tells the band number of the i-th element in the cell array.

Displaying the at Satellite Brightness Temperature for band 10 (Left), band 11 (right).

figure
imshow([SatBT.SatBT{1} SatBT.SatBT{2}],[],'InitialMagnification',8);
colormap(hot);
colorbar;
title('At Satellite''s Brightness Temperature [K]; Band 10 (left), Band 11(Right)');
drawnow();

Displaying the difference between these two temperature:

figure
imshow(SatBT.SatBT{1}-SatBT.SatBT{2},[],'InitialMagnification',8);
colormap(jet);
colorbar;
title('The difference between the two Brightness temperature')
drawnow();

Calculating "Top Of the Atmosphere Reflectance (TOA)"

To calculate TOA reflectance set operationList to 'TOARef'. To process only certain bands you can set bandList accordingly. In this example we will set bandList=[4,5,10,11]; However, note that TOA reflectance can not be calculated for band 10 and 11. Even if band 10 and 11 are listed in bandList , they are going to be ignored for this calculation.

TOARef=ToR_LandSat8(Data,'TOARef',[4,5,10,11])
Calculating Top Of Atmosphere (TOA) Reflectance ...

TOARef = 

    TOARef_bandList: [4 5]
             TOARef: {2x1 cell}

TOARef is a structure with two fields.

  • TOARef_bandList: that shows the band number that was processed.
  • TOARef: which is a cell array and each element is a processed band. To check what band number it was, you can check the TOARef_bandList field. The first number in TOARef_bandList field tells the original band number for the first element in the cell array. In general the i-th number in the TOARef_bandList field tells the band number of the i-th element in the cell array.

Calculating "At Surface Reflectance Using DOS1 Method"

To calculate surface reflectance using DOS1 method you must set operationList to 'DOS1'. Note that here the bandList is not provided. This forces the function to process all loaded bands elligible for DOS1 operation.

Reflectance_DOS1=ToR_LandSat8(Data,'DOS1')
Performing Dark-Object Subtraction 1 (DOS1) ...
The following Bands are not loaded:
     1     6     7     8     9

These bands are ignored.

Reflectance_DOS1 = 

    SurfaceRef_DOS1_bandList: [2 3 4 5]
             SurfaceRef_DOS1: {4x1 cell}

Reflectance_DOS1 is a structure with two fields.

  • SurfaceRef_DOS1_bandList: that shows the band number that was processed.
  • SurfaceRef_DOS1: which is a cell array and each element is a processed band. To check what band number it was, you can check the SurfaceRef_DOS1_bandList field. The first number in SurfaceRef_DOS1_bandList field tells the original band number for the first element in the cell array. In general the i-th number in the SurfaceRef_DOS1_bandList field tells the band number of the i-th element in the cell array.

Calculating "At Surface Reflectance Using DOS2 (COST) Method"

To calculate surface reflectance using DOS2 or COST method you must set operationList to 'DOS2' or 'COST'. Note that these are the same algorithm. Note that here the bandList is not provided. This forces the function to process all loaded bands elligible for this operation.

Reflectance_COST=ToR_LandSat8(Data,'COST')
% Alternative you could use
% Reflectance_COST=ToR_LandSat8(Data,'DOS2')
Performing Dark-Object Subtraction 2 (DOS2) (also known as COST) ...
The following Bands are not loaded:
     1     6     7     8     9

These bands are ignored.

Reflectance_COST = 

    SurfaceRef_DOS2_COST_bandList: [2 3 4 5]
             SurfaceRef_DOS2_COST: {4x1 cell}

Reflectance_COST is a structure with two fields.

  • SurfaceRef_DOS2_COST_bandList: that shows the band number that was processed.
  • SurfaceRef_DOS2_COST: which is a cell array and each element is a processed band. To check what band number it was, you can check the SurfaceRef_DOS2_COST field. The first number in SurfaceRef_DOS2_COST field tells the original band number for the first element in the cell array. In general the i-th number in the SurfaceRef_DOS2_COST field tells the band number of the i-th element in the cell array.

Calculating NDVI Using DN, TOA, Surface Reflectance

Now calculating Normalized Digital Vegetation Index or NDVI. The equation for this index is:

$NDVI=\frac{\rho_{NIR}-\rho_{Red}}{\rho_{NIR}+\rho_{Red}}$

note that $\rho$ is the surface reflectance. However, in this example we calculate NDVI using Digital Numbers (DN), TOA Reflectance, and Surface Reflectance obtained using DOS1 and DOS2/COST algorithm. Note that NDVI should never be calculated using DN. This is quite a common mistake; It is widely not noticed, because what you get as NDVI using DN resembles very much the DN that you get using reflectance. In fact sometimes it is visually even more appealing.

NDVI_Fun=@(NIR,Red) (NIR-Red)./(NIR+Red);
NDVI_DN=NDVI_Fun(double(Data.Band{5}), ...
                 double(Data.Band{4}));
NDVI_TOA=NDVI_Fun(TOARef.TOARef{TOARef.TOARef_bandList==5}, ...
                  TOARef.TOARef{TOARef.TOARef_bandList==4});
NDVI_DOS1=NDVI_Fun(Reflectance_DOS1.SurfaceRef_DOS1{Reflectance_DOS1.SurfaceRef_DOS1_bandList==5}, ...
                   Reflectance_DOS1.SurfaceRef_DOS1{Reflectance_DOS1.SurfaceRef_DOS1_bandList==4});
NDVI_COST=NDVI_Fun(Reflectance_COST.SurfaceRef_DOS2_COST{Reflectance_COST.SurfaceRef_DOS2_COST_bandList==5}, ...
                   Reflectance_COST.SurfaceRef_DOS2_COST{Reflectance_COST.SurfaceRef_DOS2_COST_bandList==4});

Note the use of:

TOARef.TOARef{TOARef.TOARef_bandList==4}

Using the approach shown above gaurantess that we are loading TOA reflectance corresponding to 4th band of LandSat8. In this example we processed only band 4 and 5 for TOA reflectance. Therefore, Band 4 is the first element in the cell array, and band 5 is the second element. However, if later you change your code and include band 2, 3, 4, and 5 for TOA then TOA reflectance corresponding to band 4 of LandSat8 is the 3rd element of the cell array. Therefore, retrieveing the bands as shown above is recommended relative to retrieving them as TOARef.TOARef{1}. Although both command retrieve the same data in this example; however, the second form is quite confusing for retrieving band 4 (while using index 1) and if later you change your code it might fail or retrieve a wrong band.

Displaying All NDVI's next to each other

figure
imshow([NDVI_DN  , NDVI_TOA; ...
        NDVI_DOS1, NDVI_COST], [-1,1],'InitialMagnification',4);
colormap(winter);
colorbar;
title('NDVI calculated using: DN (upper left), TOA (upper right), DOS1 (Lower lLeft), DOS2/COST (lower right)')
drawnow();

showing the difference between NDVI calculated using DOS1 and DOS2

figure,
imshow(NDVI_DOS1-NDVI_COST,[],'InitialMagnification',8);
colormap(jet);
colorbar;
title('Difference between NDVI calculation using DOS1 and DOS2/COST');
drawnow();

Focusing on Balboa park, San Diego Displaying All NDVI's next to each other

figure
imshow([NDVI_DN(5470:5643,3773:3943)  , NDVI_TOA(5470:5643,3773:3943); ...
        NDVI_DOS1(5470:5643,3773:3943), NDVI_COST(5470:5643,3773:3943)], [-1,1], ...
        'InitialMagnification',200);
colormap(winter);
colorbar;
title('NDVI - Balboa Park San Diego using: DN (upper left), TOA (upper right), DOS1 (Lower lLeft), DOS2/COST (lower right)')
drawnow();

As you can see NDVI is considerably increased over vegetations, when reflectances are used, relative to when DN is used.

Showing difference between NDVI calculated using DN or surface reflectance.

figure,
imshow( [NDVI_DN-NDVI_DOS1, NDVI_DN-NDVI_COST],[],'InitialMagnification',8);
colormap(jet);
colorbar;
title('Difference between NDVI calculated using DN and DOS1 (right) or DOS2/COST (left)');
drawnow();

Showing RGB image

finally combining Red, Green, and Blue band to create a true color composite of the LandSat8 scene. Note the use of imadjust() and stretchlim() to properly adjust the colors. LandSat8 RGB images usually turn out to be very dark and they need some histogram adjustment for better visualizations.

% Combining the Red, Green, and Blue band in one variable.
Irgb_COST=cat(3,Reflectance_COST.SurfaceRef_DOS2_COST{Reflectance_COST.SurfaceRef_DOS2_COST_bandList==4}, ...
                Reflectance_COST.SurfaceRef_DOS2_COST{Reflectance_COST.SurfaceRef_DOS2_COST_bandList==3}, ...
                Reflectance_COST.SurfaceRef_DOS2_COST{Reflectance_COST.SurfaceRef_DOS2_COST_bandList==2});
% Setting the no data to zero so it shows black.
Irgb_COST(isnan(Irgb_COST))=0.0;
% Adjusting the histogram and showing the results.
imshow(imadjust(Irgb_COST,stretchlim(Irgb_COST)),'InitialMagnification',8);
drawnow();

Performing multiple operation in one call

In order to perform multiple operation in one call use the function as follows:

output=ToR_LandSat8(Data,{'SatBT','TOARad'});

In this case output would be:

output =

    SatBT_bandList: [10 11]
             SatBT: {2x1 cell}
   TOARad_bandList: [2 3 4 5 10 11]
            TOARad: {6x1 cell}

References

  • GRASS GIS, http://grass.osgeo.org/grass65/manuals/i.landsat.toar.html
  • http://semiautomaticclassificationmanual.readthedocs.org/en/latest/Landsat_conversion.html
  • Chander G., B.L. Markham and D.L. Helder, 2009: Remote Sensing of Environment, vol. 113
  • Chander G.H. and B. Markham, 2003.: IEEE Transactions On Geoscience And Remote Sensing, vol. 41, no. 11.
  • Chavez P.S., jr. 1996. Image-based atmospheric corrections - Revisited and Improved. Photogrammetric Engineering and Remote Sensing 62(9): 1025-1036.
  • Huang et al: At-Satellite Reflectance, 2002: A First Order Normalization Of Landsat 7 ETM+ Images.
  • R. Irish: Landsat 7. Science Data Users Handbook. February 17, 2007; 15 May 2011.
  • Markham B.L. and J.L. Barker, 1986: Landsat MSS and TM Post-Calibration Dynamic Ranges, Exoatmospheric Reflectances and At-Satellite Temperatures. EOSAT Landsat Technical Notes, No. 1.
  • Moran M.S., R.D. Jackson, P.N. Slater and P.M. Teillet, 1992: Remote Sensing of Environment, vol. 41.
  • Song et al, 2001: Classification and Change Detection Using Landsat TM Data, When and How to Correct Atmospheric Effects? Remote Sensing of Environment, vol. 75.
  • Sobrino, J.; Jimenez-Munoz, J. C. & Paolini, L. 2004. Land surface temperature retrieval from LANDSAT TM 5 Remote Sensing of Environment, Elsevier, 90, 434-440