Nothing Special   »   [go: up one dir, main page]
Next Article in Journal
Evaluation of MODIS Gross Primary Production across Multiple Biomes in China Using Eddy Covariance Flux Data
Previous Article in Journal
Long Term Validation of Land Surface Temperature Retrieved from MSG/SEVIRI with Continuous in-Situ Measurements in Africa
You seem to have javascript disabled. Please note that many of the page functionalities won't work as expected without javascript enabled.
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Remote Sensing
Volume 8
Issue 5
10.3390/rs8050413
remotesensing-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Application of Open Source Coding Technologies in the Production of Land Surface Temperature (LST) Maps from Landsat: A PyQGIS Plugin

by
Milton Isaya Ndossi
* and
Ugur Avdan
Institute of Space and Earth Sciences, Anadolu University, Iki Eylul Campus, Eskisehir 26555, Turkey
*
Author to whom correspondence should be addressed.
Remote Sens. 2016, 8(5), 413; https://doi.org/10.3390/rs8050413
Submission received: 29 March 2016 / Revised: 13 April 2016 / Accepted: 9 May 2016 / Published: 13 May 2016
Browse Figures
Graphical abstract
"> Figure 1
<p>Application’s Interface for calculating radiance.</p> "> Figure 2
<p>Application’s Interface for calculating brightness temperature.</p> "> Figure 3
<p>Application’s interface for Land Surface Emissivity (LSE) estimation based on the NDVI image.</p> "> Figure 4
<p>Application’s interface for the Planck Equation.</p> "> Figure 5
<p>Application’s interface for Mono-Window algorithm.</p> "> Figure 6
<p>Application’s interface for the RTE.</p> "> Figure 7
<p>Application’s interface for the SCA.</p> "> Figure 8
<p>LST maps derived from Landsat 5 TM: (<b>a</b>) LST extracted using the SCA; (<b>b</b>) LST extracted using the RTE; (<b>c</b>) LST extracted using the MWA; and (<b>d</b>) LST extracted using the Planck equation.</p> "> Figure 8 Cont.
<p>LST maps derived from Landsat 5 TM: (<b>a</b>) LST extracted using the SCA; (<b>b</b>) LST extracted using the RTE; (<b>c</b>) LST extracted using the MWA; and (<b>d</b>) LST extracted using the Planck equation.</p> "> Figure 8 Cont.
<p>LST maps derived from Landsat 5 TM: (<b>a</b>) LST extracted using the SCA; (<b>b</b>) LST extracted using the RTE; (<b>c</b>) LST extracted using the MWA; and (<b>d</b>) LST extracted using the Planck equation.</p> "> Figure 9
<p>LST maps derived from Landsat 7 ETM+: (<b>a</b>) LST extracted using the SCA; (<b>b</b>) LST extracted using the RTE; (<b>c</b>) LST extracted using the MWA; and (<b>d</b>) LST extracted using the Planck equation.</p> "> Figure 9 Cont.
<p>LST maps derived from Landsat 7 ETM+: (<b>a</b>) LST extracted using the SCA; (<b>b</b>) LST extracted using the RTE; (<b>c</b>) LST extracted using the MWA; and (<b>d</b>) LST extracted using the Planck equation.</p> "> Figure 10
<p>LST maps derived from Landsat 7 ETM+: (<b>a</b>) LST extracted using the SCA; (<b>b</b>) LST extracted using the RTE; (<b>c</b>) LST extracted using the MWA; and (<b>d</b>) LST extracted using the Planck equation.</p> "> Figure 10 Cont.
<p>LST maps derived from Landsat 7 ETM+: (<b>a</b>) LST extracted using the SCA; (<b>b</b>) LST extracted using the RTE; (<b>c</b>) LST extracted using the MWA; and (<b>d</b>) LST extracted using the Planck equation.</p> "> Figure 11
<p>LST maps derived from Landsat 7 ETM+: (<b>a</b>) LST extracted using the SCA; (<b>b</b>) LST extracted using the RTE; (<b>c</b>) LST extracted using the MWA; and (<b>d</b>) LST extracted using the Planck equation.</p> "> Figure 11 Cont.
<p>LST maps derived from Landsat 7 ETM+: (<b>a</b>) LST extracted using the SCA; (<b>b</b>) LST extracted using the RTE; (<b>c</b>) LST extracted using the MWA; and (<b>d</b>) LST extracted using the Planck equation.</p> "> Figure 12
<p>LST maps derived from Landsat 8 TIRS: (<b>a</b>) LST extracted using the SCA; (<b>b</b>) LST extracted using the RTE; (<b>c</b>) LST extracted using the MWA; and (<b>d</b>) LST extracted using the Planck equation.</p> "> Figure 12 Cont.
<p>LST maps derived from Landsat 8 TIRS: (<b>a</b>) LST extracted using the SCA; (<b>b</b>) LST extracted using the RTE; (<b>c</b>) LST extracted using the MWA; and (<b>d</b>) LST extracted using the Planck equation.</p> "> Figure 13
<p>LST maps derived from Landsat 8 TIRS: (<b>a</b>) LST extracted using the SCA; (<b>b</b>) LST extracted using the RTE; (<b>c</b>) LST extracted using the MWA; and (<b>d</b>) LST extracted using the Planck equation.</p> "> Figure 13 Cont.
<p>LST maps derived from Landsat 8 TIRS: (<b>a</b>) LST extracted using the SCA; (<b>b</b>) LST extracted using the RTE; (<b>c</b>) LST extracted using the MWA; and (<b>d</b>) LST extracted using the Planck equation.</p> ">
Versions Notes

Abstract

:
This paper presents a Python QGIS (PyQGIS) plugin, which has been developed for the purpose of producing Land Surface Temperature (LST) maps from Landsat 5 TM, Landsat 7 ETM+ and Landsat 8 TIRS, Thermal Infrared (TIR) imagery. The plugin has been developed purposely to ease the process of LST extraction from Landsat Visible, Near Infrared (VNIR) and TIR imagery. It has the ability to estimate Land Surface Emissivity (LSE), calculating at-sensor radiance, calculating brightness temperature and performing correction of brightness temperature against atmospheric interference though the Plank function, Mono Window Algorithm (MWA), Single Channel Algorithm (SCA) and the Radiative Transfer Equation (RTE). Using the plugin, LST maps of Moncton, New Brunswick, Canada have been produced for Landsat 5 TM, Landsat 7 ETM+ and Landsat 8 TIRS. The study put much more emphasis on the examination of LST derived from the different algorithms of LST extraction from VNIR and TIR satellite imagery. In this study, the best LST values derived from Landsat 5 TM were obtained from the RTE and the Planck function with RMSE of 2.64 °C and 1.58 °C, respectively. While the RTE and the Planck function produced the best results for Landsat 7 ETM+ with RMSE of 3.75 °C and 3.58 °C respectively and for Landsat 8 TIRS LST retrieval, the best results were obtained from the Planck function and the SCA with RMSE of 2.07 °C and 3.06 °C, respectively.

Graphical Abstract">

Graphical Abstract

1. Introduction

Land Surface Temperature (LST) is the temperature of the surface of the ground [1]. It is one of the most vital data recorded by satellites in the recent decades. It is widely used in a variety of fields including but not limited to evapotranspiration, climate change, hydrological cycle, vegetation monitoring, urban climate and environmental studies [2]. LST is an important component of the Earth’s heat budget and is commonly required in applications of hydrology, meteorology and climatology [3].
As a result of the complexities of surface temperature over land, ground measurements cannot practically provide temperature values over wide areas. With the development of satellite technologies and the availability of satellite imagery with a high spatial resolution, satellite data remains to be the only way that can be used to measure LST over the entire globe with sufficiently high spatial resolution and with complete spatially averaged rather than point values [2].
Modern Landsat space crafts have been mounted with thermal infrared sensors which have the ability to acquire thermal images and as a result providing users with large volumes of temperature data. These sensors work basing on the physics principle that everybody on Earth radiates infrared radiation and the radiation is proportional to the temperature of the body.
Many algorithms have been developed for the purpose of extracting LST from satellite imagery, some of them include the Split Window Algorithm (SWA) [4], Single Channel Algorithm (SCA) [5], Mono-Window Algorithm (MWA) [6], and the Radiative Transfer Equation (RTE) [2]. As a result of difficulties arising from the complexities of implementing these algorithms, many users have failed to benefit enough from these data. Not only that but also due to the unavailability of these algorithms in the commonly used Remote Sensing (RS) and Geographic Information Systems (GIS) software. Most of the RS and GIS software available in the market are also very expensive to acquire and do not provide researchers with their source codes which can also play a great role in algorithm improvement and research arguments.
In this study, the Mono-Window Algorithm (MWA), Single Channel Algorithm (SCA), the Radiative Transfer Equation (RTE), the Plank Equation have been implemented to work as a plugin in QGIS. In addition to the algorithms, the plugin has an ability to estimate Land Surface Emissivity (LSE), calculating the Normalized Difference Vegetation Index (NDVI), calculating Top of Atmosphere (TOA) radiance and brightness temperature.
A total of six scenes; two for each sensor have been used in the study. These scenes cover Moncton, New Brunswick, Canada. The scenes used in the study cover different seasons in order to enable the examination of how the algorithms behave in different seasons of the year. The RTE, MWA, Planck function and the SCA have been used to derive the LST of the study area. To perform accuracy assessment of the LST obtained from the algorithms and meteorological data obtained from the Canadian weather and meteorology website which was measured in hourly basis have been used.
The developed plugin is free to download from the official QGIS repository [7] in order for the users of the Landsat sensors to make use of it. As a result of the code being available to be viewed, modified and distributed freely, it is expected that it will further promote debates, discussions, application and improvement of the algorithms.

