Calculating top of atmosphere (TOA) reflectance with QGIS

admin

The data in Landsat-8 images is represented by Digital Numbers (DN). For the calculation of environmental indices such as the Normalized Difference Vegetation Index, reflectance values are often used instead of DN. The purpose of the following exercise is to calculate the top of atmosphere reflectance (TOA reflectance) with a solar angle correction for band 4 and band 5 in QGIS 3+. We choose band 4 (red) and band 5 (near infrared or NIR) in this example. Band 4 and band 5 are needed for the calculation of NDVI.


Table of contents for this post:


Import data

We use a Landsat-8 image of Portugal (path: 204, row: 32). Landsat satellite data is available on https://earthexplorer.usgs.gov/.

First of all, we need to import the data of band 4 and band 5.

In the toolbar:
Layer >> Add Layer >> Add Raster Layer…

Click […] to navigate to your folder and select both files, band 4 and band 5. Click [Add] to import your Landsat-8 image file into the map canvas.

Notice that for band 4 the range of values goes from 0 to 19,375 and for band 5 from 0 to 22,081? These high DN numbers are due to the fact that Landsat-8 data is delivered in 16-bit unsigned integer format which allows values from 0 to 65,536 (216 or 65,536 unique values), while data for Landsat 1 – 7 are delivered in 8-bit unsigned integer format which ranges from 0 to 255 (28 or 256 unique values).

Set 0 values to nodata

0 is assigned to the outer borders of the satellite image. We do not want to calculate the border, therefore we assign nodata to the 0 values. We can use the raster conversion function translate in QGIS to do so. Behind this functionality in QGIS is the GDAL translate function (gdal_translate), which allows us to convert raster formats. In our case it is used for replacing the 0 values for nodata.

In the toolbar:
Raster >> Conversion >> Translate (Convert Format)…

Simply set assign a specified nodata value to output bands (optional) to 0. This is what will transform our 0 values to nodata.

Save your output file by clicking […] after Converted >> Save to File… Navigate to your folder of choice and give a name accordingly.

In our example your window for conversion should look like this:

Click [Run] to import the band 4 image with nodata assigned to the 0 values to the map canvas.
Repeat the process for band 5.

Conversion to TOA reflectance

For the conversion to TOA reflectance, we use the guidelines as indicated by the USGS .


DN’s can be converted to TOA reflectence using the following formula:

\rho_\lambda´= M_\rho * Q_{cal} * A_\rho


Where:

\rho_\lambda´ = TOA planetary reflectance, without correction for solar angle (this will be corrected later);
M_\rho = Band-specific multiplicative rescaling factor from the MTL file (REFLECTANCE_ADD_BAND_x, where x is the band number);
Q_{cal} = Quantized and calibrated standard product pixel values (DN);
A_\rho = Band-specific additive rescaling factor from the MTL file (REFLECTANCE_ADD_BAND_x, where x is the band number).


Retrieving values for the parameters

The MTL file is a .txt that is packed together with your data download. It is a text file containing the metadata of the image. For the image used here, the MTL filename is: LC08_L1TP_204032_20191230_20200111_01_T1_MTL.txt.

Here you find the values for the parameters M_\rho and A_\rho for each band:

In our specific case:
M_\rho = reflectance_mult_band_4 = 2.0000E-05
A_\rho = reflectance_add_band_4 = -0.100000

Be careful, these values do not necessarily have to be the same for every image.


Calculating TOA reflectance in the raster calculator

To calculate \rho_\lambda´ in QGIS, open te raster calculator: Raster >> Raster Calculator…

Insert in the raster calculator the formula with the values you found for M_\rho and A_\rho:

2.00000E-5 * "B4\_nodata@1" - 0.100000

Save your output file by clicking […] after Output layer >> Navigate to your folder of choice and give a name accordingly. Here we call the output file TOA_B4_notcorrected.

Click [OK] to process the image and calculate the reflectance values. Repeat for band 5.

The result should generate values between 0 and 1:


Correct TOA reflectance for solar angle

We still need to correct the TOA reflectance values for the solar angle. We are using local sun elevation angle to correct the TOA reflectance:

\rho_\lambda = \frac{\rho_\lambda´}{sin(\theta_{SE})}


Where:

\rho_\lambda = planetary TOA reflectance
\theta_{SE} = local sun elevation angle from the MTL file (SUN_ELEVATION).


Retrieving the value for sun elevation

In the MTL file you find the value for \theta_{SE}:


\theta_{SE} = SUN ELEVATION = 23.71021788


Before we calculate this formula in the raster calculator, we have to be aware that the MTL file gives us the local sun elevation angle in degrees and QGIS calculates in radians. To convert degree to radian we have to multiply the degrees with \frac{\pi}{180}.

With the conversion, the formula to correct for solar angle looks like this:

\rho_\lambda = \frac{\rho_\lambda´}{sin(\theta_{SE} * \frac{\pi}{180})}


Calculating corrected TOA reflectance in the raster calculator

Open te raster calculator: Raster >> Raster Calculator…

Insert in the raster calculator the formula with the value you found for \theta_{SE}:


\frac{"TOA\_B4\_notcorrected@1"}{sin(23.71021788 * \frac{3.141593}{180})}


Note that for \pi we use the value: 3.141593.

Save your output file by clicking […] after Output layer >> Navigate to your folder of choice and give a name accordingly. Here we call the output file TOA_B4_corrected.

Click [OK] to process the image and calculate the reflectance values corrected for solar angle. Repeat for band 5.

The result should generate values between 0 and 1:


In the QGIS map canvas the output looks as follows:

Close-up of solar angle corrected TOA reflectance for band 5 (NIR), location Figueira da Foz


Some thoughts…

You have succesfully converted DN’s into solar angle corrected TOA reflectance values. Though, keep in mind that this is the procedure for when you analyse one image. When you intend to compare, for example, NDVI values over multiple images, it is advised to also apply an atmospheric correction to the bands used. Thus, you calculate surface reflectance and not TOA reflectance. An interesting read on this topic is A Survival Guide to Landsat Preprocessing.

Young, N.E., Anderson, R.S., Chignell, S.M., Vorster, A.G., Lawrence, R. and Evangelista, P.H., 2017. A survival guide to Landsat preprocessing. Ecology, 98(4), pp.920-932.

Image courtesy of the U.S. Geological Survey
Landsat-8, sensor: OLI_TIRS, path 204, row 32, date acquired: 2019-12-30
Landsat product ID: LC08_L1TP_204032_20191230_20200111_01_T1