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.
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).
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:
- Conversion of raster bands from DN to reflectance and At Surface Temperature
- Land cover classification of study area
- Reclassification of the land cover classification to emissivity values
- Conversion from At Surface 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
- 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.
|Landsat bands converted to TOA reflectance and brightness temperature|
2. Land cover classification of the study area
The following are the main land cover classes that we are going to identify in the image:
- 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.
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 valuesNow we are going to reclassify the classification raster using the land surface emissivity values.
The emissivity (ε) 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):
- From the Processing toolbox navigate to SAGA > Grid - Tools > Reclassify grid values; in the tool window, under Grid select classification.tif; under method select  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.
4. Conversion from At Satellite Temperature to Land Surface Temperature
- 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).
- Under Formula write:
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.
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):
- K1 = Band-specific thermal conversion constant (in watts/meter squared * ster * μm)
- K2 = Band-specific thermal conversion constant (in kelvin)
- 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)
K1 (watts/meter squared * ster * μm)
** 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 (ε) values of various land cover types are provided in the following table (from Mallick, et al. 2012).
- Land surfaceEmissivity εSoil0.928Grass0.982Asphalt0.942Concrete0.937
- λ = wavelength of emitted radiance
- ρ = h * c / σ (1.438 * 10^-2 m K)
- h = Planck’s constant (6.626 * 10^-34 Js)
- σ = Boltzmann constant (1.38 * 10^-23 J/K)
- c = velocity of light (2.998 * 10^8 m/s)
- SatelliteBandCenter wavelength (μm)Landsat 4, 5, and 7611.45Landsat 81010.8Landsat 81112
-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