Surface temperature estimating with Landsat-8 images

The thermal infrared bands of Landsat-8 is permit to estimate land surface temperature. In this case we will calculate the land surface temperature in Kryvyi Rih and nearest areas. For calculation we will use 10th band of Landsat’s image. For this example I use image from July 15th 2016. In that day the sky was very clearly above the Kryvyi Rih. Also in that day the meteorologycal station near the Kryvyi Rih reported that air temperature is 33 Celsium degrees (in time of imaging). Where is practical examples of using similar methodics? For example post “Hell shopping areas” on my blog “DataStory”.

For image processing we will use raster and RStollbox packages. The Kryvyi Rih area polygon I get from OpenStreetMap shapefiles.

library(raster)
library(RStoolbox)
#set a temporary data directory for large graphical objects
rasterOptions(tmpdir = "/home/geka/TEMP_R/")
#! You must have a large free space on your hard disk
#because we will get several large raster objects

metaData <- readMeta("./LC81790272016197LGN00/LC81790272016197LGN00_MTL.txt")
lsat8 <- stackMeta("./LC81790272016197LGN00/LC81790272016197LGN00_MTL.txt")

hazeDN    <- estimateHaze(lsat8, hazeBands = c("B3_dn", "B4_dn"),
                          darkProp = 0.01)

lsat8_sdos <- radCor(lsat8, metaData = metaData,
                     hazeValues = hazeDN,
                     hazeBands = c("B3_dn", "B4_dn"),
                     method = "sdos")

lsat8_sdos

class       : RasterStack 
dimensions  : 7991, 7861, 62817251, 10  (nrow, ncol, ncell, nlayers)
resolution  : 30, 30  (x, y)
extent      : 359085, 594915, 5134785, 5374515  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=36 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
names       :   B1_sre,   B2_sre,   B3_sre,   B4_sre,   B5_sre,   B6_sre,   B7_sre,  B9_sre,   B10_bt,   B11_bt 
min values  :   0.0000,   0.0000,   0.0000,   0.0000,   0.0000,   0.0000,   0.0000,       ?, 147.5171, 141.6706 
max values  :   1.0000,   1.0000,   1.0000,   1.0000,   1.0000,   1.0000,   1.0000,       ?, 368.0307, 366.6090 

KrR_shape <- shapefile("KrR_32636.shp")

plotRGB(lsat8_sdos, r = 4, g = 3, b = 2, stretch = "hist")
plot(KrR_shape, lwd = 2, col = "red", add = TRUE)

lsat8_sdos.KrR <- crop(lsat8_sdos, KrR_shape)

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -InfNAs introduced by coercionNAs introduced by coercion

plotRGB(lsat8_sdos.KrR, r = 4, g = 3, b = 2, stretch = "hist")
plot(KrR_shape, lwd = 2, col = "red", add = TRUE)

For calculating of surface temperature we use formula:

\(LST(°C) = Bt/[1 + (w * Bt/p) * ln(e)] - 273.15\)

Where:

• \(Bt\) is At satellite temperature

• \(w\) is wavelength of emitted radiance

• \(p\) is \(h * c/s\), where \(h\) is Planc’s constant, \(c\) is velocity of light, \(s\) is Boltzmann constant. \(p\) is equal to 14388

• \(e\) is electromagnetic emissivity. For urban areas I use formula from Stathopolou (2007):

\(e = 0.017 * PV + 0.963\)

In formula of electromagnetic emissivity \(PV\) is “Proportion of Vegetation”. For \(PV\) calculation is used the NDVI.

\(PV = [(NDVI - NDVI_{min}) / (NDVI_{max} - NDVI_{min})]^2\)

lsat8.ndvi <- spectralIndices(lsat8_sdos.KrR, red = "B4_sre", nir = "B5_sre", indices = "NDVI")

plot(lsat8.ndvi)

lsat8.ndvi

class       : RasterLayer 
dimensions  : 2029, 1139, 2311031  (nrow, ncol, ncell)
resolution  : 30, 30  (x, y)
extent      : 510465, 544635, 5277495, 5338365  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=36 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : NDVI 
values      : -0.3736169, 0.9435266  (min, max)

lsat8.PV <- ((lsat8.ndvi - (-0.3736169)) / (0.9435266 + (-0.3736169)))^2

plot(lsat8.PV)

lsat8.emissiv <- 0.017 * lsat8.PV + 0.963

lsat8.emissiv

class       : RasterLayer 
dimensions  : 2029, 1139, 2311031  (nrow, ncol, ncell)
resolution  : 30, 30  (x, y)
extent      : 510465, 544635, 5277495, 5338365  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=36 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : NDVI 
values      : 0.963, 1.053804  (min, max)

lsat8.LST <- lsat8_sdos.KrR$B10_bt / (1 + 10.8 * (lsat8_sdos.KrR$B10_bt/14388) * log(lsat8.emissiv)) - 273.15

lsat8.LST

class       : RasterLayer 
dimensions  : 2029, 1139, 2311031  (nrow, ncol, ncell)
resolution  : 30, 30  (x, y)
extent      : 510465, 544635, 5277495, 5338365  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=36 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : layer 
values      : 22.11561, 53.92095  (min, max)

plot(lsat8.LST)

writeRaster(lsat8.LST$layer, "lsat8_LST.sdat", format = "SAGA", overwrite = TRUE)

class       : RasterLayer 
dimensions  : 2029, 1139, 2311031  (nrow, ncol, ncell)
resolution  : 30, 30  (x, y)
extent      : 510465, 544635, 5277495, 5338365  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
data source : /home/geka/Documents/RPubs/UrbanWarm/lsat8_LST.sgrd 
names       : layer 
values      : 22.11561, 53.92095  (min, max)