2. Materials and Methods

2.1. Data

2.1.1. Landsat

Landsat is the longest series of Earth Observation Satellites (EOS) in use today. These sensors have been monitoring the Earth for more than four decades providing a continuity of data though out their lifetime. The first series of these satellites was launched in 1972 and it was named as the Earth Resources Technology Satellite; it was later on renamed to Landsat 1. Since then, there has been a total of 8 Landsat satellites. Out of the eight, Landsat 6 failed to attain orbit and fell down to Earth in 1993. The remaining satellites have proven to be successful and have been providing researchers with massive volumes of data which have been used in many studies. Landsat sensors are equipped with instruments which allow them to detect electromagnetic radiation ranging from the visible to thermal infrared region of the electromagnetic spectrum. Landsat data is freely available to download from the United States Geological Survey (USGS) website [8]. Table 1 shows the Landsat scenes used in the study.

2.1.2. Meteorological Data

Historical data were used to model the differences between the LST obtained from the MWA, SCA, RTE and the Planck equation. The meteorological data were obtained from the Canadian weather and meteorology website [9] in hourly basis. As a result of the effects of day light saving, one hour was added to the local time whenever daylight saving was observed (according to environment Canada, one hour should be added whenever day light saving is observed). Day light saving was observed in the scenes belonging to the months of November, December and February, and as a result, the meteorological data for 12:00 p.m. was used. The meteorological data for the scenes belonging to the months of June and May were compared to meteorological data measured at 11:00 a.m., as daylight saving was not observed in these months.

2.1.3. Plugin Development Tools

The plugin has been developed using the Python programming language. It has been chosen because it is a language which is supported by the QGIS Application Programming Interface (API), it is independent of platform and it is supported by a variety of open source geospatial libraries. To further enable raster processing, the Geospatial Data Abstraction Library (GDAL) and the PyQt4 framework have been used to create the open source graphical based QGIS plugin.

2.2. Methodology

2.2.1. Conversion of Digital Numbers to at-Sensor Radiance

The thermal data in satellite imagery of Landsat sensors are stored in Digital Numbers (DN). DNs are used as a way of representing pixels which have not yet been calibrated into meaningful units. They are a representation of different levels of radiance in a raster image. After the satellite images have been obtained, the first process is to convert the Digital Numbers (DN) to radiance. Equation (1) shows the equation used to convert DN to spectral radiance [10] in the Landsat 8 TIRS sensor.
Lλ = MLQcal + AL − Oi
where Lλ is the Top-of-Atmosphere (TOA) spectral radiance in W/(m2·sr·μm), ML is the band specific multiplicative rescaling factor from the metadata, AL is the band-specific additive rescaling factor from the metadata, Qcal is the quantized and calibrated standard product pixel values (DN), and Oi is the offsets issued by USGS for the calibration of the TIRS bands. The Oi offset calibration has been obtained from the USGS. A value of −0.29 has been used as a calibration factor for Landsat 8 TIRS band 10 [11]. This value is subtracted from the radiances in order to calibrate the radiances recorded in Landsat 8’s band 10. Although Landsat 8 TIRS has two Thermal Infrared (TIR) bands, only data from band 10 are suitable for use in LST retrieval at the moment due to uncertainty in the band 11 values [11,12].
For Landsat 5 TM and Landsat 7 ETM+, to convert DNs to radiance, two equations are available, which are the “gain and bias” method and the spectral radiance scaling method. However, the spectral radiance equation is just another way of representing the gain and bias equation. The “gain and bias equation” and the spectral radiance rescaling equations are shown in Equations (2) and (3), respectively [13].
Lλ = gain × Qcal + bias
Lλ = ((LMAXλ − LMINλ)/(QCALMAX − QCALMIN)) * (QCAL − QCALMIN) + LMINλ
where Lλ is the spectral radiance at the sensor’s aperture in W/(m2·sr·μm), gain is the rescaled gain (the data product “gain” contained in the Level 1 product header or ancillary data record) in W/(m2·sr·μm), bias is the rescaled bias in W/(m2·sr·μm), QCAL is the quantized calibrated pixel value in DN. It corresponds to the thermal bands of the sensors i.e., Landsat TM/ETM+ band 6. LMINλ is the spectral radiance that is scaled to QCALMIN in W/(m2·sr·μm), LMAXλ is the spectral radiance that is scaled to QCALMAX in W/(m2·sr·μm), QCALMIN is the minimum quantized calibrated pixel value and QCALMAX is the maximum quantized calibrated pixel value (corresponding to LMAXλ) in DN = 255. With an exception of QCAL, the rest of the variables are provided in the metadata file of a scene.
The gain and bias rescaling factors of Landsat 7 ETM+ scenes must come from the header file since they change for every scene, however, for Landsat 5 TM scenes, the gain and bias rescaling factors are constant [14]. Figure 1 shows the interface of converting DN values to radiance. The screen is accessible through the “Radiance” tab of the plugin. The interface allows the user to browse for the metadata file, the thermal infrared band and specify the output location of the calculated radiance.

2.2.2. Computation of Brightness Temperature

Brightness temperature is the temperature that is needed by a blackbody in order for it to be able to emit the same amount of radiation per unit of surface area as the body being observed [15]. When the radiance of an object equals to that of a blackbody, the physical temperature of this blackbody is known as the brightness temperature of the object. Brightness temperature has the ability to represent temperature measurements but lacks the physical meaning of temperature [16]. After the DNs have been converted to radiance, the next step is to convert the radiance to brightness temperature.
In order to convert the radiance to brightness temperature, Equation (4) has been used in the study [10].
T s e n =   K 2 l n ( K 1 L λ + 1 )
In Equation (4), Tsen stands for the at-sensor brightness temperature (K), Lλ is the TOA spectral radiance, K1 is the band specific thermal conversion constant from the metadata and K2 band specific thermal conversion constant from the metadata. The same equation is used for all sensors. The K1 and K2 values may vary depending on the sensor and the wavelengths by which the thermal bands operate. The values of K1 and K2 can be obtained from the metadata file of a scene. Figure 2 shows the interface of computing the brightness temperature. The screen is accessible through the “Brightness Temperature” tab of the plugin.

2.2.3. Determination of Land Surface Emissivity (LSE)

Emissivity ε is the ratio of the emitted radiance ελ at wavelength λ and the blackbody emission Bλ(T) at wavelength λ and temperature T [17]. It is the ratio which compares the radiating capacity of a surface to that of a blackbody [15]. In the real world, a material which satisfies the properties of a perfect blackbody does not exist. As a result of this, there is a need of doing LSE correction when deriving LST from space. In order to be able to associate the temperature of thermal infrared energy radiated by a given object, it is necessary to know the emissivity of that object.
The emissivity of a substance is determined by its thermo-physical properties [18]. The composition of a surface is the determinant of the surface’s emissivity. The emissivity of surfaces is as well variable with wavelength. Through the use of Normalized Difference Vegetation Index (NDVI) it is possible to estimate LSE of different terrestrial materials in the 10–12 µm wavelength range [18]. The composition of the land surface changes over time, it is therefore important to estimate the LSE of an area during the satellite overpass time. To calculate the NDVI of the surface, Equation (5) has been used in the study [19]. Atmospheric parameters of scattering and absorption may also affect the estimation of land surface emissivity from NDVI [20]. In this study, the effects of scattering and absorption to NDVI were not taken into consideration.
N D V I =   N I R − R N I R + R
In Equation (5), NIR stands for the Near Infrared band and R is the Red band. To estimate the LSE using the Landsat sensors, two NDVI based algorithms have been implemented in the application.
  • Algorithm based on the NDVI image
