Deriving temperature from Landsat 8 thermal bands (TIRS)

You can use the thermal bands from Landsat 8 to calculate at-satellite brightness temperature. There are quite a few steps that I’ll walk you th­­rough to do this. What we’re going to do is:

-        Convert the raw bands into Top of Atmosphere Radiance (TOAr).

-        Convert TOAr into degrees kelvin.

-        Convert degrees kelvin into degrees Fahrenheit.

-        Save and export this workflow as a template that you can apply it to any other image.

Step 1 – The mise en place of GIS

A – Luck is where preparation meets opportunity

This is a fairly lengthy blog, so I’d recommend topping off that cup of coffee or perhaps stretching the legs before getting started.

B – Raster Function Template Editor

Open the Raster Function Template Editor from the toolbar. Right click on the raster and insert a Composite Band Function. Click on the plus sign and select Add copy of selected input. You have to do this because there are two thermal bands for Landsat 8. If you’re working with Landsat 5 or 7, you can omit this step.

[If you don’t have the Raster Function Template Editor loaded, you can add it through the customize window, and it will be under the Commands tab, in the raster section. Drag and drop it into the toolbar.]

raster function editor interface

C – Preview of the final output

Here’s what the final output will look like once you have all of the steps completed that are shown below. Note that the first branch is for Band 10 and the second is Band 11. This is to match the structure of the Thermal Raster Product for Landsat 8. The type of operation is shown inside of the brackets. While you can label the steps however you wish you do need to follow this structure, so that the chain works when you apply it to other Landsat 8 scenes.

preview of the final output

D – Reading the metadata

All of the inputs that you need to do this will be in the metadata file. In fact, the inputs that you will need to perform this conversion are the same for all of the Landsat 8 images. Here’s a table that you can use rather than combing through the metadata files.

Band 10

Band 11

Radiance Multiplier

0.0003342

0.0003342

Radiance Add

0.1

0.1

K1

774.89

480.89

K2

1321.08

1201.14

 

Step 2 – Conversion from DN to Radiance

The Radiance Multiplier and Radiance Add are used to convert the DN back into radiance. Remember from back in the day, that the equation of a line is y = mx + b? Well, that’s all you’re doing here.

-        y is going to be the TOAr that you want.

-        m is the Radiance Multiplier

-        x is the raw band

-        b is the Radiance Add.

Right click on Raster and insert Band Arithmetic. Set the method to user defined, and use b1 to refer to Input Raster. Repeat this for the other band.

conversion from DN to radiance with band arithmetic

Step 3 – Convert Radiance into degrees kelvin

From TOAr, we now need to convert to temperature in degrees kelvin. This is where you use the K1 and K2 inputs. The equation we will be using functions to create is:

K2 / ln(k1/TOAr + 1)

The easiest way to do this is to create the (k1/TOAr + 1) part first using the Band Arithmetic function. (Right click on the Band Arithmetic that you created in the previous step and insert another Band Arithmetic function. This will appear above the previous step.) Repeat this for each band, but keep in mind that K1 is different for bands 10 and 11.

converting radiance to kelvin - natural log input

To calculate the natural log, you need to use local function. (Right click on the Band Arithmetic from the previous step and insert a local function. The dialog below will appear. Select Ln under the Math section, and click OK.)

conversion from radiance to kelvin - applying ln with local function

In a separate step, divide that output into the K2 constant, using Band Arithmetic. (Right click on the local function you created in the previous step, and insert a Band Arithmetic Function.) Repeat this for each band, but keep in mind that K2 is different for bands 10 and 11.

conversion to kelvin from radiance - final step

Step 4 – Convert degrees kelvin into Fahrenheit

Now you have converted the TOAr into degrees kelvin. If you want to convert this into Celsius, you need to subtract 273.15 from the degrees kelvin. For Fahrenheit, multiply the degrees Celsius by 1.8 and then add 32. You can do this as one step within the Band Arithmetic function.

Conversion to degrees Fahrenheit with Band Arithmetic

Step 5 – Applying this function chain to any Landsat 8 image

