Algorithm for Automated Mapping of Land Surface Temperature Using LANDSAT 8 Satellite Data

Abstract

Land surface temperature is an important factor in many areas, such as global climate change, hydrological, geo-/biophysical, and urban land use/land cover. As the latest launched satellite from the LANDSAT family, LANDSAT 8 has opened new possibilities for understanding the events on the Earth with remote sensing. This study presents an algorithm for the automatic mapping of land surface temperature from LANDSAT 8 data. The tool was developed using the LANDSAT 8 thermal infrared sensor Band 10 data. Different methods and formulas were used in the algorithm that successfully retrieves the land surface temperature to help us study the thermal environment of the ground surface. To verify the algorithm, the land surface temperature and the near-air temperature were compared. The results showed that, for the first case, the standard deviation was 2.4°C, and for the second case, it was 2.7°C. For future studies, the tool should be refined with in situ measurements of land surface temperature.

Academic Editor:Guiyun Tian

Research Institute of Earth and Space Sciences, Anadolu University, Iki Eylul Campus, 26555 Eskisehir, Turkey

Received 25 November 2015; Revised 21 January 2016; Accepted 4 February 2016

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Land surface temperature (LST) is defined as the temperature felt when the land surface is touched with the hands or the skin temperature of the ground [1]. As one of the most important aspects of the land surface, LST has been a main topic for developing methodologies to be measured from space. LST is an important factor in many areas of studies, such as global climate change, hydrological and agricultural processes, and urban land use/land cover. Calculating LST from remote sensed images is needed since it is an important factor controlling most physical, chemical, and biological processes of the Earth [2]. There is a growing awareness among environmental scientists that remote sensing can and must play a role in providing the data needed to assess ecosystems conditions and to monitor change at all special scales [3]. The tool developed in this paper is simple and does not require any background knowledge so scientists can use it very easy in their researches.

The algorithm introduced in this paper has been developed using ERDAS IMAGINE 2014, with the Model Maker allowing us to create a model that will repeat the process automatically, and it is easy to develop a simple tool useful for doing pixel calculations. Without the tool, the process of retrieving LST is very long, and it is prone to many mistakes. The tool also can be developed in any software supporting pixel calculations from a given image, following step by step this paper. Although an LST retrieval method for LANDSAT 8 has been developed [1, 4], a tool is needed for the complicated process of obtaining the LST. A similar study for retrieving LST in ERDAS IMAGINE has been conducted for LANDSAT 7 data [5] but not for LANDSAT 8. The tool presented in this paper is used for calculating the LST of a given LANDSAT 8 image with the input of the fourth (red wavelength/micrometres, 0.64-0.67), fifth (near infrared (NIR) wavelength/micrometres, 0.85-0.88), and tenth (thermal infrared sensor (TIRS) wavelength/micrometres, 10.60-11.19) bands. Following January 6, 2014, recommendations of USGS of not using TIRS Band 11 due to its larger calibration uncertainty, only Band 10 was included in the algorithm.

2. Data and Methods

The algorithm was created in ERDAS IMAGINE 2014, and it can only be used to process LANDSAT 8 data because of the data complexity. The LST of any Landsat 8 satellite image can be retrieved following the steps of Figure 1. The data of Landsat 8 is available at the Earth Explorer website free of charge. In this study, the TIR band 10 was used to estimate brightness temperature and bands 4 and 5 were used for calculating the NDVI. The metadata of the satellite images used in the algorithm is presented in Table 1.

Table 1: Metadata of the satellite images.

Thermal constant, Band 10
[figure omitted; refer to PDF]	1321.08
[figure omitted; refer to PDF]	777.89
Rescaling factor, Band 10
[figure omitted; refer to PDF]	0.000342
[figure omitted; refer to PDF]	0.1
Correction, Band 10
[figure omitted; refer to PDF]	0.29

Figure 1: Flowchart for LST retrieval.

[figure omitted; refer to PDF]

