In this third post of the Landsat 8 series we will implement, step by step, the procedures explained in the Landsat web site for converting the digital numbers in the satellite images to geophysical measurements of spectral radiance, as measured at the sensor and without atmospheric correction, called ‘top-of-atmosphere’ or simply TOA.

The images that we download as TIF files are coded into 16-bit unsigned integer images. These are referred to as digital numbers. For some simple purposes, as for example, measuring the extent of features identifiable by eye, these digital numbers are perfectly sufficient. However, when we need to assess the geophysical properties of the terrain or make studies of land cover change over time, we might find convenient to convert these raw digital numbers to their original values of radiance, mostly to be able to compute reflectance, and compare surfaces on a solid geophysical basis.

Converting to radiance is extremely easy because the images are coded using a linear relationship between radiance measurement and digital numbers. For all bands, we retrieve these values by calculating:

L_{\gamma} = M_{L} Qcal + A_{L}

where L_{\gamma} is the spectral, or band-wise, TOA radiance, M_{L} is a multiplicative term, the slope of the line in the linear relationship, Qcal is the digital number, and A_{L} is the additive linear term, the value at the origin. This will give us measurements in:

\frac{W}{m^{2} ~ srad ~ \mu m}

or Watts per square meter (surface) per steradian (solid angle) per micrometer (wavelength), which is the actual radiometric information that sensors collect.

We can accomplish this for each band by searching the metainformation text file *_MTL.txt and extracting for each band the values of RADIANCE_MULT_BAND_X and RADIANCE_ADD_BAND_X, corresponding to the parameters M_{L} and A_{L} in the above equation, respectively. I perform these operations via a Perl script, extracting the parameters using regular expressions, storing them in variables ($RADMULBAND and $RADADDBAND below), and then using the powerful NetCDF calculator grdmath as a raster calculator, in Generic Mapping Tools, and GDAL tools:

gdal_translate -of NetCDF $band.TIF $band.grd
grdmath $band.grd $RADMULBAND MUL $RADADDBAND ADD = $band-radTOA.grd
gdal_translate -ot Float32 $band-radTOA.grd $band-radTOA.tif

In this snippet, I first convert the band file to NetCDF using gdal_translate, then perform the operation with grdmath, and finally export to GeoTIFF using gdal_translate. The results, for band 8, the panchromatic, in the area around the Vesuvius volcano looks like this:

Subset of a Landsat 8, band 8 (panchromatic) image, depicting radiance at top-of-atmosphere of the Vesuvius area, on 29 July 2013.
Subset of a Landsat 8, band 8 (panchromatic) image, depicting radiance at top-of-atmosphere of the Vesuvius area, on 29 July 2013.

The image looks very similar to the uncorrected one because what we made is just a linear transformation, so the relationship between nearby pixels has not been changed. However, the image now is in 32-bit floating point format, and has double the size in memory.

For the particular case of the thermal bands (bands 10 and 11), we can proceed to calculate the brightness temperature, again at TOA, without correcting for the atmospheric distortions, by applying:

T = \frac{K_{2}}{ln(\frac{K_1}{L_{\gamma}} + 1)}

where K1 and K2 are constants called K1_CONSTANT_BAND_x and K2_CONSTANT_BAND_x, respectively, in the metainformation file. As before, we extract these using a regular expression within a script, and proceed to generate a corrected TIF using:

grdmath $K1 $band-radTOA.grd DIV 1 ADD LOG INV $K2 MUL = $band-temp.grd
gdal_translate -ot Float32 $band-temp.grd $band-temp.tif

The results are in degrees Kelvin [K], which we can easily transform to Celsius by subtracting the value at 0°C, 273.15 °K. The brightness temperature in band 10, together with the difference between band 10 and band 11 look like this:

Brightness temperatures (°C) calculated for the Vesuvius area on 29 July 2013, at satellite pass time, using band 10 (left), and difference with temperature in band 11 (right).
Brightness temperatures (°C) calculated for the Vesuvius area on 29 July 2013, at satellite pass time, using band 10 (left), and difference with temperature in band 11 (right).

As you see, the difference is very small, with mean -1.79 °C, and outliers on the top of the volcano and on some urban areas. Note how the northern slopes are slightly colder than southern ones. Pay also attention to the fact that these are uncorrected for atmospheric effects and might therefore depart from the real temperatures measured on surfaces at the same time. The above figures where drawn using R, which can also be successfully used to perform the calculations described in this page.

Finally, the Landsat use page also points to an illumination correction by which one can calculate reflectance. The procedure is the same as above, applying a linear transformation using the digital numbers and two parameters (REFLECTANCE_MULT_BAND_X, REFLECTANCE_ADD_BAND_X). Landsat 8 data comes in L1T processing level, meaning that they have already been corrected for topographic effects. According to the Landsat website this is done using a combination of SRTM, ASTER DEM, etc. These conversions work fairly well in most cases, but might not be entirely correct for polar areas, for example, where the existing DEM might be too coarse.

Advertisements