When you’re finished, save this as a raster function template.xml file. Then, open the Edit Raster Function geoprocessing tool. Select the Landsat 8 thermal raster product that you want to convert to degrees Fahrenheit as the input raster. Choose the Raster Function Template that you just saved, and leave the Function Name blank. Repeat this for as many Landsat 8 scenes as you want.

Edit Raster Function Geoprocessing Tool interface

Here’s an example of two images taken on December 16, 2013. The first is from sunny southern California where we’ve been experiencing unseasonably warm temperatures. The other is from “rustic” New Brunswick, where some Esri employees decide to vacation, and which is experiencing the coldest winter in 20 years. The deep blues in each image are clouds, and have negative values. Note that we have not derived the land surface temperature, only the at-satellite brightness temperature. land surface temperature redlands, californialand surface temperature new brunswick, canada

This entry was posted in Analysis & Geoprocessing, Imagery and tagged , , , , . Bookmark the permalink.

Leave a Reply

12 Comments

  1. miki628 says:

    I failed.
    Apply function when execute Edit Raster Function tool without add the Landsat 8 thermal raster product to Content window. ArcMap crash when execute Edit Raster Function tool from add the Landsat 8 thermal raster product to Content window.
    Although I went to the same way, why can not ?

    • alexkostin says:

      Instal first ArcGIS 10.2 for (Desktop, Engine, Server) Landsat 8 Patch from http://support.esri.com/en/downloads/patches-servicepacks/view/productid/66/metaid/2012
      I got, but the result is showing RGB data, not temperature. If switch from classified to stretched mode – it looks like temperature.
      I tried 2 ways: 1 – that is written above, 2 – using formulas for BAND 10: 1321.08/ln(774.89/(0.0003342*b1+0.1))-273 and for BAND 11: 1201.14/ln(480.89/(0.0003342*b1+0.1))-273. Effect is similar

      • johndawes_ says:

        alexkostin,
        Sorry to be dense but are the formulas you have outlined in step two for Bands 10 and 11, formulas you could put in the raster calculator? Also is there a way to convert the satellite’s reflectance temp to surface temp? Cool stuff and really wanting to learn as much as possible.
        -J

        • kevin_butler says:

          Yeah, you could definitely do that in raster calculator. The advantage to doing it as a function is the responsiveness. Because it’s not actually writing out a new file (it works in memory) it’s faster to render. Here’s a paper outlining the approach used to calculate land surface temperature for MODIS and AVHRR. Not sure if USGS has something in the works for Landsat as well…
          http://www.geo.utexas.edu/courses/387h/lectures/1-a%20generalized%20split-window%20algorithm%20for%20retrieving%20land-surface%20temperature%20from%20space.pdf

          • johndawes_ says:

            Kevin,
            Thanks for the reply and sorry last question…If I fire on of the equations above into raster calculator, I would assume b1 is is the band for the respective formula correct? So then with the two newly created thermal surfaces, would I just make a composite and look at the stretched values for temperature in F or K?
            Thanks again for your help.
            -J

            John,
            Wordpress is being weird and won’t let me reply to your comment directly. Yes is the answer to your question. B1 is just a place holder for the band in question. You could make a composite or keep them as separate bands.

  2. johndawes_ says:

    Kevin,
    I think I got a good handle on it…I was able to run the equation for bands 10 and 11. The only thing I’m unsure of is that i get a range from -118 C to – 113 C . Does that seem right for at sensor brightness temp or is it too cold?
    -J

  3. poorbaby says:

    I come to step 5. But it’s error in Edit raster function. Input raster error
    you can see error in link:
    http://i.upanh.com/vtcyel
    Can you help me? Thanks

  4. wikikhtn says:

    Please show me How to insert second raster (band 11), I only have first raster (band 10) same as the original default

  5. stepankhach says:

    Hi everyone

    I can’t find “Local function” (I use ArcGIS 10.1). How can I figure out this problem?

    • kevin_butler says:

      It could be that you don’t have the spatial analyst extension, which is required for local function. There is a work around. You can use band arithmetic to make the conversion. Also, at 10.2 there was a patch released which incorporated the layers generated from this workflow into the raster product. If you’re using 10.2 or above, your raster products will have the thermal layers for Landsat.