2014/01/19

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

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 Surface Temperature
  2. Land cover classification of study area
  3. Reclassification of the land cover classification to emissivity values
  4. 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

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 (ε) 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 ε
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: minimummaximum, and new. Minimum and maximum define the reclassification range for the new value, according to the following table. 
 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).
  • 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 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 (ε) values of various land cover types are provided in the following table (from Mallick, et al. 2012).
Land surface
Emissivity ε
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 / ρ) lnε ]
where:
  • λ = 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)
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

23 comments:

  1. wow!!!, i will try this, thanks Luca

    ReplyDelete
  2. Hi Sr Congedo,

    I have a small issue with your plugin and Dufour version of Qgis running on Windows 7 32bits.
    Conversion stopped after first view.
    Do you know why? Thanks.

    My best,
    Julien

    ReplyDelete
    Replies
    1. Hi Julien,
      please check if the plugin is correctly installed and configured as explained here http://fromgistors.blogspot.com/2013/09/install-semi-automatic-classification-windows.html

      Also, could you provide more details about the issue and send me the log file (as described here http://fromgistors.blogspot.com/p/ask-for-help.html).

      Thank you.

      Delete
  3. It would be great to add this kind of thing to the QGIS training manual!

    docs.qgis.org/2.0/en/docs/training_manual/

    ReplyDelete
    Replies
    1. Thank you for the suggestion.
      I'll talk to the QGIS developers and see if it can be added.

      Cheers!

      Delete
  4. Hello Luca!,

    I saw you have a few ROI per class and I am wondering if it is possible to combine them? With an aim to have one curve per class (which include all the ROI of that class) on the signature plot. This could be great.

    By the way, is it normal to have a python error when I click on the "Edit curves line and axes parameters" button on the ROI tools (with virtual machine semi-auto QGIS).

    Cheers!

    Mathieu

    ReplyDelete
  5. Hello Mathieu,

    at the moment the spectral signatures are calculated for each ROI. However, I am planning to add something like your request in the next version of the plugin.

    About the error, yes, it is a known issue.

    Cheers!

    ReplyDelete
    Replies
    1. Hi Luca!

      I think it will be really great to introduce a separability index for the classes. Maybe something like a divergence index or the bhattacharyya distance (like in PCI-Focus) : it is so helpful when you have a lot a classes. This is just an idea...thanks!

      Delete
    2. Hi Mathieu,

      I think it is a good idea.
      I'll try to add such a function to the plugin.
      Thank you for your feedback.

      Delete
  6. I saw that you corrected the LST using the emissivity. But what about atmospheric corrections ? See the Qin (2001) and Sobrino (2004) algorithms.

    ReplyDelete
    Replies
    1. Hello Niros,
      thank you for the references. As you can see in the yellow frame at the end of the post Sobrino (2004) was already mentioned.
      The atmospheric corrections described in the articles (mono-window, single-channel) require some additional information such as the total atmospheric water vapor, which are not always available.
      However, the purpose of this post is to show a rapid method for estimating surface temperature, based only on the image. Also, the DOS1 atmospheric correction is applied for the land cover classification.
      Fortunately, the Landsat 8 images will allow for the use of split window algorithms although it is not suggested yet (http://landsat.usgs.gov/calibration_notices.php).

      Delete
    2. Yes Luca.

      Estimating LST with no atmospheric corrections is an incomplete job, according to scientific literature ( see Voogh & Oke paper here: http://www.gisknowledge.net/topic/eo_by_remote_sensing/voogt_oke_rse_03.pdf ).

      This is why I was searching to find some steps or explanations about how to apply the atmospheric corrections, using the Qin (2001) or the Jimenez-Munoz & Sobrino (2004) algorithms. And this is how I found your post here, wich is very helpful !

      But how to obtain the total atmospheric water vapor, or how to calculate it ?
      Can you give me more details about this ?

      Thank you.

      Delete
    3. I would not say that no atmospheric correction is an incomplete job. I would rather say that it can lead to lower accuracy of measures.
      The need to perform an atmospheric correction depends on the level of accuracy that is required for your study.
      I am glad you found this post helpful.
      About the atmospheric water vapor and other parameters required for atmospheric correction you can find a very useful calculator here http://atmcorr.gsfc.nasa.gov/

      Delete
  7. Hello, i have a problem when i try to pick up Rois. I get a python error. Kindly check the message and advice.

    An error has occured while executing Python code:

    Traceback (most recent call last):
    File "C:/Users/user1/.qgis2/python/plugins\SemiAutomaticClassificationPlugin\semiautomaticclassificationplugin.py", line 1022, in pointerClickROI
    self.createROI(self.pntROI)
    File "C:/Users/user1/.qgis2/python/plugins\SemiAutomaticClassificationPlugin\semiautomaticclassificationplugin.py", line 1235, in createROI
    self.regionGrowing(dictBands, point.x(), point.y(), self.rngRad, self.minROISz, tmpShp)
    File "C:/Users/user1/.qgis2/python/plugins\SemiAutomaticClassificationPlugin\semiautomaticclassificationplugin.py", line 1295, in regionGrowing
    refRstrCols = refRstrDt.RasterXSize
    AttributeError: 'NoneType' object has no attribute 'RasterXSize'

    Python version:
    2.7.4 (default, Apr 6 2013, 19:54:46) [MSC v.1500 32 bit (Intel)]


    QGIS version:
    2.2.0-Valmiera Valmiera, c3a2817

    Python path: ['C:/PROGRA~1/QGISVA~1/apps/qgis/./python/plugins\\processing', 'C:/PROGRA~1/QGISVA~1/apps/qgis/./python', u'C:/Users/user1/.qgis2/python', u'C:/Users/user1/.qgis2/python/plugins', 'C:/PROGRA~1/QGISVA~1/apps/qgis/./python/plugins', 'C:\\PROGRA~1\\INTERG~1\\ERDASI~1\\usr\\lib\\Win32Release\\python', 'C:\\PROGRA~1\\QGISVA~1\\bin\\python27.zip', 'C:\\PROGRA~1\\QGISVA~1\\apps\\Python27\\DLLs', 'C:\\PROGRA~1\\QGISVA~1\\apps\\Python27\\lib', 'C:\\PROGRA~1\\QGISVA~1\\apps\\Python27\\lib\\plat-win', 'C:\\PROGRA~1\\QGISVA~1\\apps\\Python27\\lib\\lib-tk', 'C:\\PROGRA~1\\QGISVA~1\\bin', 'C:\\PROGRA~1\\QGISVA~1\\apps\\Python27', 'C:\\PROGRA~1\\QGISVA~1\\apps\\Python27\\lib\\site-packages', 'C:\\PROGRA~1\\QGISVA~1\\apps\\Python27\\lib\\site-packages\\PIL', 'C:\\PROGRA~1\\QGISVA~1\\apps\\Python27\\lib\\site-packages\\win32', 'C:\\PROGRA~1\\QGISVA~1\\apps\\Python27\\lib\\site-packages\\win32\\lib', 'C:\\PROGRA~1\\QGISVA~1\\apps\\Python27\\lib\\site-packages\\Pythonwin', 'C:\\PROGRA~1\\QGISVA~1\\apps\\Python27\\lib\\site-packages\\Shapely-1.2.18-py2.7-win32.egg', 'C:\\PROGRA~1\\QGISVA~1\\apps\\Python27\\lib\\site-packages\\wx-2.8-msw-unicode', 'C:\\PROGRA~1\\QGISVA~1\\apps\\qgis\\python\\plugins\\fTools\\tools']

    ReplyDelete
    Replies
    1. Hello James,

      please report the log file as explained here http://fromgistors.blogspot.com/p/ask-for-help.html.

      Delete
  8. Can I automate the process using python script e.g. select the particular band/s file from a directory and apply atmospheric correction.
    Is there any python package available for doing this?
    Thank You.

    ReplyDelete
    Replies
    1. Hello,
      I do not know of any python package for atmospheric correction.
      However, you can use the Semi-Automatic Classification Plugin (which is written in Python) for a DOS1 correction of bands. Also, you could take portion of the code and adapt the script to your needs with Python.

      Delete
    2. Thanks for your reply. "Semi-Automatic" gr8 work it will really benefit the QGIS user.
      Regards,
      Suryakant

      Delete
  9. hello luca, it is ilyas from the same field

    ReplyDelete
  10. i want to know about the general variation which is possible between lst derived through r/s technique and field based estimation

    ReplyDelete
    Replies
    1. Hello Ilyas,

      the estimation of LST depends on the satellite, the method of estimation and the atmospheric correction applied.
      Generally speaking the variation between estimated LST (e.g. with the method of this tutorial) and field temperature can be of + - 3 or 4K. With particularly advanced techniques the variation can be + - 1K, but they require atmospheric data for accurate atmospheric correction.

      Delete