Module 7
Caribbean Regional Node

To demonstrate how to carry out the radiometric correction of digital imagery, using two Landsat Thematic Mapper images obtained at different seasons under different atmospheric conditions as examples.
 
 

This lesson relates to material covered in Chapter 7 of the Remote Sensing Handbook for Tropical Coastal Management and readers are recommended to consult this for further details of the techniques involved. The lesson is rather a specialist one designed to guide practitioners in radiometric and atmospheric correction; it is advanced and is quite hard work (be warned!).

Atmospheric correction will be carried out on two Landsat Thematic Mapper images of the Caicos Bank obtained at different seasons and under somewhat different atmospheric conditions. The first Landsat TM image was acquired in November 1990 whilst the second image (simulated) is for the rather different atmospheric conditions and sun elevation of June 1990. At the time of the November overpass horizontal visibility was estimated at 35 km whilst for the June one it was only 20 km. The sun elevation angle for the winter overpass was 39- but that for the summer overpass was 58-. The DN values recorded for the same areas of the Earth's surface thus differ considerably between the two images.

Familiarity with Bilko for Windows 2.0 is required to carry out this lesson. In particular, you will need experience of using Formula documents to carry out mathematical manipulations of images. Some calculations need to be performed independently; these can either be carried out on a spreadsheet such as Excel or using a calculator.

The first image was acquired by Landsat 5 TM on 22nd November 1990 at 14.55 hours Universal Time (expressed as a decimal time and thus equivalent to 14:33 GMT). The Turks & Caicos are on GMT - 5 hours so the overpass would have been at 09:33 local time. You are provided with bands 1 (blue), 2 (green) and 3 (red) of this image as the files TMNOVDN1.GIF, TMNOVDN2.GIF, and TMNOVDN3.GIF. These images are of DN values but have been geometrically corrected. The second Landsat-5 TM image has been simulated for the rather different atmospheric conditions and sun elevation of 22 June 1990 at 14.55 hours Universal Time by the reverse of the process you are learning to carry out in this lesson (i.e. surface reflectance values have been converted to DN values at the sensor). You are provided with bands 1 (blue), 2 (green) and 3 (red) of this image as the files TMJUNDN1.GIF, TMJUNDN2.GIF and TMJUNDN3.GIF. These images are also of DN values and have been geometrically corrected so that pixels can be compared between seasons. The centre of each scene is at 21.68- N and 72.29- W.

Digital sensors record the intensity of electromagnetic radiation (ER) from each spot viewed on the Earth's surface as a digital number (DN) for each spectral band. The exact range of DN that a sensor utilises depends on its radiometric resolution. For example, a sensor such as Landsat MSS measures radiation on a 0-63 DN scale whilst Landsat TM measures it on a 0-255 scale. Although the DN values recorded by a sensor are proportional to upwelling ER (radiance), the true units are W m^-2 ster^-1 m m^-1 (Box 3.1).

The majority of image processing has been based on raw DN values in which actual spectral radiances are not of interest (e.g. when classifying a single satellite image). However, there are problems with this approach. The spectral signature of a habitat (say seagrass) is not transferable if measured in digital numbers. The values are image specific - i.e. they are dependent on the viewing geometry of the satellite at the moment the image was taken, the location of the sun, specific weather conditions, and so on. It is generally far more useful to convert the DN values to spectral units.

This has two great advantages:

1) a spectral signature with meaningful units can be compared from one image to another. This would be required where the area of study is larger than a single scene or if monitoring change at a single site where several scenes taken over a period of years are being compared.

2) there is growing recognition that remote sensing could make effective use of "spectral libraries" - i.e. libraries of spectral signatures containing lists of habitats and their reflectance (see Box 3.1).

While spectral radiances can be obtained from the sensor calibration, several factors still complicate the quality of remotely sensed information. The spectral radiances obtained from the calibration only account for the spectral radiance measured at the satellite sensor. By the time ER is recorded by a satellite or airborne sensor, it has already passed through the Earth's atmosphere twice (sun to target and target to sensor).

During this passage (Figure 3.1), the radiation is affected by two processes: absorption which reduces its intensity and scattering which alters its direction. Absorption occurs when electromagnetic radiation interacts with gases such as water vapour, carbon dioxide and ozone. Scattering results from interactions between ER and both gas molecules and airborne particulate matter (aerosols). These molecules and particles range in size from the raindrop (>100 m m) to the microscopic (<1 m m). Scattering will redirect incident electromagnetic radiation and deflect reflected ER from its path (Figure 3.1).

