Calculation of Land Surface Temperature (LST) from Landsat 8 using R

Lake Baikal

Hey there!
I was recently asked by a user how to calculate the Land Surface Temperature (LST) from Landsat 8 imagery and decided to write an article on this topic. I hope you will enjoy it.

Let’s start with the basics – What is LST?


The Land Surface Temperature (LST) is the radiative skin temperature of ground. It depends on the albedo, the vegetation cover and the soil moisture. In most cases, LST is a mixture of vegetation and bare soil temperatures. Because both respond rapidly to changes in incoming solar radiation due to cloud cover and aerosol load modifications and diurnal variation of illumination, the LST displays quick variations too. In turn, the LST influences the partition of energy between ground and vegetation, and determines the surface air temperature.

Definition by Copernicus Global Land Service

There are various methods to calculate LST. I am going to show you the method from the official USGS Webpage using  Bands 10 and 11 from the Thermal Infrared Sensor (TIRS) of the Landsat 8 satellite. It is recommended to use Band 10 in quantitative analysis because Band 11 is significantly more contaminated by stray light than Band 10.

Step 1: Raw image

Let’s start with a look on our raw image downloaded from one of the Landsat mirrors. My image is from July 25th 2015 and shows the majestic Baikal Lake in Siberia, Russia with its surrounding vegetation.(RGB 432):

Lake Baikal

Nature is very generous to Lake Baikal. A necklace of coniferous forests, decorating the shore mountain chains, creates a unique beauty of landscapes, shelters animals and birds, and enlivens the surroundings. (BWW.irk.ru)

Lake Baikal is the largest freshwater lake by volume in the world, containing roughly 20% of the world’s unfrozen surface fresh water. It is characterised by a very heterogenous environment and is surrounded by vast forest areas, bush and grasslands and a very dry steppe areas. People might think that it’s cold in Siberia, but the annual mean temperature at Lake Baikal is actually 11.8 °C. July is the hottest months with an average temperature of 23°C. (Let’s keep those values for our LST analysis in mind.)

Step 2: Convert Digital Numbers (DN) to Top Of Atmosphere (TOA) Reflectance

For our analysis and calculation of the LST we need Band 10 and Band 11 only. The values in the the downloaded raw Landsat 8 images are so called Digital Numbers of the sensor. First, these have to be converted to TOA reflectance values using the following equation:

DN to TOA Reflectance

from http://landsat.usgs.gov/Landsat8_Using_Product.php

To solve this equation you need the RADIANCE_MULT_BAND_x and RADIANCE_ADD_BAND_x from the Metafile that is included in every Landsat 8 download. These values (can) differ from image to image, so make sure you look them up for every scene!! Just open the metafile in a text editor and look for the corresponding values for Band 10 and Band 11:
Metafile

Now we can convert the values in R as follows:

#Values from Metafile
RADIANCE_MULT_BAND_10 <- 3.3420E-04
RADIANCE_MULT_BAND_11 <- 3.3420E-04

RADIANCE_ADD_BAND_10 <- 0.10000
RADIANCE_ADD_BAND_11 <- 0.10000

#Load raster package and load band 10 & 11 into R (navigate to your image directory first)
library(raster)
band_10 <- raster("band10.tif") #change image name accordingly
band_11 <- raster("band11.tif") #change image name accordingly

#Calculate TOA from DN:
toa_band10 <- calc(band_10, fun=function(x){RADIANCE_MULT_BAND_10 * x + RADIANCE_ADD_BAND_10})
toa_band11 <- calc(band_11, fun=function(x){RADIANCE_MULT_BAND_11 * x + RADIANCE_ADD_BAND_11})

Step 3:  Conversion to At-Satellite Brightness Temperature

In the next step we convert our reflectance values into Satellite Brightess Temperature, which in this case is the LST in Kelvin. This is done with the following equation and also needs some data from the metafile:

Reflectance to LST

http://landsat.usgs.gov/Landsat8_Using_Product.php

Once you found the corresponding values for K1_CONSTANT_BAND_10, K2_CONSTANT_BAND_10, K1_CONSTANT_BAND_11 and K2_CONSTANT_BAND_11 in your metafile, you can calculate the LST with R as follows:

#Values from Metafile
K1_CONSTANT_BAND_10 <- 774.8853
K1_CONSTANT_BAND_11 <- 480.8883
K2_CONSTANT_BAND_10 <- 1321.0789
K2_CONSTANT_BAND_11 <- 1201.1442

