Estimation of Land Surface Temperature with Landsat Thermal Infrared Band: a Tutorial Using the Semi-Automatic Classification Plugin for QGIS


 Updated tutorial at https://fromgistors.blogspot.com/2016/09/estimation-of-land-surface-temperature.html


This post is a tutorial for the estimation of Land Surface Temperature using a Landsat image acquired over Paris (France), using the Semi-Automatic Classification Plugin for QGIS, which allows for supervised classifications.

Before the tutorial, please watch the following video that illustrates the study area and provides very useful information about thermal infrared images, and their application (footage courtesy of European Space Agency/ESA). Also, a brief description of the area that we are going to classify is available here.
As shown in the previous video, the study area is covered by the urban surfaces, vegetation and agricultural fields. The thermal infrared band is particularly useful for assessing the temperature difference between the city and the surrounding rural areas, and studying the urban heat island phenomenon.
We are going to classify a Landsat 8 image acquired in September 2013 (available from the U.S. Geological Survey), using a semi-automatic approach (i.e. the entire image is classified automatically basing on several samples of pixels "ROIs" that are used as reference classes), and therefore estimate Land Surface Temperature using the thermal band and the land cover classification (used for the definition of surface emissivity).

The Semi-Automatic Classification Plugin allows for the automatic conversion of Landsat thermal band to At-Satellite Brightness Temperature (as illustrated in the yellow frame at the end of this post), measured at the sensor.
There are several studies about the calculation of the Land Surface Temperature. For instance, Landsat 8 provides two thermal bands (bands 10 and 11) that could be used with split-window methods. However, "given the larger uncertainty in the Band 11 values, users should work with TIRS Band 10 data as a single spectral band (like Landsat 7 Enhanced Thematic Mapper Plus (ETM+)) and should not attempt a split-window correction using both TIRS Bands 10 and 11" (from http://landsat.usgs.gov/calibration_notices.php).

Alternatively, there are methods that use the NDVI for the estimation of land surface emissivity, or using a land cover classification for the definition of the land surface emissivity of each class (see the yellow frame at the end of this post).

In this tutorial we are going to use a land cover classification for the definition of surface emissivity, which is required for the calculation of the Land Surface Temperature. The following phases are required:
  1. Conversion of raster bands from DN to reflectance and At Satellite Temperature
  2. Land cover classification of study area
  3. Reclassification of the land cover classification to emissivity values
  4. Conversion from At Satellite Temperature to Land Surface Temperature


 First, download the image from here (available from the U.S. Geological Survey). The image is a subset of the entire scene, including the following Landsat bands (each band is a single 16 bit raster):
2 - Blue;
3 - Green;
4 - Red;
5 - Near-Infrared;
6 - Short Wavelength Infrared 1;
7 - Short Wavelength Infrared 2;
10 - Thermal Infrared (TIRS) 1.

1. Conversion of raster bands from DN to reflectance and At Satellite Temperature

In this step, we are going to convert the DN (i.e. Digital Numbers) to the physical measure of Top Of Atmosphere reflectance (TOA) and the thermal band to At-Satellite Brightness Temperature. Also, we are going to apply a simple atmospheric correction to the bands from 2 to 7 (i.e. the DOS1 method - Dark Object Subtraction 1 - which is an image-based technique).
  • Open QGIS and start the Semi-Automatic Classification Plugin; in the main interface select the tab Pre processing > Landsat;
  • Select the directory that contains the Landsat bands (and also the required metafile MTL.txt), and select the output directory where converted bands are saved;
  • Check the option Apply DOS1 atmospheric correction, and click Perform conversion to convert Landsat bands to reflectance.

At the end of the process, all the converted bands are loaded in QGIS.
Landsat bands converted to TOA reflectance and brightness temperature

2. Land cover classification of the study area

This is a required step for the conversion from At-Satellite Brightness Temperature to Land Surface Temperature, because we are going to use a land cover classification for the definition of the land surface emissivity.
The following are the main land cover classes that we are going to identify in the image:
  1. Built-up
  2. Soil
  3. Vegetation
  4. Water
In order to perform the land cover classification we need to:
  • Define the classification inputs, which are Landsat bands from 2 to 7 (excluding the thermal band 10);
  • Create the ROIs;
  • Classify the study area, using the Spectral Angle Mapping Algorithm.
For detailed instructions about the supervised classification using the Semi-Automatic Classification Plugin, see my previous post here.

 At the end of the process, we have a land cover classification of the area, as shown in the following image (you can download the ROIs and the land cover classification from here).
Supervised classification of Landsat image

3. Reclassification of the land cover classification to emissivity values

Now we are going to reclassify the classification raster using the land surface emissivity values.
The emissivity (e) values for the land cover classes are provided in the following table (these values are only indicative, because they should be obtained from field survey):
Land surface
Emissivity e
Soil
0.93
Vegetation
0.98
Built-up
0.94
Water
0.98

 In QGIS:
  • From the Processing toolbox navigate to SAGA > Grid - Tools > Reclassify grid values; in the tool window, under Grid select classification.tif; under method select [2] simple table;
  • Under Lookup table click the button ... to open the window Fixed Table; the lookup table has 3 columns: minimum, maximum, and new. Minimum and maximum define the reclassification range for the new value, according to the following table. 
#UPDATE 09/05/2015: alternatively it is possible to use the tool Reclassification of the Semi-Automatic Classification Plugin  v.4 . See http://semiautomaticclassificationmanual-v4.readthedocs.org/en/latest/main_interface_window.html#reclassification
 minimum  maximum  new 
000Unclassified
11190.94Built-up
21290.93Soil
31390.98Vegetation
41490.98Water
  • Select where to save the emissivity raster, and click "OK".

After a few seconds the reclassified grid (emissivity) will be loaded in QGIS.
Emissivity raster

4. Conversion from At Satellite Temperature to Land Surface Temperature

Now we are ready to convert the At-Satellite Brightness Temperature to Land Surface Temperature.
In QGIS:
  • From the Processing toolbox navigate to SAGA > Grid - Calculus > Raster calculator; under Raster layers select the emissivity raster and the brightness temperature raster (the emissivity raster must be above the brightness temperature raster).
#UPDATE 09/05/2015: alternatively it is possible to use the tool Band calc of the Semi-Automatic Classification Plugin v.4 . See http://semiautomaticclassificationmanual-v4.readthedocs.org/en/latest/main_interface_window.html#band-calc
  • Under Formula write:
b / ( 1 + ( 10.8 * b / 14380 ) * ln(a) )

 where a is the emissivity raster and b is the brightness temperature raster (for the details about this formula, see the yellow frame at the end of this post).

 After a few seconds, the Land Surface Temperature (in kelvin) will be displayed, as shown in the following figure. You can download the temperature raster from here.
Land Surface Temperature (in kelvin)

 Please notice that the aim of this tutorial is to show how surface temperature can be estimated using open source programs and free images, and the results are for demonstration purposes. In order to achieve better results, one should perform field survey for improving the land cover classification of the area and the measurement of surface emissivities.

 The following is the video of this tutorial.




Conversion to At-Satellite Brightness Temperature
For Landsat thermal bands, the conversion of DN to At-Satellite Brightness Temperature is given by (from https://landsat.usgs.gov/Landsat8_Using_Product.php):
TB = K2 / ln( (K1/Lλ)+ 1)
where:
  • K1 = Band-specific thermal conversion constant (in watts/meter squared * ster * µm)
  • K2 = Band-specific thermal conversion constant (in kelvin)
and Lλ is the Spectral Radiance at the sensor's aperture, measured in watts/(meter squared * ster * µm); for Landsat images it is given by (from https://landsat.usgs.gov/Landsat8_Using_Product.php)
Lλ = ML * Qcal + AL
where:
  • ML = Band-specific multiplicative rescaling factor from Landsat metadata (RADIANCE_MULT_BAND_x, where x is the band number)
  • AL = Band-specific additive rescaling factor from Landsat metadata (RADIANCE_ADD_BAND_x, where x is the band number)
  • Qcal = Quantized and calibrated standard product pixel values (DN)
The K1 and K2 constant for Landsat sensors are provided in the following table:
Constant
Landsat 4*
Landsat 5*
Landsat 7**
K1 (watts/meter squared * ster * µm)
671.62
607.76
666.09
K2 (Kelvin)
1284.30
1260.56
1282.71
* from Chander & Markham (2003)
** from NASA (2011)

 For Landsat 8, the K1 and K2 values are provided in the image metafile.

 Estimation of Land Surface Temperature
There are several studies about the calculation of land surface temperature. For instance, using NDVI for the estimation of land surface emissivity (Sobrino, et al., 2004), or using a land cover classification for the definition of the land surface emissivity of each class (Weng, et al. 2004).
For instance, the emissivity (e) values of various land cover types are provided in the following table (from Mallick, et al. 2012).
Land surface
Emissivity e
Soil
0.928
Grass
0.982
Asphalt
0.942
Concrete
0.937
Therefore, the land surface temperature can be calculated as (Weng, et al. 2004):

T = TB / [ 1 + (λ * TB / p) ln(e) ]
where:
  • λ = wavelength of emitted radiance
  • p
     = h * c / s     (1.438 * 10^-2 m K)
  • h = Planck’s constant (6.626 * 10^-34 Js)
  • s = Boltzmann constant (1.38 * 10^-23 J/K)
  • c = velocity of light (2.998 * 10^8 m/s)
The values of λ  for the thermal bands of Landsat setellites are listed in the following table:
Satellite
Band
Center wavelength (µm)
Landsat 4, 5, and 7
6
11.45
Landsat 8
10
10.8
Landsat 8
11
12
References
-Chander, G. & Markham, B. 2003. Revised Landsat-5 TM radiometric calibration procedures and postcalibration dynamic ranges Geoscience and Remote Sensing, IEEE Transactions on, 41, 2674 – 2677
-Mallick, J.; Singh, C. K.; Shashtri, S.; Rahman, A. & Mukherjee, S. 2012. Land surface emissivity retrieval based on moisture index from LANDSAT TM satellite data over heterogeneous surfaces of Delhi city International Journal of Applied Earth Observation and Geoinformation,, 19, 348 - 358
-NASA (Ed.) 2011. Landsat 7 Science Data Users Handbook Landsat Project Science Office at NASA's Goddard Space Flight Center in Greenbelt, 186
-Sobrino, J.; Jiménez-Muñoz, J. C. & Paolini, L. 2004. Land surface temperature retrieval from LANDSAT TM 5 Remote Sensing of Environment, Elsevier, 90, 434-440
-Weng, Q.; Lu, D. & Schubring, J. 2004. Estimation of land surface temperature–vegetation abundance relationship for urban heat island studies. Remote Sensing of Environment, Elsevier Science Inc., Box 882 New York NY 10159 USA, 89, 467-483
Newer posts Older posts