Converts Landsat ETM+ thermal imagery to degrees Celsius. See associated blog post here: http://www.hakimabdi.com/20111030/estimating-land-surface-temperature-from-landsat-thermal-imagery/

ManhattanLST.r

#########################################################################################
# title         : ManhattanLST.R
# purpose       : Converts Landsat ETM+ thermal imagery to degrees Celsius
# author        : Abdulhakim Abdi (@HakimAbdi)
# input         : Landsat TM / ETM+ band 6
# output        : Land surface temperature map in degrees Celsius
#########################################################################################

# Install and load the required packages
install.packages(c("rgdal" "sp", "raster"), repos='http://cran.r-project.org')
library(rgdal, sp, raster)

#1 Set the working directory
setwd("g:/data")

#2 Calculate NDVI
b3 = raster("Landsat_band_3.img")
b4 = raster("Landsat_band_4.img")
ndvi = (b4-b3)/(b4+b3)

#3 Calculate emissivity as per Sobrino et. al. (2004)
pv = (ndvi - 0.2)/(0.5 - 0.2)
em = (0.004*pv) + 0.986

# Note: The following step (#4) is only necessary if you are working 
# with Landsat imagery acquired before February 25th, 2010. All
# Landsat data processed after that date will have a thermal band
# resolution of 30 meters. See the announcement here: 
# http://landsat.usgs.gov/about_LU_Vol_4_Issue_Special_Edition.php

#4 Import Landsat thermal band and resample to 30m
thermal = raster("Landsat_band_61.img")
b61 = resample(thermal, b4, method="bilinear")

#5 Calculate radiance and black body temperature as per NASA (2009)
# Note: See Chander et. al. (2009) for applicable values.

rad = (17.04/254)*(b61-1)
Tb= 1282.71/log((666.09/rad)+1)

#6 Compute temperature as per Weng et. al. (2004) and convert to Celsius
# Note: See Chander et. al. (2009) for applicable values.
kelvin = Tb/(1+(11.5*(Tb/0.01438))*log(em))
celsius = kelvin-273

#7 Write to file
writeRaster(celsius, filename="celsius.tif", format="GTiff")

# End of script

#########################################################################################
#########################################################################################
#########################################################################################

# The "landsat" package's thermalband function summarizes all the above steps:

function (x, band)
{
    if (band == 6)
        band.coefs <- c(0.055376, 1.18, 607.76, 1260.56)
    if (band == 61)
        band.coefs <- c(0.067087, -0.07, 666.09, 1282.71)
    if (band == 62)
        band.coefs <- c(0.037205, 3.16, 666.09, 1282.7)
    results <- x
    x <- as.vector(as.matrix(x))
    x <- x * band.coefs[1] + band.coefs[2]
    x <- band.coefs[4]/log(band.coefs[3]/x + 1)
    if (class(results) == "SpatialGridDataFrame")
        results@data[, 1] <- x
    else if (is.data.frame(results))
        results <- data.frame(matrix(x, nrow = nrow(results),
            ncol = ncol(results)))
    else if (is.matrix(results))
        results <- matrix(x, nrow = nrow(results), ncol = ncol(results))
    else results <- x
    results
}