According to Zhang, Wang et al. [21], when the NDVI of an area is known, the LSE can be estimated. The LSE of a pixel is estimated by classifying the pixels according to the class that they fall into. When a pixel has an NDVI value that is below −0.185, then the LSE value of 0.995 is assigned to the pixel, when the NDVI value is greater than or equal to −0.185 and less than 0.157, a LSE value of 0.985 is assigned to the pixel, when the NDVI value is greater than or equal to 0.157 and less than or equal to 0.727, a logarithmic relationship between NDVI and LSE is used [19], and finally, when the NDVI is greater than 0.727, the pixel is assigned a value of 0.990 as shown in Table 2.
  • NDVI threshold LSE estimation algorithm
The NDVI threshold emissivity estimation algorithm uses certain NDVI values to distinguish between soil pixels (NDVI < NDVIs) and pixels containing vegetation cover (NDVI > NDVIs). when estimating LSE from NDVI, authors have proposed the use of 0.2 and 0.5 as NDVI of soil and NDVI of vegetation in global scales [22]. According to the NDVI threshold algorithm, three possible classes can be formed for LSE estimation.
  • (NDVI < NDVIs) in this scenario, a pixel is considered to be composed of bare soil or rock, therefore, the LSE of soil is assigned to the pixel.
  • (NDVI > NDVIs) in this scenario, a pixel is considered to be composed of full vegetation cover, therefore, the LSE of vegetation is assigned to the pixel.
  • (NDVIs ≤ NDVI ≤ NDVIv) in this scenario, a pixel is considered to be composed of a mixture of vegetation and rocks/soil; therefore, the authors in [22] introduced Equation (6) to represent the relationship between NDVI and LSE.
ε λ = ε v λ P v + ε s λ ( 1 − P v ) + C λ
where, εv stands for the emissivity of vegetation, εs is the emissivity of soil, Pv is the proportion of vegetation (sometimes known as Fractional Vegetation Cover (FVC)) and C is a value which represents the effect of the geometry of a surface; (C = 0 for flat surfaces). The value of C takes into consideration the cavity effect due to surface roughness. The authors in [22] proposed the cavity term for a mixed area and near-nadir view is given by Equation (7).
C λ = ( 1 − ε s λ ) ε v λ F ′ ( 1 − P v )
The value of F′ is a geometrical factor ranging from 0 to 1, depending on the geometrical distribution of the surface. In the assumption of different geometric distributions of surface roughness, a mean value of 0.55 is chosen [23].
To calculate the proportion of vegetation cover at a pixel level using NDVI, Equation (8) has been used in the plugin [24,25]. The emissivity values shown in Table 3 have been used in the study [18].
P v = [ N D V I − N D V I m i n N D V I m a x − N D V I m i n ] 2
Figure 3 shows the application interface of LSE estimation in the plugin. The option is accessible through the “Land Surface Emissivity tab”.

2.2.4. Estimation of Atmospheric Transmittance, Upwelling and Down-Welling Radiance

In this study, the estimation of the atmospheric parameters (transmittance, upwelling and down-welling radiance) of Landsat 5 TM, Landsat 7 ETM+ and Landsat 8 TIRS has been done through the use of the NASA’S atmospheric parameters calculator [26]. The tool uses the National Centers for Environmental Prediction (NCEP) modeled atmospheric global profiles interpolated at a particular date, time, and location as inputs [26]. Through this tool, the atmospheric profile of the study area has been simulated to the conditions during the satellite overpass time. The tool supports the simulation of atmospheric parameters as from January 2000. The obtained atmospheric transmittance, upwelling and down welling radiance used in the study is as shown in Table 4.

2.2.5. Determination of Near Surface Air Temperature

The near surface air temperatures used in the study have been obtained from the temperature measurements from the meteorological stations. During the determination of the effective mean atmospheric temperature, the mean of the near surface temperatures which have been obtained from the Canadian meteorological stations were used. For Landsat 5 TM, a total of 10 and 9 meteorological stations have been used for the imagery dated 31 December 2010 and 13 November 2010, respectively. For Landsat 7 ETM+, a total of 9 and 8 meteorological stations have been used for the imagery dated 21 November 2010 and 23 February 2016, respectively. For Landsat 8 TIRS, a total of 7 and 9 meteorological stations have been for the imagery dated 4 June 2015 and 29 May 2013, respectively.

2.2.6. Determination of Atmospheric Water Vapor

In order to apply the SCA, the knowledge of the amount of the atmospheric water vapor during the satellite overpass time is crucial. In this study, Equation (9) has been used in the estimation of atmospheric water vapor [27,28,29]. The relative humidity has been obtained from the meteorological data.
w i = 0.0981   × { 10 × 0.6108 × e x p [ 17.27 × ( T 0 − 273.15 ) 237.3 + ( T 0 − 273.15 ) ]   × R H } + 0.1679
In Equation (9), Wi stands for atmospheric water vapor, T0 stands for near surface air temperature in Kelvin and RH stands for relative humidity.

2.2.7. Brightness Temperature Emissivity/Atmospheric Correction

It is important to correct brightness temperature against LSE and atmospheric parameters. There are many algorithms which have been designed for the purpose of estimating LST from thermal infrared satellite imagery. These algorithms vary from one sensor to another and each one of them has its strength, limitations and difference in level of accuracies. In this study, four algorithms have been implemented in the plugin.

Inversion of Planck’s Function

After the land surface emissivity is estimated, the next step is to perform LSE correction. To do so, the Plank function has been programmed in the plugin. Equation (10) shows the Planck function [30,31]. The Planck function corrects the emission of a substance in comparison to a blackbody.
T s = B T { 1 + [ λ · B T ρ ] · l n ε }
where Ts is the land surface temperature (K); BT is the at-sensor brightness temperature (K), λ is the wavelength of the emitted radiance; ρ is the (h * c/σ) = 1.438 · 10−2 mK; and ε is the spectral emissivity. Figure 4 shows the application’s interface for the Plank equation.

Mono-Window Algorithm (MWA)

The brightness temperature obtained from Equation (4) is supposed to be corrected against LSE as well as atmospheric parameters which may affect the data obtained by the sensors in space as a result of absorption and scattering of electromagnetic radiation. To obtain a reasonably high quality LST, Liu and Zhang [29] used the MWA [6]. The algorithm requires the knowledge of LSE, atmospheric transmittance and effective mean atmospheric temperature. Sensitivity analysis on the MWA revealed that the possible error of the ground emissivity which is difficult to estimate has relatively insignificant effect to the probable LST estimation error which is more sensitive to the possible error of transmittance and mean atmospheric temperature [6]. Validation of the simulated data for various situations of seven typical atmospheres points out that the algorithm is able to provide accurate LST retrieval from Landsat 5 TM and Landsat 7 ETM+ data. The difference between the retrieved and simulated ones is less than 0.4 °C for most situations [6]. The MWA is shown by Equation (11) [6].
T s = a i ( 1 − C i − D i ) + [ b i ( 1 − C i − D i ) + C i + D i ] T i − D i T a C i
where TS is the Land Surface Temperature (LST), Ti is the brightness temperature from Equation (4), Ta is the effective mean atmospheric temperature, ai = −67.355351 and bi = 0.458606. The values of Ci and Di can be calculated using the Equations (12) and (13), respectively [6].
C i = ε i τ i
D i = ( 1 − τ i ) [ 1 + ( 1 − ε i ) τ i ]
where εi is the ground surface emissivity and τi is the atmospheric transmittance. Ta, εi and τi are the three parameters needed to convert the brightness temperature to LST. In order to acquire the effective mean atmospheric temperature of an area, Qin et al. [6] introduced linear relations for the approximation depending on the location of the atmosphere of an area of study. These are as shown in Table 5. In this study, the mid-latitude winter was considered to be the most suitable atmospheric profile in accordance to the time and location when the Landsat 5 TM imagery dated 31 December 2010, Landsat 5 TM imagery dated 13 November 2010, Landsat 7 ETM+ imagery dated 21 November 2010 and for Landsat 7 ETM+ imagery dated 23 February 2016. For Landsat 8 TIRS, the mid-latitude summer was considered to be the suitable atmospheric profile for the imagery dated 4 June 2015 and 29 May 2013.
The temperatures Ta and T0 in the equations are both measured in Kelvin. These formulas mean that, in standard atmospheric conditions with a clear sky and absence of turbulence, the effective mean atmospheric temperature Ta is a linear relation of near surface air temperature T0. This is due to the fact that the water vapor distribution and atmospheric temperature distribution on Ta are assumed to be constant for the standard distributions. With the introduction of National Aeronautics and Space Administration (NASA)’s atmospheric parameters calculator [26], it is now possible to calculate the atmospheric transmittance, the up-welling and down-welling radiances. Figure 5 shows the interface of the MWA in the developed application.