2.1. Top of Atmospheric Spectral Radiance

The first step of the algorithm is the input of Band 10. After inputting band 10, in the background, the tool uses formulas taken from the USGS web page for retrieving the top of atmospheric (TOA) spectral radiance ( [figure omitted; refer to PDF] ): [figure omitted; refer to PDF] where [figure omitted; refer to PDF] represents the band-specific multiplicative rescaling factor, [figure omitted; refer to PDF] is the Band 10 image, [figure omitted; refer to PDF] is the band-specific additive rescaling factor, and [figure omitted; refer to PDF] is the correction for Band 10 [6].

2.2. Conversion of Radiance to At-Sensor Temperature

After the digital numbers (DNs) are converted to reflection, the TIRS band data should be converted from spectral radiance to brightness temperature (BT) using the thermal constants provided in the metadata file. The following equation is used in the tool's algorithm to convert reflectance to BT [7]: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] and [figure omitted; refer to PDF] stand for the band-specific thermal conversion constants from the metadata.

For obtaining the results in Celsius, the radiant temperature is revised by adding the absolute zero (approx. -273.15°C) [8].

2.3. NDVI Method for Emissivity Correction

2.3.1. Calculating NDVI

Landsat visible and near-infrared bands were used for calculating the Normal Difference Vegetation Index (NDVI). The importance of estimating the NDVI is essential since the amount of vegetation present is an important factor and NDVI can be used to infer general vegetation condition [9]. The calculation of the NDVI is important because, afterward, the proportion of the vegetation ( [figure omitted; refer to PDF] ) should be calculated, and they are highly related with the NDVI, and emissivity ( [figure omitted; refer to PDF] ) should be calculated, which is related to the [figure omitted; refer to PDF] : [figure omitted; refer to PDF] where NIR represents the near-infrared band (Band 5) and [figure omitted; refer to PDF] represents the red band (Band 4).

2.3.2. Calculating the Proportion of Vegetation

[figure omitted; refer to PDF] is calculated according to (4). A method for calculating [figure omitted; refer to PDF] [4] suggests using the NDVI values for vegetation and soil ( [figure omitted; refer to PDF] and [figure omitted; refer to PDF] ) to apply in global conditions [10]: [figure omitted; refer to PDF] However, since the NDVI values differ for every area, the value for vegetated surfaces, 0.5, may be too low. Global values from NDVI can be calculated from at-surface reflectivities, but it would not be possible to establish global values in the case of an NDVI computed from TOA reflectivities, since [figure omitted; refer to PDF] and [figure omitted; refer to PDF] will depend on the atmospheric conditions [11].

2.3.3. Calculating Land Surface Emissivity