The unit of electromagnetic radiation is W m^-2 ster^-1 m m^-1. That is, the rate of transfer of energy (Watt, W) recorded at a sensor, per square metre on the ground, for one steradian (three dimensional angle from a point on Earth's surface to the sensor), per unit wavelength being measured. This measure is referred to as the spectral radiance. Prior to the launch of a sensor, the relationship between measured spectral radiance and DN is determined. This is known as the sensor calibration. It is worth clarifying terminology at this point. The term radiance refers to any radiation leaving the Earth (i.e. upwelling, toward the sensor). A different term, irradiance, is used to describe downwelling radiation reaching the Earth from the sun (Figure 3.1). The ratio of upwelling to downwelling radiation is known as reflectance. Reflectance does not have units and is measured on a scale from 0 to 1 (or 0-100%).

Absorption and scattering create an overall effect of "haziness" which reduces the contrast in the image. Scattering also creates the "adjacency effect" in which the radiance recorded for a given pixel partly incorporates the scattered radiance from neighbouring pixels.

In order to make a meaningful measure of radiance at the Earth's surface, the atmospheric interferences must be removed from the data. This process is called "atmospheric correction". The entire process of radiometric correction involves three steps (Figure 3.2).

The spectral radiance of features on the ground are usually converted to reflectance. This is because spectral radiance will depend on the degree of illumination of the object (i.e. the irradiance). Thus spectral radiances will depend on such factors as time of day, season, latitude, etc. Since reflectance represents the ratio of radiance to irradiance, it provides a standardised measure which is directly comparable between images.

A considerable amount of additional information is needed to allow you to carry out the radiometric correction of an image. Much of this is contained in header files which come with the imagery. Two tables of information relating to the Landsat TM imagery are included here whilst other information is introduced in the lesson as needed. Table 3.1 is extracted from the November 1990 Landsat TM image header, whilst Table 3.2 contains some satellite specific information you will need.

Question: 3.1. Why do you think the Landsat TM3 DN values for the deep water and for the coral reef and seagrass areas are the same? Answers

This is a fairly straightforward process which requires information on the gain and bias of the sensor in each band (Figure 3.3). The transformation is based on a calibration curve of DN to radiance which has been calculated by the operators of the satellite system. The calibration is carried out before the sensor is launched and the accuracy declines as the sensitivity of the sensor changes over time. Periodically attempts are made to re-calibrate the sensor.

Figure 3.3. Calibration of 8-bit satellite data. Gain represents the gradient of the calibration. Bias defines the spectral radiance of the sensor for a DN of zero.

The calibration is given by the following expression for at satellite spectral radiance, Ll :

The method for calculating gain and bias varies according to when the imagery was processed (at least, this is the case for imagery obtained from EOSAT). Gains and biases for each band l are calculated from the lower (Lminl ) and upper (Lmaxl ) limits of the post-calibration spectral radiance range (Table 3.1 and Figure 3.3). Since the imagery was processed after October 1st 1991 we can use Equation 3.2 to calculate the gains and biases for each waveband which are needed to solve Equation 3.1 and convert the DN values to at satellite spectral radiances.

Where data have been processed after October 1st 1991 (as in this case), the values of Lmax and Lmin for Landsat TM can be obtained from the header file which accompanies the data. The header file is in ASCII format and the gains/biases are stored in fields 21-33 in band sequence. The gains and biases have been extracted from the header file for you and are displayed in Table 3.1. These values are given as in-band radiance values (mW cm^-2 ster^-1) and need to be converted to spectral radiances across each band (mW cm^-2 ster^-1m m^-1). This is done by dividing each value by the spectral band width. Band widths for the TM sensors carried on Landsats 4 and 5 are listed in Table 3.2. The Landsat TM images used here were taken from the Landsat-5 satellite.

The apparent reflectance, which for satellite images is termed exoatmospheric reflectance, r , relates the measured radiance, L (which is what the formulae above will output), to the solar irradiance incident at the top of the atmosphere and is expressed as a decimal fraction between 0 and 1:

A detailed discussion of the methods available for atmospheric correction is available in Kaufman (1989). Atmospheric correction techniques can be broadly split into three groups:

1. Removal of path radiance (e.g. dark pixel subtraction which will be carried out in Lesson 4),
2. Radiance-reflectance conversion, and
3. Atmospheric modelling (e.g. 5S radiative transfer code, which will be used here).

Using the 5S Radiative Transfer Code: The model predicts the apparent (exoatmospheric) reflectance at the top of the atmosphere using information about the surface reflectance and atmospheric conditions (i.e. it works in the opposite way that one might expect). Since the true apparent reflectance has been calculated from the sensor calibration and exoatmospheric irradiance (above), the model can be inverted to predict the true surface reflectance (i.e. the desired output). In practice, some of the model outputs are used to create inversion coefficients which may then be applied to the image file. We cannot run the model programme here but will introduce you to the inputs needed and provide you with the outputs from a Unix version of the model which we have run for this Landsat TM data. There are three stages in using the 5S code.

Stage 1 - Run the 5S code for each band in the imagery

The inputs of the model are summarised in Table 3.6. Note that it is possible to input either a general model for the type of atmospheric conditions or, if known, specific values for atmospheric properties at the time the image was taken. In the event that no atmospheric information is available, the only parameter that needs to be estimated is the horizontal visibility in kilometres (meteorological range) which for our November image was estimated at 35 km (a value appropriate for the humid tropics in clear weather), but for the June image was poorer at only 20 km.

Stage 2 - Calculate the inversion coefficients and spherical albedo from the 5S code output

The 5S code provides a complex output but only some of the information is required by the user. The following example of the output for Landsat TM1 (Box 3.2) highlights the important information in bold.

  Equation 3.4

Equation 3.5

The coefficients A_I and B_I and the exoatmospheric reflectance data derived in Step 2 can then be combined using Formula documents to create give new images Y for each waveband using Equation 3.6. [Note the minus sign in front of the Reflectance value in Equation 3.5, which means that B_I will always be negative].

where p are the values stored in your new connected images (@1, @2 and @3).

You now have the information needed to carry out Stages 2 and 3 of the atmospheric correction which will convert your exoatmospheric reflectances to at surface reflectances. It is best to carry this out as two stages with intermediate formulae to create the new images Y for each waveband (Equation 3.6) and final formulae, incorporating these, to carry out Equation 3.7.

Once you have set up the formulae to create the new images (Y) for each waveband using Equation 3.6, you can take the appropriate spherical albedo (S) values from Table 3.8 and use the following equation to obtain surface reflectance r_s on a scale of 0-1:

At the start of the lesson you compared the DN values of the November and June images at certain row and column coordinates. Now you will apply a pre-prepared Formula document to the June images to correct them and then compare the surface reflectance values at five coordinates in the corrected June images with those in the corrected November images.

1. The red band penetrates poorly so that by about 5 m depth even bright sand is not likely to reflect any more light back than deep water. The seagrass area has a very low albedo (is very dark) and being thick vegetation will have a very low reflectance in the red waveband anyway. Thus this area reflects no more red light than deepwater. The coral reef area will also have low reflectance in the red band due to the presence of zooxanthellae in the coral tissues (with similar pigments to other photosynthetic organisms) and is also probably too deep to reflect more red light than deep water areas.

The predecence of operators ( * / + - ) means that only those brackets which have been included are needed. Thus division and multiplication always precede addition and subtraction. The brackets are needed to make sure that the addition of the Lmaxl /254 and Lminl /255 to calculate the gain are carried out before the DN values are multiplied by the resultant gain.

2. The Julian Day corresponding to 22 November 1990 is 326.
3. The square of the Earth-Sun distance in astronomical units (d?) on 22 November 1990 is 0.975522 AU (Astronomical Units)? [1 -0.01674 cos(0.9856x322)]?
4. The sun zenith angle in degrees (SZ) is 90--39- = 51- (where 39- is the sun elevation angle). In radians this is 51 x p /180 = 0.89012.
5. The values of ESUN for Landsat-5 TM1, TM2 and TM3 are: 195.7, 182.9 and 155.7 respectively.
6. If you needed to convert a SPOT XS1 solar exoatmospheric spectral irradiance (ESUN) of 1855 W m^-2 ster^-1 m m^-1 to units of mW cm^-2 ster^-1 m m^-1, you need to consider change in two of the four units. W are changing to mW (x 1000) and units of power per m? are changing to units of power per cm?. There are 100 cm in a metre and thus 100 x 100 = 10,000 cm? in one square metre, so we need also to divide by 10,000. So the net result is x 1000/10,000 which is equivalent to divide by 10. So 10 W m^-2 ster^-1 m m^-1 = 1 mW cm^-2 ster^-1 m m^-1. Thus 1855 W m^-2 ster^-1 m m^-1 = 185.5 mW cm^-2 ster^-1 m m^-1.

7. If the output image were an 8-bit integer image all pixel values would either be 0 or 1 (i.e. they would be rounded to the nearest integer).
8. Because of the extra haze in the atmosphere in June more of the incoming sunlight is reflected back at the sensor making the whole image brighter with higher DN values in each band at all pixels. Once the atmospheric conditions have been corrected for then the images appear remarkably similar with pixel values almost identical. In reality correction would not be this good but this is what it aims to do.
9. Because raw DN values recorded from given sites on the Earth's surface will vary with atmospheric conditions, time of year, etc. Thus unless radiometric and atmospheric correction have been carried out any comparison of pixel values at given locations is fairly meaningless. Any differences could be due as much to changes in sun angle, the distance of the Earth from the sun or atmospheric conditions as to changes in the habitat.
10. The average absolute difference between the November and June TM band 2 DN values expressed as a percentage of the average of the ten DN values is 28.28%.
11. The average absolute difference between the November and June TM band 2 surface reflectance values at the five coordinates expressed as a percentage of the average of the ten surface reflectance values is 0.94%.

# Start of Bilko for Windows formula document to radiometrically correct Landsat-5 TM bands

# 1-3 collected over the Turks and Caicos on 22 November 1990.

#

# Formula document 1.

# =================

#

# Step 1. Converting DN to at satellite spectral radiance (L) using formulae of the type:

#

# Lmin + (Lmax/254 - Lmin/255) * @n ;

#

# Input values

# ==========

# Lmin: TM1 = -0.116, TM2 = -0.183, TM3 = -0.159

# Lmax: TM1 = 15.996; TM2 = 31.776; TM3 = 24.394

#

const Lmin1 = -0.116 ;

const Lmin2 = -0.183 ;

const Lmin3 = -0.159 ;

const Lmax1 = 15.996 ;

const Lmax2 = 31.776 ;

const Lmax3 = 24.394 ;

# Intermediate formulae for L for each TM band:

#

# Lmin1 + (Lmax1/254 - Lmin1/255)*@1;

# Lmin2 + (Lmax2/254 - Lmin2/255)*@2;

# Lmin3 + (Lmax3/254 - Lmin3/255)*@3;

#

# Step 2. Converting at satellite spectral radiance (L) to exoatmospheric reflectance

#

# Input values

# ==========

# pi = 3.141593

# d? = 0.975522 JD = 326 for 22/11/90 image. (call dsquared)

# SZ = 90-39 = 51- = 0.89012 radians

# ESUN: TM1 = 195.7, TM2 = 182.9, TM3 = 155.7

const pi =3.141593 ;

const dsquared = 0.97552 ;

const SZ = 0.89012 ;

const ESUN1 = 195.7 ;

const ESUN2 = 182.9 ;

const ESUN3 = 155.7 ;

#

# Let at satellite spectral radiance = L (see intermediate formulae above)

#

# Converting L to exoatmospheric reflectance (on scale 0-1) with formulae of the type:

#

# pi * L * dsquared / (ESUN * cos(SZ)) ;

#

pi * (Lmin1 + (Lmax1/254 - Lmin1/255)*@1) * dsquared / (ESUN1 * cos(SZ)) ;

pi * (Lmin2 + (Lmax2/254 - Lmin2/255)*@2) * dsquared / (ESUN2 * cos(SZ)) ;

pi * (Lmin3 + (Lmax3/254 - Lmin3/255)*@3) * dsquared / (ESUN3 * cos(SZ)) ;
 
 

# Start of Bilko for Windows formula document to atmospherically correct Landsat-5 TM bands

# 1-3 collected over the Turks and Caicos on 22 November 1990.

#

# Formula document 2.

# =================

# Stage 2 of atmospheric correction using 5S radiative transfer model outputs

#

# Input values

# ==========

# AI = 1 / (Global gas transmittance * Total scattering transmittance)

# TM1 = 1.3056, TM2 = 1.2769, TM3 = 1.1987

#

# BI = - Reflectance / Total scattering transmittance

# TM1 = -0.0992, TM2 = -0.0515, TM3 = -0.0301

#

const AI1 = 1.3056 ;

const AI2 = 1.2769 ;

const AI3 = 1.1987 ;

const BI1 = -0.0992 ;

const BI2 = -0.0515 ;

const BI3 = -0.0301 ;

# Let exoatmospheric reflectance = @n (i.e. images output by first formula document)

#

# Converting exoatmospheric reflectance (scale 0-1) to intermediate image Y with formulae of the type:

# AI * @n + BI;

#

# Intermediate formulae for Y:

#

# AI1 * @1 + BI1;

# AI2 * @2 + BI2;

# AI3 * @3 + BI3;

#

# Stage 3 of atmospheric correction using 5S radiative transfer model outputs

#

# Input values

# ==========

# S = Spherical albedo: TM1 = 0.156, TM2 = 0.108, TM3 = 0.079

#

const S1 = 0.156 ;

const S2 = 0.108 ;

const S3 = 0.079 ;

# Let intermediate image = Y (see intermediate formulae above)

#

# Converting Y to surface reflectance (on scale 0-1) with formulae of the type:

#

# Y / (1 + S * Y) ;

#

(AI1 * @1 + BI1) / (1 + S1 * (AI1 * @1 + BI1) ) ;

(AI2 * @2 + BI2) / (1 + S2 * (AI2 * @2 + BI2) ) ;

(AI3 * @3 + BI3) / (1 + S3 * (AI3 * @3 + BI3) ) ;
 
 

Back to main page

Back to Module 7

USDOC | NOAA | NESDIS | CoastWatch