Radiative Transfer Equation (RTE)

To obtain LST data from the TIR bands of Landsat 5 TM, Landsat 7 ETM+ and Landsat 8 TIRS sensor, the radiative transfer equation has been implemented in the plugin. The equation uses one thermal infrared channel to extract land surface temperature. The Landsat TM and Landsat ETM+ sensors are fitted out with only one TIR channel. This makes the application of split window algorithms impossible, however, the high spatial resolution of these sensors makes them usable for local and regional studies [32]. Single channel algorithms have been successfully used in several studies for LST retrieval [33,34,35,36]. The RTE is shown in Equation (14).
L λ = [ ε λ L λ ( T s ) + ( 1 − ε λ ) L d o w n ] τ + L λ u p
where Lλ is the thermal infrared radiance received by satellite a sensor which is mainly comprised of surface radiance, down-welling radiance from the atmosphere and up-welling radiance to the atmosphere [37]; τ is the atmospheric transmittance; ε is the land surface emissivity; Lλ(Ts) is the radiance of a blackbody target of kinetic temperature Ts; Lλup is the upwelling or atmospheric path radiance; and Lλdown is the down-welling or sky radiance.
Due to the reason that in most of the conditions, the meteorological data of study areas during the satellite overpass time is usually missing, the NASA’s Atmospheric Correction Parameter Calculator [26] has been used simulate the atmospheric conditions. With the required variables in our hands, the land surface radiance, Lλ(Ts) of kinetic temperature Ts can be calculated by Equation (15).
L λ ( T s ) = L λ − L λ u p τ ε λ − 1 − ε λ ε λ L λ d o w n
With the radiance calculated, the radiances are converted to LST using the relationship which is similar to the Plank equation with two prelaunch calibration constants of K1 and K2. To do this, Equation (16) has been used in the study [31,32].
T s = K 2 l n ( K 1 L λ ( T s ) + 1 )
where Ts is the LST in Kelvin which is changed to degree Celsius (°C) by subtracting 273.15; and K1 and K2 are thermal constants whose values are obtained from the metadata file [32]. Figure 6 shows the application’s interface for the RTE.

Single Channel Algorithm (SCA)

The generalized SCA [38] was developed was developed for the purpose of extracting LST from a single thermal infrared band. According to the algorithm, the LST is expressed as shown in Equations (17)–(19).
T s = γ [ ε − 1 ( ψ 1 L s e n s o r + ψ 2 ) + ψ 3 ] + δ
γ = { C 2 L s e n s o r T s e n s o r 2 [ λ 4 C 1 L s e n s o r + λ − 1 ] } − 1
δ = − γ L s e n s o r + T s e n s o r
From the equations, Ts stands for LST, Tsensor is the at sensor brightness temperature in Kelvin, λ is the effective wavelength of a thermal infrared band in use, C1 = 1.19104 × 108 W·m-2·sr-1·μm4 and C2 = 14387.7 µm·K. The ψ1, ψ2, and ψ3 atmospheric parameters can be estimated from Equations (20)–(22), respectively. Figure 7 shows the application’s interface for the SCA.
ψ 1 = 0.14714 w 2 − 0.15583 w + 1.1234
ψ 2 = − 1.1836 w 2 − 0.3760 w − 0.52894
ψ 3 = − 0.04554 w 2 + 1.8719 w − 0.39071

3. Results

In this study, LST maps for New Brunswick, Canada have been produced using Landsat 5 TM, Landsat 7 ETM+ and Landsat 8 TIRS sensors. The LST has been derived by using the MWA, RTE, Plank function and the SCA. In the study, Zhang, Wang et al.’s LSE [21] algorithm was considered to be the most suitable in the estimation of LSE as it produced the best results in LST extraction. Meteorological data were used to model the accuracy assessment of the LST values derived from the sensors.
Figure 8 shows the LST maps produced by Landsat 5 TM the imagery dated 13 November 2010 from the study area through the use of the SCA, RTE, MWA and the Planck function. The meteorological stations used in the modeling of the accuracy of the derived LST have also been shown on the maps.
Figure 9 shows the LST maps produced by the imagery acquired from Landsat 5 TM on 31 December 2010 from the study area through the use of the SCA, RTE, MWA and the Planck function. The meteorological stations used in the modeling of the accuracy of the derived LST have also been shown on the maps.
Figure 10 shows the LST maps produced by the imagery acquired from Landsat 7 ETM+ on 21 November 2010 from the study area through the use of the SCA, RTE, MWA and the Planck function. The meteorological stations used in the modeling of the accuracy of the derived LST have also been shown on the maps.
Figure 11 shows the LST maps produced by the imagery acquired from Landsat 7 ETM+ on 23 February 2016 from the study area through the use of the SCA, RTE, MWA and the Planck function. The meteorological stations used in the modeling of the accuracy of the derived LST have also been shown on the maps.
Figure 12 shows the LST maps produced by the imagery acquired from Landsat 8 TIRS on 4 June 2015 from the study area through the use of the SCA, RTE, MWA and the Planck function. The meteorological stations used in the modeling of the accuracy of the derived LST have also been shown on the maps.
Figure 13 shows the LST maps produced by the imagery acquired from Landsat 8 TIRS on 29 May 2013 from the study area through the use of the SCA, RTE, MWA and the Planck function. The meteorological stations used in the modeling of the accuracy of the derived LST have also been shown on the maps.
The accuracy analysis of the values obtained from the sensors has been done at pixel level by observing the differences between the near surface temperatures obtained from the meteorological data and the LST stored in a pixel after the use of a land surface temperature extraction algorithm. Table 6 and Table 7 show the values obtained from the algorithms and the near surface temperatures for Landsat 5 TM. Table 8 and Table 9 show the values obtained from the algorithms and the near surface temperatures for Landsat 7 ETM+. Table 10 and Table 11 show the values obtained from the algorithms and the near surface temperatures for Landsat 8 TIRS. All temperatures are in degrees Celsius.
The coordinates of the meteorological stations were converted to shape files and thereafter the land surface temperature existing in the pixel that the meteorological station was located was compared to the value of the near surface temperature which was measured by the meteorological station. In this way, the possible errors were reduced as temperature varies from one place to another due to changes in topography and land cover. The locations of the meteorological stations have been shown on the maps (Figure 8, Figure 9, Figure 10, Figure 11, Figure 12 and Figure 13).

4. Discussion

One of the challenges of deriving land surface temperatures from space is doing the accuracy assessment of the results using values which have been measured on the ground. This is because satellite imagery covers a large area (the approximate scene size is 170 km north to south by 183 km east to west for Landsat 5 TM, Landsat 7 ETM+ and Landsat 8 TIRS) and the area varies with topography, and land cover. Not only that but also due to the fact that the spatial resolution of satellite imagery is usually different from the measurements which are made on the ground.
Atmospheric conditions of a study area can also seriously hamper the accuracy analysis of data, especially the presence of cloud cover. For this reason, the Landsat imagery that has been used in this study had a cloud cover of less than 10%. After LST obtained from Landsat 5 TM, Landsat 7 ETM+ and Landsat 8 TIRS have been resampled, the spatial resolution of the LST values represents 30 meters on the ground. This poses a challenge as a pixel value is an average of an area with 900 m2 on the ground while the data used for the assessment is just a single point (location of a measurement made on the ground) falling within the pixel being assessed. In this study, among the data used included the ones obtained from the Landsat 7 ETM+ satellite sensor; this sensor has been affected by the presence of the Scan Line Corrector (SLC) as from 31 May 2003. Striped imagery was used in the study as a result of the absence of enough meteorological data of the study area before 2003. Some of the meteorological stations used in the accuracy analysis of the Landsat 7 ETM+ image dated 23 February 2016 were affected by stripes. These stations have been shown in Table 9.
The derived LST from the satellite imagery used in this study was compared to meteorological data due to the absence of in situ data during the satellite overpass time. According to studies done by previous researchers [29] there exists a relationship between land surface temperature and near surface temperatures. Land surface temperature can even be used to estimate near surface air temperature [39,40,41]. Authors in [29] used near surface temperatures in the validation of satellite derived land surface temperatures. Due to these reasons, meteorological data was considered as a tool to model the accuracy of the LST derived from the sensors involved in this study.
In this study, the LST values which were obtained from the algorithms produced results which were very close. For Landsat 5 TM, the best results were obtained from the RTE and the Planck function which produced RMSE of 2.64 °C and 1.58 °C, respectively. For Landsat 7 ETM+ LST retrieval, the best results were obtained from the RTE and the Planck function with RMSE of 3.75 °C and 3.58 °C respectively. For Landsat 8 TIRS LST retrieval, the best results were obtained from the Planck function and the SCA with RMSE of 2.07 °C and 3.06 °C respectively. This implies that all the algorithms are suitable for the retrieval of LST from the sensors used in the study. The differences in results may be due to the errors which may arise from the simulation of atmospheric parameters, estimation of LSE as well as errors which may arise from the estimation of atmospheric water vapor. It should also be noted that the temperature values used in the accuracy assessment were near surface temperatures measured from the meteorological stations and not actual temperatures of the ground. These temperatures are measured at a height of up to two meters above the ground.
Generally, the LST derived from the algorithms had differences of less than 1 °C between each other. The SCA produced LST which were generally lower in comparison to the other algorithms. The SCA has been found to be less accurate in LST estimation by other authors as a result of being highly dependent in atmospheric [24,42] parameters. The algorithm is also affected by errors that may arise from the estimation of land surface emissivities, which can contribute to an error between 0.2 K and 0.7 K [42].
The atmospheric water vapor used in the SCA was obtained from meteorological data. As a result of the data being obtained from points and thereafter an average value being used for a wide area, this may result to bigger discrepancies between the near surface temperatures and the LST values from space. For future studies, it is important to use atmospheric water vapor values which have been obtained from satellites for LST estimation. Atmospheric water vapor is the most variable component of the atmosphere and varies both horizontally and vertically and it may affect the NDVI of an area [20]. The atmospheric effects arising from clouds, water vapor and aerosols may reduce the contrast between the visible and near infrared reflections and as a result lead to smaller NDVI values [20] hence affecting the LSE which is derived from NDVI.