#Calculate LST in Kelvin for Band 10 and Band 11
temp10_kelvin <- calc(toa_band10, fun=function(x){K2_CONSTANT_BAND_10/log(K1_CONSTANT_BAND_10/x + 1)})
temp11_kelvin <- calc(toa_band11, fun=function(x){K2_CONSTANT_BAND_11/log(K1_CONSTANT_BAND_11/x + 1)})

#Convert Kelvin to Celsius for Band 10 and 11
temp10_celsius <- calc(temp10_kelvin, fun=function(x){x - 273.15})
temp11_celsius <- calc(temp11_kelvin, fun=function(x){x - 273.15})

#Export raster images
writeRaster(temp10_celsius, "temp10_c.tif")
writeRaster(temp11_celsius, "temp11_c.tif")

Step 4: Show final result in QGIS

And this is how the result looks like. The following image depicts the LST values in Celsius for the Band 10 in the area around the Baikal Lake in Siberia from July 25th, 2015:

LST in Celsius

You can see a cloud mask (white) and the cold blue water body of the deep Baikal Lake with an average water surface temperature of about 5-6 °C. The surrounding bushlands and forest have a LST of about 18-22 °C.  The dry steppe and sand areas in the southern part of the image have a very high temperature (since the image was taken around noon on a hot summer day.)

I hope this tutorial was useful to some of you. If you have any questions regarding the code or the conversion, please dont hesitate to ask in the comment section. For more details you can visit the following websites:

 

Cheers Martin

About This Author

Martin was born in Czech Republic and studied at the University of Natural Resources and Life Sciences, Vienna. He is currently working at an Earth Observation Company in Austria, specialised in Land Monitoring. His main interests are: Open-source applications like R, (geospatial) statistics and data-management, web-mapping and visualization. He loves travelling, geocaching, photography and TV series.

13 Comments

You can post comments in this post.


  • Hi,
    I want to calculate LST for Landsat8 data using the above codes in R but I have faced with the following error.
    I would be grateful if you help me.
    #Calculate TOA from DN:
    toa_band10 <- calc(band_10, fun=function(x){RADIANCE_MULT_BAND_10 * x + RADIANCE_ADD_BAND_10})
    ERROR:
    unable to find an inherited method for function 'calc' for signature '"character", "function"'

    maliheh 1 year ago Reply


    • I would need to see more of your code. It seems that either your raster wasn’t read into R properly, or your variabls are not integers but characters. Can you post the output of the following variables?
      toa_band10
      RADIANCE_MULT_BAND_10
      RADIANCE_ADD_BAND_10

      Martin 1 year ago Reply


      • Thanks. It was right. My raster data wasn’t read into R properly. I could create it at the end.

        maliheh 1 year ago Reply


  • It is possible to use this code in order to convert Thermal band of Landsat TM5 to LST product or it is only exclusive to OLI data.
    thanks

    maliheh 1 year ago Reply


  • The calculations done correctly. How do I get the image in QGIS ?
    What color combination used ?

    Thank you

    Ronaldo

    Ronaldo 11 months ago Reply


  • Tif fo Lake Titicaca in Peru.
    The following link generated tif ( temp10_c.tif ) .

    https://drive.google.com/file/d/0B8q_146s7i16VUo1clpfQTY1a2M/view?usp=sharing

    False-Colour:

    https://drive.google.com/file/d/0B8q_146s7i16bnJFNV9iaVVheUE/view?usp=sharing

    Can you help me get a tif similar to this blog

    Thank you

    Ronaldo 11 months ago Reply


  • Hi,Martin. Could you please tell me how to generate cloud mask? Thank you very much.

    Jimmy 8 months ago Reply


  • In step 3 you state, “Satellite Brightess Temperature, which in this case is the LST in Kelvin”. Aren’t these different things? It looks like you are calculating at-sensor temperature, but I thought this is different than land surface temperature (LST) because of emissivity and atmospheric effects.

    Mankoff 4 months ago Reply


    • Exactly my point, the value being calculated here is at sensor brightness temperature. For calculating Land Surface Temperature other methods such as split window algorithm, Temperature emissivity separation etc. are to be used.

      Reema 2 months ago Reply


  • Dear Sir,
    thanks for helping us.
    could you please help me to find out LAI leaf area index using R.?
    and if you could help me to get the Albedo calculation using R?

    Regards

    Ahmad shah 3 weeks ago Reply


  • Thank you for showing us your code, it is very helpful!

    Edmond Lai 3 weeks ago Reply


Post A Reply

*