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 GeoVille - an Earth Observation Company based 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 sports.

33 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 8 years 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 8 years ago Reply


    • USgs already releases the LST tiff in their analysis ready data (ard) and can be downloaded already processed

      Aslan Varoqua 5 years 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 8 years ago Reply


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

    Thank you

    Ronaldo

    Ronaldo 7 years 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 7 years ago Reply


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

    Jimmy 7 years 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 7 years 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 7 years 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 7 years ago Reply


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

    Edmond Lai 7 years ago Reply


  • Hi
    I got the city LST from Landsat 8, however, the temperature of houses is low while it should not be so.
    would you please help me?

    freshte 6 years ago Reply


    • I don’t know where you are but the first time I produced a thermal map from Landsat 8 band 10, my assumption was that grassland would be cool and urban areas hot. So I was very surprised when it came out the other way around. However on thinking it through, the satellite goes over my site at 10 am which is about 4 hours after sunrise. At this point, only the side of any buildings facing the sun have received any direct warming. The other half is in shade. So any pixel over urban areas is going to an average of how much is in shade and how much has been directly warmed by the sun. Grasslands on the other hand, essentially provide a fairly flat but exposed boundary layer of air which is increasingly exposed to direct sunlight as the morning progresses. Other factors such as soil colour may also be affecting these grassland values. Unfortunately, I wasn’t in a position to test these explanations but they make sense of the output and were another good indication of never trust assumptions till they’re tested. I suspect if Landsat went over after midday or in the afternoon, urban areas would come out much hotter than these grasslands but unfortunately Landsat data doesn’t allow us to test this assumption.

      Dave Furniss 6 years ago Reply


      • very interesting insight!

        Martin 6 years ago Reply


  • Hi Martin,
    I would like to cite your explanations and I would like to have your last name for MLA citation.
    It would be pleasure for your cooperation.
    Thank you

    Jake 6 years ago Reply


    • Hey Jake,
      what an honour. You can use the following reference:
      Siklar, M. (2016). Calculation of Land Surface Temperature (LST) from Landsat 8 using R. [Blog] GIS-Blog.com. Available at: http://www.gis-blog.com/wp-admin/post.php?post=872 [Accessed 17 Nov. 2017].

      Let me know when you work gets published. I am sure it’s a good read.

      Cheers Martin

      Martin 6 years ago Reply


  • This is at brightness temperature…. To calculate land surface temperature, emissivity to be added .. right??

    Yunus 6 years ago Reply


  • it seems that you calculated Brightness Temperature.Land surface temperature would be calculated with emissivity and other constants:
    LST= Brightness Temp / 1 + Wevelenghth (10 Micron) * ( Brightness Temp / P) * Ln (e)
    P= h*c/s = 14380
    h= Plank constant
    c= velocity of light
    s = Boltzmann constant
    e= emissivity
    e=0.004*Pv+0.986
    Pv=proportion of vegetation= (NDVI-NDVImin / NDVImax-NDVImin)^2
    LST calculation by NDVI Threshold Method (Sobrino et al, 2014)

    Ramiin 6 years ago Reply


    • Wavelength number depends on which band do you use. For band 10 is 10.6 and for band 11 is 11.5

      Ramin 6 years ago Reply


  • Hello Martin,
    We are working on the LANDSAT 7&8, we want to download the data for the 10 years. the area under cover is 16 tiles of LANDSAT. I need a code that will help me to download the bulk images with higher downloading speed. Really appreciated your help

    Mohsin Waqas 6 years ago Reply


  • thanks a lot Martin,

    the above code is definitely brilliant if you want to know what is happening in the background. Also definitely broadens one’s coding skills.

    However, if you want to perform the above task on multiple LANDSAT files and want it done a lot quicker with minimum typos in the code, you might want to look into the option below.
    http://terradue.github.io/rLandsat8/section3/bt.html
    This approach gets rid of the background pixels. However, the result is in kelvin so you’ll want to add another line of code that converts the data to celcius and exports a .tif file

    Peter Kabano 6 years ago Reply


    • Hey Peter!
      Thanks for sharing the link. I wasn’t aware that there is a rLandsat package. I am sure that is is very useful for many readers.
      Best regards

      Martin

      Martin 6 years ago Reply


  • hello martin, thank you very much for the information, it’s great.
    Could you show how to do it with LANDSAT 5?
    Best Regards!

    Simón 5 years ago Reply


  • Hi Martin,
    I hope you are doing well.
    Could you please tell me how to calculate the Urban-rural temperature difference. Please also share any tutorial if it is available.

    Thank you

    Best Regards!
    Nayab

    Nayab 4 years ago Reply


  • Many thanks for this !!!

    Amao 4 years ago Reply


  • Thanks a lot Martin,
    But could any help me with R script to calculate Temperature Vegetation Dryness Index(TVDI) OR Vegetation Temperature Condition Index (VTCI )with LST and NDVI/EVI

    THANKS

    Abdulrashid Hassan 4 years ago Reply


  • Could you please tell me is LST depends upon month?

    Anupam Das 4 years ago Reply


  • Hi, I calculating LST for Landsat 8 but when I use export query at the end a warning message is occurring
    Warning message:
    In .gd_SetProject(object, …) : NOT UPDATED FOR PROJ >= 6
    can you help me to sort this issue?

    Mustafa526 4 years ago Reply


Post A Reply

*