5. Conclusions

The LST maps in this study have been produced using a QIS plugin which was written using the Python programming language. This has demonstrated how powerful open source tools can be in terms of remote sensing and geographic information systems data processing and analysis. As a result of this study being done through the use of open source tools, it encourages the use of Landsat 5 TM, Landsat 7 ETM+ and Landsat 8 TIRS data in the production of land surface temperature maps, the improvement of algorithms and also provides a platform for debates concerning LST extraction from satellite imagery. This plugin enables users from other disciplines such a hydrology, landscape engineering, urban planning and other related fields to easily extract land surface temperature from visible and near infrared satellite imagery. This plugin also enables users to make use of more than one algorithm of land surface temperature estimation and as a result managing to benefit from possible advantages of them while reducing the possible disadvantages which may arise from the use of a single algorithm in land surface temperature estimation.
This study also revealed that the MWA, the SCA, the RTE and the Planck function can all be used in the estimation of LST from Landsat 5 TM, Landsat 7 ETM+ and Landsat 8 TIRS. The availability of meteorological data during the satellite overpass time can, however, play a great role towards determining the algorithm to be used in the estimation of LST. This therefore makes the Planck function easier to use in comparison to the other algorithms as it does not require atmospheric variables during the satellite overpass time.
In order to obtain better results in the future, researchers should do accuracy analysis of the LST derived from satellite imagery with actual in situ data collected during the satellite overpass time. It is also important to estimate the atmospheric water vapor from satellites instead of using meteorological data as satellite derived atmospheric water vapor provide better average values over spatially large areas than meteorological data. To perform atmospheric correction of NDVI and emissivity, atmospheric correction algorithms such as 6S [43] and Moderate Resolution Atmospheric Transmission (MODTRAN) [44] may be applied.

Acknowledgments

This study was supported by Anadolu University Scientific Research Projects Commission under the grant number 1601F031. We sincerely appreciate the research funding, it was provided by Anadolu University for this article. We express high gratitude to the United States Geological Survey (USGS) for making the Landsat data freely available. We would also like to give our special thanks go to the Government of Canada (Environment Canada) for making the historic meteorological data available for use. This has made this study successful in a great extent. Last but not least, we would like to express our deepest gratitude to the QGIS development team who made it possible for us to make use of the QGIS application programming interface for data processing though the development of the plugin.

Author Contributions

Both authors contributed in the design and methodology of the study. The Python Plugin was written by Milton Isaya Ndossi while Ugur Avdan worked on the methodology, accuracy assessment, interpretation of the results and the final version of this contribution.

Conflicts of Interest

The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

Abbreviations

The following abbreviations are used in this manuscript:
NSTNear Surface Temperature
RHRelative Humidity
VNIRVisible and Near Infrared
TIRThermal Infrared
LSTLand Surface Temperature
LSELand Surface Emissivity
SCASingle Channel Algorithm
MWAMono Window Algorithm
RTERadiative Transfer Equation
RSRemote Sensing
GISGeographic Information Systems
NDVINormalized Difference Vegetation Index
TOATop of Atmosphere
EOSEarth Observation Satellites
USGSUnited States Geological Survey
UTCCoordinated Universal Time
GDALGeospatial Data Abstraction Library
DNDigital Numbers
KKelvin
OLIOperational Land Imager
TMThematic Mapper
ETM+Enhanced Thematic Mapper Plus
NASANational Aeronautics and Space Administration
NCEPNational Centers for Environmental Prediction