The land surface emissivity (LSE ( [figure omitted; refer to PDF] )) must be known in order to estimate LST, since the LSE is a proportionality factor that scales blackbody radiance (Planck's law) to predict emitted radiance, and it is the efficiency of transmitting thermal energy across the surface into the atmosphere [12]. The determination of the ground emissivity is calculated conditionally as suggested in [10]: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are the vegetation and soil emissivities, respectively, and [figure omitted; refer to PDF] represents the surface roughness ( [figure omitted; refer to PDF] = 0 for homogenous and flat surfaces) taken as a constant value of 0.005 [13]. The condition can be represented with the following formula and the emissivity constant values shown in Table 1 [4]: [figure omitted; refer to PDF] When the NDVI is less than 0, it is classified as water, and the emissivity value of 0.991 is assigned. For NDVI values between 0 and 0.2, it is considered that the land is covered with soil, and the emissivity value of 0.996 is assigned. Values between 0.2 and 0.5 are considered mixtures of soil and vegetation cover and (6) is applied to retrieve the emissivity. In the last case, when the NDVI value is greater than 0.5, it is considered to be covered with vegetation, and the value of 0.973 is assigned.

The last step of retrieving the LST or the emissivity-corrected land surface temperature [figure omitted; refer to PDF] is computed as follows [14]: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the LST in Celsius (°C, (2)), BT is at-sensor BT (°C), [figure omitted; refer to PDF] is the wavelength of emitted radiance (for which the peak response and the average of the limiting wavelength ( [figure omitted; refer to PDF] ) [15] will be used), [figure omitted; refer to PDF] is the emissivity calculated in (6), and [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the Boltzmann constant (1.38 × 10^-23 J/K), [figure omitted; refer to PDF] is Planck's constant (6.626 × 10^-34 J s), and [figure omitted; refer to PDF] is the velocity of light (2.998 × 10⁸ m/s) [9].

3. LST Validation

The two major LST validation models are through ground measurements or near-surface air temperature [16, 17]. The LST results comparing with the ground measurements results may have an error up to 5°C; in the case of Srivastava et al., the accuracy of the results in some area showed difference of ±2°C with actual ground temperature measurements. According to Liu and Zhang, another method using the mean near-surface air temperature to verify the retrieved LST results showed that the LST retrieving error is about 0.7°C. For the validation, six representative points have been used.

For the validation of the final retrieved LST results in the presented tool, the mean near-surface air temperature was used [18] but with bigger amount of data and taking not only the mean temperature but also the actual temperature in the given pixel at the moment of the satellite passing over the area for 27 representative points.

The comparison was made with air temperature, which is different and can sometimes result in big differences since the resolution of LANDSAT 8 for the used bands is 100 m for the thermal band and 30 m for the red and NIR bands. The LST was calculated and taken for the pixel in which the meteorological station fell. Sometimes, the differences can be very big depending on the weather condition and other factors [19]. It should also be taken into consideration that there is 1.1 to 2 meters' difference between the LST and the air temperature, which means that differences in the temperatures are normal and expected.

3.1. Application of the Algorithm to Ontario and Quebec, Canada

Hourly data were collected from the Canadian Weather and Meteorology website (http://climate.weather.gc.ca/) and used for comparison with the retrieved LST for which, according to the available data, satellite images were downloaded for 02/05/2015 (Toronto area) and 04/06/2015 (Moncton area) for the areas shown in Figure 2.

Figure 2: Application of algorithm in Ontario and Quebec, Canada. (a) Geographic location of Ontario in Canada; (b) frames of satellite images of study areas; (c) first case located between Toronto and Huntsville; (d) second case located in surrounding area of the city of Moncton.

[figure omitted; refer to PDF]

The study area was the Canadian provinces of Ontario and Quebec (Figure 2). One satellite image was downloaded from each of the two provinces. These areas were chosen because of their specifications. That is, both study areas included water, urban areas, and green areas.

3.2. Comparison of LST Validation Results

To compare the results, two different satellite images from two different dates in two different areas were chosen according to the available data. After downloading the satellite images from http://earthexplorer.usgs.gov/, LSTs were retrieved in ERDAS using the algorithm presented in this paper. In the first case, the satellite image was located between Toronto and the city of Huntsville near Lake Simcoe in Ontario, Canada. For this area, 20 meteorological stations were found, but only 16 of them were used for the accuracy assessment because of the presence of clouds or other unwanted events. The differences between the retrieved LSTs and the air temperatures and details on the stations are presented in Table 2 and Figure 3.

Table 2: Emissivity of representative terrestrial materials for LANDSAT 8 TIRS Band 10.

Terrestrial material	Water	Building	Soil	Vegetation
Emissivity	0.991	0.962	0.966	0.973

Figure 3: Retrieved LST image and meteorological stations from first study area used in accuracy assessment.

[figure omitted; refer to PDF]

In the second case, the study area was located in the area surrounding the city of Moncton and included part of New Brunswick, Prince Edward Island, and Nova Scotia in Canada. For this area, we found 11 meteorological stations, and all of them were used for the accuracy assessment. The details are presented in Table 3 and Figure 4.

Table 3: Details and differences of station from first study case.

ST name	Data 12pm	LST	Difference	Latitude	Longitude	[figure omitted; refer to PDF]
Barrie-Oro	19.9	20.9	-1.0	44°29[variant prime]00.000[variant prime][variant prime]	79°33[variant prime]00.000[variant prime][variant prime]	289.00 m
Pa Hardwood Mountain Bike Park	19.2	24.9	-5.7	44°31[variant prime]08.900[variant prime][variant prime]	79°35[variant prime]24.200[variant prime][variant prime]	334.50 m
Borden Awos	20.0	20.7	-0.7	44°16[variant prime]20.000[variant prime][variant prime]	79°54[variant prime]42.000[variant prime][variant prime]	222.50 m
Pa Udora Strong	19.8	22.3	-2.5	44°15[variant prime]03.600[variant prime][variant prime]	79°12[variant prime]18.300[variant prime][variant prime]	266.50 m
Lagoon City	10.7	9.2	1.5	44°32[variant prime]50.000[variant prime][variant prime]	79°13[variant prime]00.000[variant prime][variant prime]	220.70 m
Collingwood	13.8	18.7	-4.9	44°30[variant prime]00.000[variant prime][variant prime]	80°13[variant prime]00.000[variant prime][variant prime]	179.80 m
Mono Centre	19.1	18.4	0.7	44°01[variant prime]56.100[variant prime][variant prime]	80°01[variant prime]28.010[variant prime][variant prime]	436.00 m
Uxbridge West	19.4	22.3	-2.9	44°05[variant prime]54.000[variant prime][variant prime]	79°09[variant prime]49.020[variant prime][variant prime]	325.00 m
Pa Uxbridge Taris	19.7	23.2	-3.5	44°03[variant prime]16.000[variant prime][variant prime]	79°06[variant prime]55.000[variant prime][variant prime]	359.50 m
Pa Atmos Vaughan	20.9	17.9	3.1	43°51[variant prime]47.700[variant prime][variant prime]	79°32[variant prime]28.900[variant prime][variant prime]	254.00 m
Pa Angus Glen Golf Club	20.4	23.3	-2.9	43°54[variant prime]29.800[variant prime][variant prime]	79°19[variant prime]23.400[variant prime][variant prime]	230.50 m
Toronto Buttonville A	21.2	23.3	-2.1	43°51[variant prime]44.000[variant prime][variant prime]	79°22[variant prime]12.000[variant prime][variant prime]	198.10 m
Pa Claremont Silo Farm	19.8	21.8	-2.0	43°59[variant prime]37.900[variant prime][variant prime]	79°05[variant prime]43.900[variant prime][variant prime]	263.50 m
Pa Markham North Toyota	21.6	27.4	-5.8	43°49[variant prime]01.000[variant prime][variant prime]	79°20[variant prime]33.810[variant prime][variant prime]	187.50 m
Pa Atmos Claremont	21.5	23.5	-2.0	43°56[variant prime]09.800[variant prime][variant prime]	79°05[variant prime]05.400[variant prime][variant prime]	167.00 m
Pa Atmos Erin	19.8	21.6	-1.8	43°49[variant prime]33.620[variant prime][variant prime]	80°07[variant prime]12.840[variant prime][variant prime]	470.00 m

Figure 4: Retrieved LST image and meteorological stations from second study area used in accuracy assessment.

[figure omitted; refer to PDF]

4. Conclusion

This paper presented a new LST software tool and its algorithm created in ERDAS for calculating the LST from LANDSAT 8 TIRS. The algorithm was derived using the observed thermal radiance of the TIRS Band 10 of LANDSAT 8 TIRS. To verify the final retrieved LST results, the near-surface air temperature method was used. From the analysis of the two areas in Canada from two different dates, the standard deviation calculated for the first case based on 16 meteorological stations was 2.4°C, and that for the second case based on 11 stations was 2.7°C. It should be mentioned that sometimes, the difference between the near-surface temperature and the LST can be drastic since we are comparing two different temperatures in different places (ground temperature and 1.1 to 2.0 m off the ground). It should also be taken into consideration that the resolution of the LANDSAT 8 TIRS data is 100 m for the thermal band and 30 m for the red and NIR bands. Values smaller than -5°C in the two cases were considered to be clouds or other unwanted events on the satellite images since the data were from springtime; it was not expected. From Tables 3 and 4 and Figures 3 and 4, it can be concluded that, for the first case, the smallest difference between the LST retrieved from the presented tool and the near-air temperature was 0.7°C and the biggest was 5.8°C. In the second study case, the smallest difference was 0.2°C and the biggest was 7.8°C.

Table 4: Details and differences of station from second study case.

ST name	Data 11pm	LST	Difference	Latitude	Longitude	[figure omitted; refer to PDF]
Mechanic Settlement	8.4	8.6	-0.2	45°41[variant prime]37.040[variant prime][variant prime]	65°09[variant prime]54.040[variant prime][variant prime]	403.00 m
Moncton Intl A	15.3	14.3	1.0	46°06[variant prime]44.000[variant prime][variant prime]	64°40[variant prime]43.000[variant prime][variant prime]	70.70 m
Fundy Park (Alma) Cs	12.4	10.2	2.2	45°36[variant prime]00.000[variant prime][variant prime]	64°57[variant prime]00.000[variant prime][variant prime]	42.70 m
Buctouche Cda Cs	15.1	14.9	0.2	46°25[variant prime]49.006[variant prime][variant prime]	64°46[variant prime]05.009[variant prime][variant prime]	35.90 m
Nappan Auto	15.4	13.3	2.1	45°45[variant prime]34.400[variant prime][variant prime]	64°14[variant prime]29.200[variant prime][variant prime]	19.80 m
Doaktown Auto Rcs	15.9	11.5	4.4	46°35[variant prime]06.090[variant prime][variant prime]	66°00[variant prime]35.071[variant prime][variant prime]	43.00 m
Saint John A	15.1	12.5	2.6	45°18[variant prime]58.000[variant prime][variant prime]	65°53[variant prime]24.000[variant prime][variant prime]	108.80 m
Parrsboro	15.9	8.1	7.8	45°24[variant prime]48.000[variant prime][variant prime]	64°20[variant prime]49.000[variant prime][variant prime]	30.90 m
Gagetown A	15.0	10.2	4.8	45°50[variant prime]00.000[variant prime][variant prime]	66°26[variant prime]00.000[variant prime][variant prime]	50.60 m
Gagetown Awos A	16.3	18.6	-2.3	45°50[variant prime]20.000[variant prime][variant prime]	66°26[variant prime]59.000[variant prime][variant prime]	/
Summerside	13.9	11.7	2.2	46°26[variant prime]28.000[variant prime][variant prime]	63°50[variant prime]17.000[variant prime][variant prime]	12.20 m

The tool can be significant since it is opening an opportunity to many researchers to get to the LST values easy and they can apply them in a number of researchs. The tool presented in this paper produced quite good results considering that the accuracy assessment was made with near-air temperatures and it can be used for different types of research. For future studies, the tool should be refined with in situ measurements of LST.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

References

[1] A. Rajeshwari, N. Mani, "Estimation of land surface temperature of dindigul district using landsat 8 data," International Journal of Research in Engineering and Technology , vol. 3, no. 5, pp. 122-126, 2014.

[2] F. Becker, Z.-L. Li, "Towards a local split window method over land surfaces," International Journal of Remote Sensing , vol. 11, no. 3, pp. 369-393, 1990.

[3] S. Ustin Manual of Remote Sensing: Remote Sensing for Natural Resource Management and Environmental Monitoring , John Wiley & Sons, Hoboken, NJ, USA, 2004.

[4] F. Wang, Z. Qin, C. Song, L. Tu, A. Karnieli, S. Zhao, "An improved mono-window algorithm for land surface temperature retrieval from landsat 8 thermal infrared sensor data," Remote Sensing , vol. 7, no. 4, pp. 4268-4289, 2015.

[5] Q. Q. Sun, J. J. Tan, Y. H. Xu, "An ERDAS image processing method for retrieving LST and describing urban heat evolution: a case study in the Pearl River Delta Region in South China," Environmental Earth Sciences , vol. 59, no. 5, pp. 1047-1055, 2009.

[6] J. A. Barsi, J. R. Schott, S. J. Hook, N. G. Raqueno, B. L. Markham, R. G. Radocinski, "Landsat-8 thermal infrared sensor (TIRS) vicarious radiometric calibration," Remote Sensing , vol. 6, no. 11, pp. 11607-11626, 2014.

[7] USGS 2013, http://landsat.usgs.gov/Landsat8_Using_Product.php

[8] H.-Q. Xu, B.-Q. Chen, "Remote sensing of the urban heat island and its changes in Xiamen City of SE China," Journal of Environmental Sciences , vol. 16, no. 2, pp. 276-281, 2004.

[9] Q. H. Weng, D. S. Lu, J. Schubring, "Estimation of land surface temperature-vegetation abundance relationship for urban heat island studies," Remote Sensing of Environment , vol. 89, no. 4, pp. 467-483, 2004.

[10] J. A. Sobrino, J. C. Jimenez-Muñoz, L. Paolini, "Land surface temperature retrieval from LANDSAT TM 5," Remote Sensing of Environment , vol. 90, no. 4, pp. 434-440, 2004.

[11] J. C. Jimenez-Muñoz, J. A. Sobrino, A. Plaza, L. Guanter, J. Moreno, P. Martinez, "Comparison between fractional vegetation cover retrievals from vegetation indices and spectral mixture analysis: case study of PROBA/CHRIS data over an agricultural area," Sensors , vol. 9, no. 2, pp. 768-793, 2009.

[12] J. C. Jimenez-Munoz, J. A. Sobrino, A. Gillespie, D. Sabol, W. T. Gustafson, "Improved land surface emissivities over agricultural areas using ASTER NDVI," Remote Sensing of Environment , vol. 103, no. 4, pp. 474-487, 2006.

[13] J. A. Sobrino, N. Raissouni, "Toward remote sensing methods for land cover dynamic monitoring: application to Morocco," International Journal of Remote Sensing , vol. 21, no. 2, pp. 353-366, 2000.

[14] M. Stathopoulou, C. Cartalis, "Daytime urban heat islands from Landsat ETM+ and Corine land cover data: an application to major cities in Greece," Solar Energy , vol. 81, no. 3, pp. 358-368, 2007.

[15] B. L. Markham, J. L. Barker, "Spectral characterization of the Landsat Thematic Mapper sensors," International Journal of Remote Sensing , vol. 6, no. 5, pp. 697-716, 1985.

[16] Z.-L. Li, B.-H. Tang, H. Wu, H. Ren, G. Yan, Z. Wan, I. F. Trigo, J. A. Sobrino, "Satellite-derived land surface temperature: current status and perspectives," Remote Sensing of Environment , vol. 131, pp. 14-37, 2013.

[17] P. K. Srivastava, T. J. Majumdar, A. K. Bhattacharya, "Surface temperature estimation in Singhbhum Shear Zone of India using Landsat-7 ETM+ thermal infrared data," Advances in Space Research , vol. 43, no. 10, pp. 1563-1574, 2009.

[18] L. Liu, Y. Z. Zhang, "Urban heat island analysis using the landsat TM data and ASTER data: a case study in Hong Kong," Remote Sensing , vol. 3, no. 7, pp. 1535-1552, 2011.

[19] K. Gallo, R. Hale, D. Tarpley, Y. Yu, "Evaluation of the relationship between air and land surface temperature under clear- and cloudy-sky conditions," Journal of Applied Meteorology and Climatology , vol. 50, no. 3, pp. 767-775, 2011.