The Thermal InfraRed Sensor (TIRS) is a new instrument carried onboard the Landsat 8 satellite that is dedicated to capturing temperature-specific information. Using radiation information from the two electromagnetic spectral bands covered by this sensor, it is possible to estimate the temperature at the Earth’s surface (albeit at a 100m resolution, compared to the 30m resolution of the other instrument, the Operational Land Imager).
The city-state of Delhi, India as visualised from band 10 of Landsat 8's TIRS instrument.
I used data from the TIRS to estimate the surface temperature in the city-state of Delhi, India as of the 29th of May, 2013. The relevant tarball file containing the data was downloaded using the United States’ Geological Survey’s (USGS) EarthExplorer tool; the area of interest was encompassed by [scene identifier: path 146 row 040] in the WRS-2 scheme. I think I’ll leave the specific explanations describing WRS-2, path/row values and the other miscellaneous small data-management operations for a later post. For now, I’ll just explain that these are important things to know when in the process of actually obtaining this data. When the tarball is unpacked fully, the bands from the TIRS instrument are bands 10 and 11; the relevant .tif files are [“identifier”_B10.tif] and [“identifier”_B11.tif], and these were clipped to the administrative boundary of Delhi. There’s also a text file containing metadata: [“identifier”_MTL.txt] is essential for the math we’re going to do on these two bands.
Each pixel of processed Level 1 (L1) Landsat 8 ETM+ data is stored as a Digital Number (DN) with a value between 0 and 2^15** . This means that when you visualise the two TIRS bands in your GIS/image-manipulation program of choice, you’re going to see a grey-scale image, where dark is cool and white is hot. The bands differ very slightly from each other; this is because each registers a different band of radiation, which is useful for correcting for atmospheric distortions.
To obtain actual surface temperature information, we need to first a) convert these DNs to top-of-the-atmosphere (ToA) radiance values, and then secondly b) convert the ToA radiance values to ToA brightness temperature in Kelvin. The interesting thing is that the satellite stores information as radiance values; it’s converted to DNs when processed from L0 (raw data) to L1. We’re essentially reconverting this information to its original format.
Since the math does get a bit technical, I've posted the details at the end of this blogpost.
Since the TIRS has two bands, I conducted the mathematical operation on both, and then averaged the results. I know that there are more sophisticated ways of using both bands to accurately estimate surfacte temperature; as the USGS comment to this post has noted, I’ve only derived ToA brightness temperature here. Another step is required to calculate the actual surface temperature. In fact, a fellow Landsat 8 data user, Valerio de Luca, has informed me that there’s an equation of the split-window algorithm quadratic form that can be used for this purpose; thank you! I’ll update this post again once I test that. Finally, the conversion from degrees Kelvin to degrees anything else has been sufficiently well-documented that I’m going to leave it out; I’ve converted this map into degrees Celsius.
An optional step in the middle is to correct for atmospheric distortion and solar irradiance by converting the DNs to ToA planetary reflectance values and then converting these to ToA radiance values. This requires information such as Earth-Sun distance and the solar angle, some of which is in the metadata file. I’ll be experimenting with this later as well.
Please pitch in if you have any suggestions on improving these techniques; I would definitely welcome points on better ways to use data from the TIRS instrument!
*Update 1 (245am IST on the 9th of August) : Firstly, I’d like to thank the USGS Landsat team, who’ve commented on this post; I’m very grateful for their encouragement and have updated the post accordingly. Please read the map title as ToA brightness-temperature in Delhi; another step is required before it accurately depicts the actual surface temperature values.
**Update 2 (2105 GMT on the 2nd of November) – Landsat 8 ETM+ Digital Number range has been corrected to reflect the fact that the ETM+ sensor stores 16-bit data; thanks for the comment Ciprian!
***Update 3 (0110 IST on the 5th of November 2017) - This article has been cross-posted from geohackers.in
I’ve referred to the documents below while figuring the math and writing this blogpost; they’re very useful.
USGS equations re: Landsat 8 data
https://landsat.usgs.gov/Landsat8_Using_Product.php
Converting Digital Numbers to Kelvin
http://www.yale.edu/ceo/Documentation/Landsat_DN_to_Kelvin.pdf
Information re: ToA reflectance
http://www.pancroma.com/Landsat%20Top%20of%20Atmosphere%20Reflectance.html
GRASS methodology
http://grass.osgeo.org/grass64/manuals/i.landsat.toar.html
More information on the inverted Planck function
http://www.yale.edu/ceo/Documentation/ComputingThePlanckFunction.pdf
Markham et al. (2005): Processing data from the EO-1 ALI instrument
http://www.asprs.org/a/publications/proceedings/pecora16/Markham_B.pdf
Converting from spectral radiance to planetary reflectance
http://eo1.usgs.gov/faq/question?id=21
<math>
a) Digital Numbers to ToA radiance values
Lλ = (ML * Qcal) + AL – - equation I
where:
Lλ = TOA spectral radiance
ML = Band-specific multiplicative rescaling factor
AL = Band-specific additive rescaling factor
Qcal = Quantized and calibrated standard product pixel values (DN)
Human-readable format:
Spectral Radiance = Band-specific multiplicative rescaling factor into Digital Number plus Band-specific additive rescaling factor from the metadata1. From the Markham et. al. (2005) paper on processing data from the EO-1 ALI instrument, I’m led to understand that the multiplicative and additive scaling factors result in ‘better preservaton of the ALI’s precision and dynamic range in going from raw to calibrated data’; this probably holds true for the L8 instruments as well but do correct me if I’m wrong.
In practise, this is:
[ Spec_rad_B“X”.tif ] = RADIANCE_MULT_BAND_”X” * [“identifier”_B”X”.tif] + RADIANCE_ADD_BAND_”X”
where “X” is the specific band number, the RADIANCE values are obtained from the metadata file and Spec_rad_B”X”.tif is the output spectral radiance file name (user-determined).
b) ToA radiance values to ToA brightness-temperature in Kelvin
T = K2 / (ln (K1/ Lλ +1))[2] – equation II
(The physics-inclined will realise that equation II above is the inverted form of the Planck function.)
where:
T = ToA brightness temperature (K)
Lλ = TOA spectral radiance
K1 = Band-specific thermal conversion constant
K2 = Band-specific thermal conversion constant
Human-readable format:
Temperature (Kelvin) = Band-specific thermal conversion constant I by the natural log of (Band-specific thermal conversion constant II by spectral radiance plus 1)
In practise, this is:
[Temp_B”X”.tif] = K2_CONSTANT_BAND_”X” / ln((K1_CONSTANT_BAND_”X” / [Spec_rad_B”X”.tif]) + 1)
where “X” is the specific band number, the CONSTANT values are extracted from the metadata file and the Temp_B”X”.tif file is the output temperature file name (user-determined).
</math>