References

  1. Wu, P.; Shen, H.; Zhang, L.; Göttsche, F.-M. Integrated fusion of multi-scale polar-orbiting and geostationary satellite observations for the mapping of high spatial and temporal resolution land surface temperature. Remote Sens. Environ. 2015, 156, 169–181. [Google Scholar]
  2. Li, Z.-L.; Tang, B.-H.; Wu, H.; Ren, H.; Yan, G.; Wan, Z.; Trigo, I.F.; Sobrino, J.A. Satellite-derived land surface temperature: Current status and perspectives. Remote Sens. Environ. 2013, 131, 14–37. [Google Scholar] [CrossRef]
  3. Pandya, M.R.; Shah, D.B.; Trivedi, H.J.; Darji, N.P.; Ramakrishnan, R.; Panigrahy, S.; Parihar, J.S.; Kirankumar, A.S. Retrieval of land surface temperature from the Kalpana-1 VHRR data using a single-channel algorithm and its validation over western India. ISPRS J. Photogramm. Remote Sens. 2014, 94, 160–168. [Google Scholar] [CrossRef]
  4. Jiménez-Muñoz, J.-C.; Sobrino, J.A. Split-window coefficients for land surface temperature retrieval from low-resolution thermal infrared sensors. IEEE Geosci. Remote Sens. Lett. 2008, 5, 806–809. [Google Scholar] [CrossRef]
  5. Jiménez-Muñoz, J.C.; Sobrino, J.A. A single-channel algorithm for land-surface temperature retrieval from aster data. IEEE Geosci. Remote Sens. Lett. 2010, 7, 176–179. [Google Scholar] [CrossRef]
  6. Qin, Z.-H.; Karnieli, A.; Berliner, P. A mono-window algorithm for retrieving land surface temperature from landsat TM data and its application to the Israel-Egypt border region. Int. J. Remote Sens. 2001, 22, 3719–3746. [Google Scholar] [CrossRef]
  7. Team, Q.D. QGIS Python Plugins Repository. Available online: http://plugins.qgis.org/ (accessed on 11 May 2016).
  8. USGS. Earth Explorer. Available online: http://earthexplorer.usgs.gov/ (accessed on 5 May 2016).
  9. Historical Data. Available online: http://climate.weather.gc.ca/ (accessed on 5 May 2016).
  10. USGS. Using the USGS Landsat 8 Product. Available online: http://landsat.usgs.gov/Landsat8_Using_Product.php (accessed on 9 November 2014).
  11. USGS. Landsat 8 (L8) Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS). Available online: http://landsat.usgs.gov/calibration_notices.php (accessed on 9 February 2016).
  12. Barsi, J.A.; Schott, J.R.; Hook, S.J.; Raqueno, N.G.; Markham, B.L.; Radocinski, R.G. Landsat-8 thermal infrared sensor (TIRS) vicarious radiometric calibration. Remote Sens. 2014, 6, 11607–11626. [Google Scholar] [CrossRef]
  13. USGS. Frequently Asked Questions about the Landsat Missions. Available online: http://landsat.usgs.gov/how_is_radiance_calculated.php (accessed on 24 March 2016).
  14. Finn, M.P.; Reed, M.D.; Yamamoto, K.H. A Straight forward Guide for Processing Radiance and Reflectance for eo-1 ALI, Landsat 5 tm, Landsat 7 ETM+, and Aster; Unpublished Report; USGS/Center of Excellence for Geospatial Information Science: Washington, DC, USA, 2012. [Google Scholar]
  15. Kruse, P.W.; McGlauchlin, L.D.; McQuistan, R.B. Elements of Infrared Technology: Generation, Transmission and Detection; Wiley: New York, NY, USA, 1962; Volume 1962. [Google Scholar]
  16. Wang, S.L.L. Chapter 8—Land-surface temperature and thermal infrared emissivity. In Advanced Remote Sensing; Wang, S.L.L., Ed.; Academic Press: Boston, FL, USA, 2012; pp. 235–271. [Google Scholar]
  17. Tang, B.-H.; Wu, H.; Li, C.; Li, Z.-L. Estimation of broadband surface emissivity from narrowband emissivities. Opt. Express 2011, 19, 185–192. [Google Scholar] [CrossRef] [PubMed]
  18. Wang, F.; Qin, Z.; Song, C.; Tu, L.; Karnieli, A.; Zhao, S. An improved mono-window algorithm for land surface temperature retrieval from landsat 8 thermal infrared sensor data. Remote Sens. 2015, 7, 4268–4289. [Google Scholar] [CrossRef]
  19. Vandegriend, A.; Owe, M.; Vugts, H.; Ramothwa, G. Botswana Water and Surface Energy Balance Research Program. Part 1: Integrated Approach and Field Campaign Results; NASA Goddard Space Flight Center: Greenbelt, MD, USA, 1992. [Google Scholar]
  20. Srivastava, P.K.; Han, D.; Rico-Ramirez, M.A.; Bray, M.; Islam, T.; Gupta, M.; Dai, Q. Estimation of land surface temperature from atmospherically corrected landsat TM image using 6S and NCEP global reanalysis product. Environ. Earth Sci. 2014, 72, 5183–5196. [Google Scholar] [CrossRef]
  21. Zhang, J.; Wang, Y.; Li, Y. A C++ program for retrieving land surface temperature from the data of landsat TM/ETM+ band6. Comput. Geosci. 2006, 32, 1796–1805. [Google Scholar] [CrossRef]
  22. Sobrino, J.; Raissouni, N. Toward remote sensing methods for land cover dynamic monitoring: Application to morocco. Int. J. Remote Sens. 2000, 21, 353–366. [Google Scholar] [CrossRef]
  23. Sobrino, J.; Caselles, V.; Becker, F. Significance of the remotely sensed thermal infrared measurements obtained over a citrus orchard. ISPRS J. Photogramm. Remote Sens. 1990, 44, 343–354. [Google Scholar] [CrossRef]
  24. Yu, X.; Guo, X.; Wu, Z. Land surface temperature retrieval from landsat 8 TIRS—Comparison between radiative transfer equation-based method, split window algorithm and single channel method. Remote Sens. 2014, 6, 9829–9852. [Google Scholar] [CrossRef]
  25. Carlson, T.N.; Ripley, D.A. On the relation between NDVI, fractional vegetation cover, and leaf area index. Remote Sens. Environ. 1997, 62, 241–252. [Google Scholar] [CrossRef]
  26. Barsi, J.; Barker, J.L.; Schott, J.R. An atmospheric correction parameter calculator for a single thermal band earth-sensing instrument. In Proceedings of the 2003 IEEE International Geoscience and Remote Sensing Symposium, Toulouse, France, 21–25 July 2003.
  27. Yang, J.; Qiu, J. The empirical expressions of the relation between precipitable water and ground water vapor pressure for some areas in china. Sci. Atmos. Sin. 1996, 20, 620–626. [Google Scholar]
  28. Li, J. Estimating land surface temperature from landsat-5 TM. Remote Sens. Technol. Appl. 2006, 21, 322–326. [Google Scholar]
  29. Liu, L.; Zhang, Y. Urban heat island analysis using the landsat TM data and ASTER data: A case study in Hong Kong. Remote Sens. 2011, 3, 1535–1552. [Google Scholar] [CrossRef]
  30. Artis, D.A.; Carnahan, W.H. Survey of emissivity variability in thermography of urban areas. Remote Sens. Environ. 1982, 12, 313–329. [Google Scholar] [CrossRef]
  31. Sinha, S.; Pandey, P.C.; Sharma, L.K.; Nathawat, M.S.; Kumar, P.; Kanga, S. Remote estimation of land surface temperature for different lulc features of a moist deciduous tropical forest region. In Remote Sensing Applications in Environmental Research; Springer: Berlin, Germany; Heidelberg, Germany, 2014; pp. 57–68. [Google Scholar]
  32. Srivastava, P.; Majumdar, T.; Bhattacharya, A.K. Surface temperature estimation in singhbhum shear zone of India using landsat-7 ETM+ thermal infrared data. Adv. Space Res. 2009, 43, 1563–1574. [Google Scholar] [CrossRef]
  33. Goetz, S.; Halthore, R.; Hall, F.; Markham, B. Surface temperature retrieval in a temperate grassland with multiresolution sensors. J. Geophys. Res. 1995, 100, 25397–25410. [Google Scholar] [CrossRef]
  34. Wukelic, G.; Gibbons, D.; Martucci, L.; Foote, H. Radiometric calibration of landsat thematic mapper thermal band. Remote Sens. Environ. 1989, 28, 339–347. [Google Scholar] [CrossRef]
  35. Schott, J.R.; Volchok, W.J. Thematic mapper thermal infrared calibration. Photogramm. Eng. Remote Sens. 1985, 51, 1351–1357. [Google Scholar]
  36. Weng, Q.; Lu, D.; Schubring, J. Estimation of land surface temperature–vegetation abundance relationship for urban heat island studies. Remote Sens. Environ. 2004, 89, 467–483. [Google Scholar] [CrossRef]
  37. Barsi, J.A.; Schott, J.; Palluconi, F.D.; Helder, D.; Hook, S.; Markham, B.; Chander, G.; O’Donnell, E. Landsat tm and ETM+ thermal band calibration. Can. J. Remote Sens. 2003, 29, 141–153. [Google Scholar] [CrossRef]
  38. Jiménez-Muñoz, J.C.; Sobrino, J.A. A generalized single-channel method for retrieving land surface temperature from remote sensing data. J. Geophys. Res. 2003, 108. [Google Scholar] [CrossRef]
  39. Cresswell, M.; Morse, A.; Thomson, M.; Connor, S. Estimating surface air temperatures, from meteosat land surface temperatures, using an empirical solar zenith angle model. Int. J. Remote Sens. 1999, 20, 1125–1132. [Google Scholar] [CrossRef]
  40. Prihodko, L.; Goward, S.N. Estimation of air temperature from remotely sensed surface observations. Remote Sens. Environ. 1997, 60, 335–346. [Google Scholar] [CrossRef]
  41. Stisen, S.; Sandholt, I.; Nørgaard, A.; Fensholt, R.; Eklundh, L. Estimation of diurnal air temperature using MSG SEVIRI data in West Africa. Remote Sens. Environ. 2007, 110, 262–274. [Google Scholar] [CrossRef]
  42. Jiménez-Muñoz, J.; Sobrino, J. Error sources on the land surface temperature retrieved from thermal infrared single channel remote sensing data. Int. J. Remote Sens. 2006, 27, 999–1014. [Google Scholar] [CrossRef]
  43. Kotchenova, S.Y.; Vermote, E.F. Validation of a vector version of the 6s radiative transfer code for atmospheric correction of satellite data. Part II. Homogeneous Lambertian and anisotropic surfaces. Appl. Opt. 2007, 46, 4455–4464. [Google Scholar] [CrossRef] [PubMed]
  44. Berk, A.; Bernstein, L.S.; Robertson, D.C. Modtran: A Moderate Resolution Model for Lowtran; DTIC Document: Belvoir, VA, USA, 1987. [Google Scholar]
