Blog GIS & Territories

We have previously seen how to correct Landsat images for TOA (top-of-atmosphere)reflectance. Here is a diagram that will help us understand the following:

Therefore the reflectance we have calculated is the% of light reflected to the totality of the incident visible light. But as we can see in the diagram,the satellite sensor measures two things at the same time: the light reflected by our targets on the earth’s surface, plus the light scattered by the particles suspended in the air.

We can push the atmospheric correction of the satellite images to remove the light due to scattering.

As we have already said, this is of no interest unless if you work on separate images in time. The diffusion percentage being the same for an image, therefore it is time lost to make this adjustment.

The instructions to perform this correction on Landsat images is available at  

But if you are beginner, you will find difficult to understand the procedure. We’ll discuss the procedure by using the QGis tools.

DOS Method (Dark Object Subtraction)

The method we will discuss is the most frequently used.
The subtraction of black objects is a simple empirical atmospheric adjustment method for satellite imaging, which assumes that the reflectance of dark objects includes an important component of atmospheric scattering.

The concept is simple: if you have a completely black object, and therefore does not reflect light at all, this object should have a value of reflectance equal to zero, will have a nonzero value at the satellite sensor. This value will be the reflectance due to the particles present in the atmosphere.

The subtraction of dark objects researches in each band the darkest pixel value. The broadcast is deleted by subtracting this value of each pixel of the band.

First thing: the adjustment we are going to perform does not apply only
to the visible light and near infrared. Therefore we will work for Landsat 8
images, on the bands 2 (Blue), 3 (green), 4 (Red) and 5 (near infra-red).

Second thing to know: the principle used is to calculate the scattering for band 4 (red) and to use a correspondence formula to determine the adjustment for the other three bands.

Third thing: we will work with bands 2 , 3, 4 and 5 in TOA reflectance but, also, in gross radiance (DN) for band 4 (red). In the user manual cited above this does not appear at first glance. Let’s see the procedure step by step:

1- the TOA reflectance for the four bands with is calculated with the following formula

TOA band = ( (DN * 0.00002)-0.1) / sin ( solar elevation)

2- by analysing the gross radiance (DN) values of the red band , we determine a radiance due to atmospheric diffusion ( DNred )

3- we transform this radiance to TOA reflectance

CORred = ( ( DNred* 0.00002) -0.1) / sin ( solar elevation  )

4- using an abacus we  enter the red band diffusion and obtain the diffusion for the bands 2,3 and 5 ( CORblue ,CORgreen , CORir )

5- we subtract the TOA reflectance values to obtain the reflectance values of the four bands adjusted for the atmospheric effect.

Red = TOA ( red) – CORred

Rblue = TOA ( blue) – CORblue

Rgreen= TOA ( green) – CORgreen

Rir = TOA ( ir ) – CORir

Let’s discuss in detail steps 2 and 4.

Calculation of the red band diffusion

There are several methods for determining the adjustment to be made for the red band. The recommended method is the so-called Bin 5 – 0.008. Here you will find a detailed description of the method.

In practice, with QGis, we will use the band histogram to determine the threshold value of the atmospheric scattering. 1- Firstly, you must run the Raster Layer Statistics tool on the Band 4 (red) data. 

Once executed, you will have a result line in the Results Viewer panel 

Double-click on this line to open the results in your browser:  

Note the minimum value, here 3916.

We will look for the value to be assigned to the correction of the atmospheric diffusion in the image histogram.

Click on the image layer -> Properties -> Histogram  

Enter the minimum value found in the statistics (here 3916) in the Min box, then click on the Max box below to move the minimum value of the dotted line on the histogram.  

You can use the line format to find the base of the histogram, but for this example, the bar format will be used and it is recommended. If the histogram is in line format, click Preferences / Actions and clear the Mark as lines check box.

A histogram in bar format will appear. The QGIS histogram does not represent / show all the values as ArcGIS histogram, but it does display the format that may be appropriate to establish a scatter value that represents the base of the histogram (if increasing relatively abruptly the frequencies).  The following steps show how the QGIS histogram can be used to establish dispersion.

Hover over the histogram to view a magnifying glass. Drag the magnifying glass to the bottom end of the histogram to zoom the histogram base. If you zoom too much or if you want to do it again for any reason, click Preferences/ Actions and Recalculate Histogram. The histogram is zoomed in below.

You will see that the bar with the lowest value is still quite superior to the low image value you previously entered in the Min (3619) area; the QGIS histogram does not show many values ​​in a statistical tail (Landsat 8 usually has long statistical tails). The base of the histogram is the area used to establish the dispersion. By zooming in the three first bars of the histogram:  

In this example, the values ​​between roughly 5570 and 5595 may approach to the base of the histogram – these are very low values ​​and are not completely on the statistical tail (this process may be somewhat subjective). For the purposes of this document, the DN 5569 will be used as the base of the histogram. The histogram will be zoomed one more time for this example.

The rule to use is to take as dispersion value the minimum value of the first bar exceeding a frequency of 5 (hence the name Bin 5) and having no frequencies lower than 5 to its right. In the following example:

The first two bars exceed the value 5, but since the third one is lower,
the min value of the fourth bar will be retained.

In our example, the radiance value to retain for the atmospheric adjustment of the red band will be 5569.

We can continue to the third step,

3- we transform this radiance into TOA reflectance

CORred = ((DNred * 0.00002) -0.1)
/ sin (solar elevation)

  CORred = ((5569 *
0.00002) -0.1) /0.42631886   = 0.026694

This value corresponds to the Bin 5. Since the method is Bin 5 – 0.008,
we remove 0.008 from this value: 0.018694

It is this value that must be subtracted from the red TOA value previously calculated to obtain the adjusted values ​​of the atmospheric diffusion.

Adjustment of other bands

From the correction value of the red band, we obtain the values to be applied to each of the other bands:

Go to

You will find a window where you can enter the adjustment value for the red

Re-enter your value (here 0.018647), you will automatically obtain the values ​​of the other corrections that will display  

We thus obtain the corrections of

• 0.06309 for the Blue Band (2)
• 0.03442 for the Green band (3)
• 0.00626 for the near infrared band (5)

To obtain the corrected bands we use the raster calculator by subtracting to each TOA adjusted band the correction constant that we found.

Si cet article vous a intéressé et que vous pensez qu'il pourrait bénéficier à d'autres personnes, n'hésitez pas à le partager sur vos réseaux sociaux en utilisant les boutons ci-dessous. Votre partage est apprécié !

Atilio Francois