Remotesensing 08 00413 g001 1024
Figure 1. Application’s Interface for calculating radiance.
Figure 1. Application’s Interface for calculating radiance.
Remotesensing 08 00413 g001
Remotesensing 08 00413 g002 1024
Figure 2. Application’s Interface for calculating brightness temperature.
Figure 2. Application’s Interface for calculating brightness temperature.
Remotesensing 08 00413 g002
Remotesensing 08 00413 g003 1024
Figure 3. Application’s interface for Land Surface Emissivity (LSE) estimation based on the NDVI image.
Figure 3. Application’s interface for Land Surface Emissivity (LSE) estimation based on the NDVI image.
Remotesensing 08 00413 g003
Remotesensing 08 00413 g004 1024
Figure 4. Application’s interface for the Planck Equation.
Figure 4. Application’s interface for the Planck Equation.
Remotesensing 08 00413 g004
Remotesensing 08 00413 g005 1024
Figure 5. Application’s interface for Mono-Window algorithm.
Figure 5. Application’s interface for Mono-Window algorithm.
Remotesensing 08 00413 g005
Remotesensing 08 00413 g006 1024
Figure 6. Application’s interface for the RTE.
Figure 6. Application’s interface for the RTE.
Remotesensing 08 00413 g006
Remotesensing 08 00413 g007 1024
Figure 7. Application’s interface for the SCA.
Figure 7. Application’s interface for the SCA.
Remotesensing 08 00413 g007
Remotesensing 08 00413 g008a 1024Remotesensing 08 00413 g008b 1024Remotesensing 08 00413 g008c 1024
Figure 8. LST maps derived from Landsat 5 TM: (a) LST extracted using the SCA; (b) LST extracted using the RTE; (c) LST extracted using the MWA; and (d) LST extracted using the Planck equation.
Figure 8. LST maps derived from Landsat 5 TM: (a) LST extracted using the SCA; (b) LST extracted using the RTE; (c) LST extracted using the MWA; and (d) LST extracted using the Planck equation.
Remotesensing 08 00413 g008aRemotesensing 08 00413 g008bRemotesensing 08 00413 g008c
Remotesensing 08 00413 g009a 1024Remotesensing 08 00413 g009b 1024
Figure 9. LST maps derived from Landsat 7 ETM+: (a) LST extracted using the SCA; (b) LST extracted using the RTE; (c) LST extracted using the MWA; and (d) LST extracted using the Planck equation.
Figure 9. LST maps derived from Landsat 7 ETM+: (a) LST extracted using the SCA; (b) LST extracted using the RTE; (c) LST extracted using the MWA; and (d) LST extracted using the Planck equation.
Remotesensing 08 00413 g009aRemotesensing 08 00413 g009b
Remotesensing 08 00413 g010a 1024Remotesensing 08 00413 g010b 1024
Figure 10. LST maps derived from Landsat 7 ETM+: (a) LST extracted using the SCA; (b) LST extracted using the RTE; (c) LST extracted using the MWA; and (d) LST extracted using the Planck equation.
Figure 10. LST maps derived from Landsat 7 ETM+: (a) LST extracted using the SCA; (b) LST extracted using the RTE; (c) LST extracted using the MWA; and (d) LST extracted using the Planck equation.
Remotesensing 08 00413 g010aRemotesensing 08 00413 g010b
Remotesensing 08 00413 g011a 1024Remotesensing 08 00413 g011b 1024
Figure 11. LST maps derived from Landsat 7 ETM+: (a) LST extracted using the SCA; (b) LST extracted using the RTE; (c) LST extracted using the MWA; and (d) LST extracted using the Planck equation.
Figure 11. LST maps derived from Landsat 7 ETM+: (a) LST extracted using the SCA; (b) LST extracted using the RTE; (c) LST extracted using the MWA; and (d) LST extracted using the Planck equation.
Remotesensing 08 00413 g011aRemotesensing 08 00413 g011b
Remotesensing 08 00413 g012a 1024Remotesensing 08 00413 g012b 1024
Figure 12. LST maps derived from Landsat 8 TIRS: (a) LST extracted using the SCA; (b) LST extracted using the RTE; (c) LST extracted using the MWA; and (d) LST extracted using the Planck equation.
Figure 12. LST maps derived from Landsat 8 TIRS: (a) LST extracted using the SCA; (b) LST extracted using the RTE; (c) LST extracted using the MWA; and (d) LST extracted using the Planck equation.
Remotesensing 08 00413 g012aRemotesensing 08 00413 g012b
Remotesensing 08 00413 g013a 1024Remotesensing 08 00413 g013b 1024
Figure 13. LST maps derived from Landsat 8 TIRS: (a) LST extracted using the SCA; (b) LST extracted using the RTE; (c) LST extracted using the MWA; and (d) LST extracted using the Planck equation.
Figure 13. LST maps derived from Landsat 8 TIRS: (a) LST extracted using the SCA; (b) LST extracted using the RTE; (c) LST extracted using the MWA; and (d) LST extracted using the Planck equation.
Remotesensing 08 00413 g013aRemotesensing 08 00413 g013b
Table
Table 1. Landsat scenes used in this study.
Table 1. Landsat scenes used in this study.
SensorDateTime (UTC)PathRow
Landsat 5 TM13 November 201014:57009028
Landsat 5 TM31 December 201014:57009028
Landsat 7 ETM+21 November 201015:00009028
Landsat 7 ETM+23 February 201615:09009028
Landsat 8 TIRS/OLI29 May 201315:09009028
Landsat 8 TIRS/OLI4 June 201515:06009028
Table
Table 2. Algorithm based on the NDVI image.
Table 2. Algorithm based on the NDVI image.
NDVILSE
NDVI < −0.1850.995
−0.185 ≤ NDVI < 0.1570.985
0.157 ≤ NDVI ≤ 0.7271.009 + 0.047 × ln(NDVI)
NDVI > 0.7270.990
Table
Table 3. LSE values for various terrestrial materials [18].
Table 3. LSE values for various terrestrial materials [18].
Terrestrial MaterialSoilBuilt-Up AreasVegetationWater
Emissivity0.9660.9620.9730.991
Table
Table 4. Simulated atmospheric parameters for the study.
Table 4. Simulated atmospheric parameters for the study.
Sensor TypeDate Time (UTC)Atmospheric TransmittanceUpwelling Radiance (W/(m2·sr·µm)) Down-Welling Radiance (W/(m2·sr·µm))
Landsat 5 TM13 November 201014:570.890.721.20
Landsat 5 TM31 December 201014:570.890.580.97
Landsat 7 ETM+21 November 201015:000.960.180.31
Landsat 7 ETM+23 February 201615:090.970.110.20
Landsat 8 TIRS4 June 201515:060.870.911.52
Landsat 8 TIRS29 May 201315:090.890.721.22
Table
Table 5. Linear relations for the approximation of effective mean atmospheric temperature Ta from the near surface air temperature [6].
Table 5. Linear relations for the approximation of effective mean atmospheric temperature Ta from the near surface air temperature [6].
AtmospheresLinear Relations Equations
For USA 1976Ta = 25.9396 + 0.88045T0
Tropical modelTa = 17.9769 + 0.91715T0
Mid-latitude SummerTa = 16.0110 + 0.92621T0
Mid-latitude winterTa = 19.2704 + 0.91118T0
Table
Table 6. Comparison between LST values obtained from Landsat 5 TM using the MWA, Planck equation, RTE, SCA and meteorological data obtained on 31 December 2010.
Table 6. Comparison between LST values obtained from Landsat 5 TM using the MWA, Planck equation, RTE, SCA and meteorological data obtained on 31 December 2010.
Station NameLatitudeLongitudeElevation (Meters)NST (°C) at 12:00 RH (%)MWA (°C)PLANCK (°C)RTE (°C)SCA (°C)
Mechanic Settlement45°41′37.040″N65°09′54.040″W403−3.199−4.75−4.70−4.55−5.56
Fundy Park (Alma) Cs45°36′00.000″N64°57′00.000″W42.73.662−0.91−1.29−0.73−2.17
Buctouche Cda Cs46°25′49.006″N64°46′05.009″W35.9−0.978−8.06−7.64−7.85−8.48
Nappan Auto45°45′34.400″N64°14′29.200″W19.8−7.973−3.46−3.55−3.26−4.42
Saint John A45°19′05.000″N65°53′08.050″W108.8−0.269−0.29−0.73−0.11−1.62
Parrsboro45°24′48.000″N64°20′49.000″W30.92.168−0.29−0.73−0.11−1.62
Gagetown A45°50′00.000″N66°26′00.000″W50.6−2.279−0.91−1.29−0.73−2.17
Gagetown Awos A45°50′20.000″N66°26′59.000″W50.6−177−0.29−0.73−0.11−1.62
Summerside46°26′28.000″N63°50′17.000″W12.2−0.978−3.46−3.55−3.26−4.42
Mechanic Settlement45°41′37.040″N65°09′54.040″W403−3.199−4.75−4.70−4.55−5.56
Root Mean Square Error (RMSE) 2.672.782.643.17
NST Near Surface Temperature (Day light saving was observed); PLANCK Planck Equation; RTE Radiative Transfer Equation; SCA Single Channel Algorithm; RH Relative Humidity.
Table
Table 7. Comparison between LST values obtained from Landsat 5 TM using the MWA, Planck equation, RTE, SCA and meteorological data obtained on 13 November 2010.
Table 7. Comparison between LST values obtained from Landsat 5 TM using the MWA, Planck equation, RTE, SCA and meteorological data obtained on 13 November 2010.
Station NameLatitudeLongitudeElevation (Meters)NST (°C) at 12:00RH (%)MWA (°C)PLANCK (°C)RTE (°C)SCA (°C)
Mechanic Settlement45°41′37.040″N65°09′54.040″W40315.44015.1514.9715.3213.28
Fundy Park (Alma) Cs45°36′00.000″N64°57′00.000″W42.713.15315.7015.8015.7310.86
Buctouche Cda Cs46°25′49.006″N64°46′05.009″W35.912.43912.7012.5812.9010.36
Nappan Auto45°45′34.400″N64°14′29.200″W19.812.15214.1714.2614.3010.36
Saint John A45°19′05.000″N65°53′08.050″W108.812.24114.8314.2515.1413.27
Parrsboro45°24′48.000″N64°20′49.000″W30.913.24814.5614.5714.7010.85
Gagetown A45°50′00.000″N66°26′00.000″W50.611.63610.1110.2710.408.37
Gagetown Awos A45°50′20.000″N66°26′59.000″W50.612.03610.9610.8211.309.86
Summerside46°26′28.000″N63°50′17.000″W12.213.64012.6312.5512.8810.36
Root Mean Square Error (RMSE) 1.641.581.682.33
NST Near Surface Temperature (Day light saving was observed); PLANCK Planck Equation; RTE Radiative Transfer Equation; SCA Single Channel Algorithm; RH Relative Humidity.
Table
Table 8. Comparison between LST values obtained from Landsat 7 ETM+ using the MWA, Planck equation, RTE, SCA and meteorological data obtained on 21 November 2010.
Table 8. Comparison between LST values obtained from Landsat 7 ETM+ using the MWA, Planck equation, RTE, SCA and meteorological data obtained on 21 November 2010.
Station NameLatitudeLongitudeElevation (Meters)NST (°C) at 12:00 RH (%)MWA (°C)PLANCK (°C)RTE (°C)SCA (°C)
Mechanic Settlement45°41′37.040″N65°09′54.040″W403−6.878−11.75−11.52−11.43−12.34
Fundy Park (Alma) Cs45°36′00.000″N64°57′00.000″W42.7−1.958−5.50−5.53−5.21−6.39
Buctouche Cda Cs46°25′49.006″N64°46′05.009″W35.9−4.464−9.62−9.48−9.31−10.31
Nappan Auto45°45′34.400″N64°14′29.200″W19.8−4.470−8.23−8.14−7.92−8.98
Saint John A45°19′05.000″N65°53′08.050″W108.8−3.655−6.18−6.18−5.88−7.03
Parrsboro45°24′48.000″N64°20′49.000″W30.9−2.660N/A N/A N/A N/A
Gagetown Awos A45°50′20.000″N66°26′59.000″W50.6−2.955−3.51−3.62−3.22−4.49
Summerside46°26′28.000″N63°50′17.000″W12.2−368−8.23−8.14−7.92−8.98
Mechanic Settlement45°41′37.040″N65°09′54.040″W403−6.878−11.75−11.52−11.43−12.34
Root Mean Square Error (RMSE) 4.033.943.754.73
NST Near Surface Temperature (Day light saving was observed); PLANCK Planck Equation; RTE Radiative Transfer Equation; SCA Single Channel Algorithm; RH Relative Humidity; N/A—Data not available as a result of Scan Line Corrector (SLC).
Table
Table 9. Comparison between LST values obtained from Landsat 7 ETM+ using the MWA, Planck equation, RTE, SCA and meteorological data obtained on 23 February 2016.
Table 9. Comparison between LST values obtained from Landsat 7 ETM+ using the MWA, Planck equation, RTE, SCA and meteorological data obtained on 23 February 2016.
Station NameLatitudeLongitudeElevation (Meters)NST (°C) at 12:00 RH (%)MWA (°C)PLANCK (°C)RTE (°C)SCA (°C)
Mechanic Settlement45°41′37.040″N65°09′54.040″W403−8.135−14.46−14.31−14.04−15.12
Fundy Park (Alma) Cs45°36′00.000″N64°57′00.000″W42.7−6.950−10.17−10.15−9.77−10.98
Buctouche Cda Cs46°25′49.006″N64°46′05.009″W35.9−8.945−5.40−5.53−5.03−6.39
Nappan Auto45°45′34.400″N64°14′29.200″W19.8−7.968N/AN/AN/AN/A
Parrsboro45°24′48.000″N64°20′49.000″W30.9−5.949−7.41−7.48−7.03−8.33
Moncton lnternational Airport46°06′44.000″N64°40′43.000″W70.7−844−4.08−4.25−3.72−5.12
Gagetown A45°50′00.000″N66°26′00.000″W50.6−8.447N/AN/AN/AN/A
Summerside46°26′28.000″N63°50′17.000″W12.2−7.355N/AN/AN/AN/A
Root Mean Square Error (RMSE) 3.683.583.613.80
NST Near Surface Temperature (Day light saving was observed); PLANCK Planck Equation; RTE Radiative Transfer Equation; SCA Single Channel Algorithm; RH Relative Humidity; N/A—Data not available as a result of Scan Line Corrector (SLC).
Table
Table 10. Comparison between LST values obtained from Landsat 8 TIRS using the MWA, Planck equation, RTE, SCA and meteorological data obtained on 4 June 2015.
Table 10. Comparison between LST values obtained from Landsat 8 TIRS using the MWA, Planck equation, RTE, SCA and meteorological data obtained on 4 June 2015.
Station NameLatitudeLongitudeElevation (m)NST (°C) at 11:00 amRH (%)MWA (°C)PLANCK (°C)RTE (°C)SCA (°C)
Mechanic Settlement45°41′37.040″N65°09′54.040″W40313.1579.279.659.438.15
Fundy Park (Alma) Cs45°36′00.000″N64°57′00.000″W42.712.45613.4213.3613.5111.33
Buctouche Cda Cs46°25′49.006″N64°46′05.009″W35.915.15816.6616.1416.7314.20
Nappan Auto45°45′34.400″N64°14′29.200″W19.815.46314.8414.4015.0113.18
Parrsboro45°24′48.000″N64°20′49.000″W30.915.94211.7012.1011.699.17
Summerside46°26′28.000″N63°50′17.000″W12.213.96213.4413.5513.4510.82
Moncton International Airport46°06′44.000″N64°40′43.000″W70.716.94417.7417.5517.8414.25
Root Mean Square Error (RMSE) 2.302.072.283.65
NST Near Surface Temperature; PLANCK Planck Equation; RTE Radiative Transfer Equation; SCA Single Channel Algorithm; RH Relative Humidity.
Table
Table 11. Comparison between LST values obtained from Landsat 8 TIRS using the MWA, Planck equation, RTE, SCA and meteorological data obtained on 29 May 2013.
Table 11. Comparison between LST values obtained from Landsat 8 TIRS using the MWA, Planck equation, RTE, SCA and meteorological data obtained on 29 May 2013.
Station NameLatitudeLongitudeElevation (Meters)NST (°C) at 11:00 amRH (%)MWA (°C)PLANCK (°C)RTE (°C)SCA (°C)
Mechanic Settlement45°41′37.040″N65°09′54.040″W403172722.2421.8023.0718.85
Fundy Park (Alma) Cs45°36′00.000″N64°57′00.000″W42.718.13818.4118.4119.3115.53
Buctouche Cda Cs46°25′49.006″N64°46′05.009″W35.920.13523.0722.2823.9820.62
Nappan Auto45°45′34.400″N64°14′29.200″W19.819.33521.6621.0622.5919.29
Parrsboro45°24′48.000″N64°20′49.000″W30.916.24922.0421.6022.8818.77
Moncton International Airport46°06′44.000″N64°40′43.000″W70.720.24129.7428.7830.2423.70
Gagetown A45°50′00.000″N66°26′00.000″W50.6203115.3617.5816.3413.32
Gagetown Awos A45°50′20.000″N66°26′59.000″W50.6203122.0221.7422.7918.06
Summerside46°26′28.000″N63°50′17.000″W12.216.15020.2619.6721.2618.65
Root Mean Square Error (RMSE) 4.834.335.353.06
NST Near Surface Temperature; PLANCK Planck Equation; RTE Radiative Transfer Equation; SCA Single Channel Algorithm; RH Relative Humidity.

Share and Cite

MDPI and ACS Style

Isaya Ndossi, M.; Avdan, U. Application of Open Source Coding Technologies in the Production of Land Surface Temperature (LST) Maps from Landsat: A PyQGIS Plugin. Remote Sens. 2016, 8, 413. https://doi.org/10.3390/rs8050413

AMA Style

Isaya Ndossi M, Avdan U. Application of Open Source Coding Technologies in the Production of Land Surface Temperature (LST) Maps from Landsat: A PyQGIS Plugin. Remote Sensing. 2016; 8(5):413. https://doi.org/10.3390/rs8050413

Chicago/Turabian Style

Isaya Ndossi, Milton, and Ugur Avdan. 2016. "Application of Open Source Coding Technologies in the Production of Land Surface Temperature (LST) Maps from Landsat: A PyQGIS Plugin" Remote Sensing 8, no. 5: 413. https://doi.org/10.3390/rs8050413

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop