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

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.

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. Landsat scenes used in this study.

Sensor	Date	Time (UTC)	Path	Row
Landsat 5 TM	13 November 2010	14:57	009	028
Landsat 5 TM	31 December 2010	14:57	009	028
Landsat 7 ETM+	21 November 2010	15:00	009	028
Landsat 7 ETM+	23 February 2016	15:09	009	028
Landsat 8 TIRS/OLI	29 May 2013	15:09	009	028
Landsat 8 TIRS/OLI	4 June 2015	15:06	009	028

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_λ = M_LQ_cal + A_L − O_i

(1)

where L_λ is the Top-of-Atmosphere (TOA) spectral radiance in W/(m²·sr·μm), M_L is the band specific multiplicative rescaling factor from the metadata, A_L is the band-specific additive rescaling factor from the metadata, Q_cal is the quantized and calibrated standard product pixel values (DN), and O_i is the offsets issued by USGS for the calibration of the TIRS bands. The O_i 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

(2)

L_λ = ((LMAX_λ − LMIN_λ)/(QCALMAX − QCALMIN)) * (QCAL − QCALMIN) + LMIN_λ

(3)

where L_λ is the spectral radiance at the sensor’s aperture in W/(m²·sr·μm), gain is the rescaled gain (the data product “gain” contained in the Level 1 product header or ancillary data record) in W/(m²·sr·μm), bias is the rescaled bias in W/(m²·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/(m²·sr·μm), LMAX_λ is the spectral radiance that is scaled to QCALMAX in W/(m²·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_{sen}= \frac{K_2}{ln\left(\frac{K_1}{L_{\mathsf{\lambda }}}+1\right)}$$

(4)

In Equation (4), T_sen stands for the at-sensor brightness temperature (K), L_λ is the TOA spectral radiance, K₁ is the band specific thermal conversion constant from the metadata and K₂ band specific thermal conversion constant from the metadata. The same equation is used for all sensors. The K₁ and K₂ values may vary depending on the sensor and the wavelengths by which the thermal bands operate. The values of K₁ and K₂ 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.

$$NDVI=\frac{NIR-R}{NIR+R}$$

(5)

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. Algorithm based on the NDVI image.

NDVI	LSE
NDVI < −0.185	0.995
−0.185 ≤ NDVI < 0.157	0.985
0.157 ≤ NDVI ≤ 0.727	1.009 + 0.047 × ln(NDVI)
NDVI > 0.727	0.990

.

• NDVI threshold LSE estimation algorithm

The NDVI threshold emissivity estimation algorithm uses certain NDVI values to distinguish between soil pixels (NDVI < NDVI_s) and pixels containing vegetation cover (NDVI > NDVI_s). 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 < NDVI_s) 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 > NDVI_s) in this scenario, a pixel is considered to be composed of full vegetation cover, therefore, the LSE of vegetation is assigned to the pixel.
• (NDVI_s ≤ NDVI ≤ NDVI_v) 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.

$${\mathsf{\epsilon }}_{\mathsf{\lambda }}={\mathsf{\epsilon }}_{v\mathsf{\lambda }}{P}_v+{\mathsf{\epsilon }}_{s\mathsf{\lambda }}\left(1-P_v\right)+C_{\mathsf{\lambda }}$$

(6)

where, ε_v stands for the emissivity of vegetation, ε_s is the emissivity of soil, P_v 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_{\mathsf{\lambda }}=\left(1-{\mathsf{\epsilon }}_{s\mathsf{\lambda }}\right){\mathsf{\epsilon }}_{v\mathsf{\lambda }}{F}^{\prime }\left(1-P_v\right)$$

(7)

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. LSE values for various terrestrial materials [18].

Terrestrial Material	Soil	Built-Up Areas	Vegetation	Water
Emissivity	0.966	0.962	0.973	0.991

have been used in the study [18].

$$P_v={\left[\frac{NDVI-NDV{I}_{min}}{NDV{I}_{max}-NDV{I}_{min}}\right]}^2$$

(8)

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. Simulated atmospheric parameters for the study.

Sensor Type	Date	Time (UTC)	Atmospheric Transmittance	Upwelling Radiance (W/(m2·sr·µm))	Down-Welling Radiance (W/(m2·sr·µm))
Landsat 5 TM	13 November 2010	14:57	0.89	0.72	1.20
Landsat 5 TM	31 December 2010	14:57	0.89	0.58	0.97
Landsat 7 ETM+	21 November 2010	15:00	0.96	0.18	0.31
Landsat 7 ETM+	23 February 2016	15:09	0.97	0.11	0.20
Landsat 8 TIRS	4 June 2015	15:06	0.87	0.91	1.52
Landsat 8 TIRS	29 May 2013	15:09	0.89	0.72	1.22

.

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\times \left\{10\times 0.6108\times exp\left[\frac{17.27\times \left(T_0-273.15\right)}{237.3+\left(T_0-273.15\right)}\right] \times RH\right\}+0.1679$$

(9)

In Equation (9), W_i stands for atmospheric water vapor, T₀ 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=\frac{BT}{\left\{1+\left[\frac{\lambda ·BT}{\rho }\right]·ln\epsilon \right\}}$$

(10)

where T_s 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⁻² 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=\frac{a_{\mathrm{i}}\left(1-C_{\mathrm{i}}-D_{\mathrm{i}}\right)+\left[b_{\mathrm{i}}\left(1-C_i-D_{\mathrm{i}}\right)+C_i+D_{\mathrm{i}}\right]T_{\mathrm{i}}-D_{\mathrm{i}}{T}_a}{C_{\mathrm{i}}}$$

(11)

where T_S is the Land Surface Temperature (LST), T_i is the brightness temperature from Equation (4), T_a is the effective mean atmospheric temperature, a_i = −67.355351 and b_i = 0.458606. The values of C_i and D_i can be calculated using the Equations (12) and (13), respectively [6].

$$C_i={\mathsf{\epsilon }}_i{\mathsf{\tau }}_i$$

(12)

$$D_{\mathrm{i}}=\left(1-{\mathsf{\tau }}_{\mathrm{i}}\right)\left[1+\left(1-{\mathsf{\epsilon }}_{\mathrm{i}}\right){\mathsf{\tau }}_{\mathrm{i}}\right]$$

(13)

where ε_i is the ground surface emissivity and τ_i is the atmospheric transmittance. T_a, ε_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. Linear relations for the approximation of effective mean atmospheric temperature T_a from the near surface air temperature [6].

Atmospheres	Linear Relations Equations
For USA 1976	Ta = 25.9396 + 0.88045T0
Tropical model	Ta = 17.9769 + 0.91715T0
Mid-latitude Summer	Ta = 16.0110 + 0.92621T0
Mid-latitude winter	Ta = 19.2704 + 0.91118T0

. 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 T_a and T₀ 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 T_a is a linear relation of near surface air temperature T₀. 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_{\lambda }=\left[{\mathsf{\epsilon }}_{\mathsf{\lambda }}{L}_{\mathsf{\lambda }}\left(T_s\right)+\left(1-{\mathsf{\epsilon }}_{\mathsf{\lambda }}\right)L_{down}\right]\mathsf{\tau }+L_{\mathsf{\lambda }up}$$

(14)

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_λ(T_s) is the radiance of a blackbody target of kinetic temperature T_s; 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_λ(T_s) of kinetic temperature T_s can be calculated by Equation (15).

$$L_{\mathsf{\lambda }}\left(T_s\right)=\frac{L_{\mathsf{\lambda }}-L_{\mathsf{\lambda }up}}{{\mathsf{\tau }\mathsf{\epsilon }}_{\mathsf{\lambda }}}-\frac{1-{\mathsf{\epsilon }}_{\mathsf{\lambda }}}{{\mathsf{\epsilon }}_{\mathsf{\lambda }}}{L}_{\mathsf{\lambda }down}$$

(15)

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 K₁ and K₂. To do this, Equation (16) has been used in the study [31,32].

$$T_s=\frac{K_2}{ln\left(\frac{K_1}{L_{\mathsf{\lambda }}\left(T_s\right)}+1\right)}$$

(16)

where T_s is the LST in Kelvin which is changed to degree Celsius (°C) by subtracting 273.15; and K₁ and K₂ 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=\mathsf{\gamma }\left[{\mathsf{\epsilon }}^{-1}\left({\mathsf{\psi }}_1{L}_{sensor}+{\mathsf{\psi }}_2\right)+{\mathsf{\psi }}_3\right]+\mathsf{\delta }$$

(17)

$$\mathsf{\gamma }={\left\{\frac{C_2{L}_{sensor}}{T_{sensor}^2}\left[\frac{{\mathsf{\lambda }}^4}{C_1}{L}_{sensor}+{\mathsf{\lambda }}^{-1}\right]\right\}}^{-1}$$

(18)

$$\mathsf{\delta }=-\mathsf{\gamma }{L}_{sensor}+T_{sensor}$$

(19)

From the equations, T_s stands for LST, T_sensor is the at sensor brightness temperature in Kelvin, λ is the effective wavelength of a thermal infrared band in use, C₁ = 1.19104 × 10⁸ W·m^-2·sr^-1·μm⁴ and C₂ = 14387.7 µm·K. The ψ₁, ψ₂, and ψ₃ atmospheric parameters can be estimated from Equations (20)–(22), respectively. Figure 7 shows the application’s interface for the SCA.

$${\mathsf{\psi }}_1=0.14714w^2-0.15583w+1.1234$$

(20)

$${\mathsf{\psi }}_2=-1.1836w^2-0.3760w-0.52894$$

(21)

$${\mathsf{\psi }}_3=-0.04554w^2+1.8719w-0.39071$$

(22)

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. 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 Name	Latitude	Longitude	Elevation (Meters)	NST (°C) at 12:00	RH (%)	MWA (°C)	PLANCK (°C)	RTE (°C)	SCA (°C)
Mechanic Settlement	45°41′37.040″N	65°09′54.040″W	403	−3.1	99	−4.75	−4.70	−4.55	−5.56
Fundy Park (Alma) Cs	45°36′00.000″N	64°57′00.000″W	42.7	3.6	62	−0.91	−1.29	−0.73	−2.17
Buctouche Cda Cs	46°25′49.006″N	64°46′05.009″W	35.9	−0.9	78	−8.06	−7.64	−7.85	−8.48
Nappan Auto	45°45′34.400″N	64°14′29.200″W	19.8	−7.9	73	−3.46	−3.55	−3.26	−4.42
Saint John A	45°19′05.000″N	65°53′08.050″W	108.8	−0.2	69	−0.29	−0.73	−0.11	−1.62
Parrsboro	45°24′48.000″N	64°20′49.000″W	30.9	2.1	68	−0.29	−0.73	−0.11	−1.62
Gagetown A	45°50′00.000″N	66°26′00.000″W	50.6	−2.2	79	−0.91	−1.29	−0.73	−2.17
Gagetown Awos A	45°50′20.000″N	66°26′59.000″W	50.6	−1	77	−0.29	−0.73	−0.11	−1.62
Summerside	46°26′28.000″N	63°50′17.000″W	12.2	−0.9	78	−3.46	−3.55	−3.26	−4.42
Mechanic Settlement	45°41′37.040″N	65°09′54.040″W	403	−3.1	99	−4.75	−4.70	−4.55	−5.56
Root Mean Square Error (RMSE)						2.67	2.78	2.64	3.17

^NST Near Surface Temperature (Day light saving was observed); ^PLANCK Planck Equation; ^RTE Radiative Transfer Equation; ^SCA Single Channel Algorithm; ^RH Relative Humidity.

and

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 Name	Latitude	Longitude	Elevation (Meters)	NST (°C) at 12:00	RH (%)	MWA (°C)	PLANCK (°C)	RTE (°C)	SCA (°C)
Mechanic Settlement	45°41′37.040″N	65°09′54.040″W	403	15.4	40	15.15	14.97	15.32	13.28
Fundy Park (Alma) Cs	45°36′00.000″N	64°57′00.000″W	42.7	13.1	53	15.70	15.80	15.73	10.86
Buctouche Cda Cs	46°25′49.006″N	64°46′05.009″W	35.9	12.4	39	12.70	12.58	12.90	10.36
Nappan Auto	45°45′34.400″N	64°14′29.200″W	19.8	12.1	52	14.17	14.26	14.30	10.36
Saint John A	45°19′05.000″N	65°53′08.050″W	108.8	12.2	41	14.83	14.25	15.14	13.27
Parrsboro	45°24′48.000″N	64°20′49.000″W	30.9	13.2	48	14.56	14.57	14.70	10.85
Gagetown A	45°50′00.000″N	66°26′00.000″W	50.6	11.6	36	10.11	10.27	10.40	8.37
Gagetown Awos A	45°50′20.000″N	66°26′59.000″W	50.6	12.0	36	10.96	10.82	11.30	9.86
Summerside	46°26′28.000″N	63°50′17.000″W	12.2	13.6	40	12.63	12.55	12.88	10.36
Root Mean Square Error (RMSE)						1.64	1.58	1.68	2.33

^NST Near Surface Temperature (Day light saving was observed); ^PLANCK Planck Equation; ^RTE Radiative Transfer Equation; ^SCA Single Channel Algorithm; ^RH Relative Humidity.

show the values obtained from the algorithms and the near surface temperatures for Landsat 5 TM.

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 Name	Latitude	Longitude	Elevation (Meters)	NST (°C) at 12:00	RH (%)	MWA (°C)	PLANCK (°C)	RTE (°C)	SCA (°C)
Mechanic Settlement	45°41′37.040″N	65°09′54.040″W	403	−6.8	78	−11.75	−11.52	−11.43	−12.34
Fundy Park (Alma) Cs	45°36′00.000″N	64°57′00.000″W	42.7	−1.9	58	−5.50	−5.53	−5.21	−6.39
Buctouche Cda Cs	46°25′49.006″N	64°46′05.009″W	35.9	−4.4	64	−9.62	−9.48	−9.31	−10.31
Nappan Auto	45°45′34.400″N	64°14′29.200″W	19.8	−4.4	70	−8.23	−8.14	−7.92	−8.98
Saint John A	45°19′05.000″N	65°53′08.050″W	108.8	−3.6	55	−6.18	−6.18	−5.88	−7.03
Parrsboro	45°24′48.000″N	64°20′49.000″W	30.9	−2.6	60	N/A	N/A	N/A	N/A
Gagetown Awos A	45°50′20.000″N	66°26′59.000″W	50.6	−2.9	55	−3.51	−3.62	−3.22	−4.49
Summerside	46°26′28.000″N	63°50′17.000″W	12.2	−3	68	−8.23	−8.14	−7.92	−8.98
Mechanic Settlement	45°41′37.040″N	65°09′54.040″W	403	−6.8	78	−11.75	−11.52	−11.43	−12.34
Root Mean Square Error (RMSE)						4.03	3.94	3.75	4.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).

and

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 Name	Latitude	Longitude	Elevation (Meters)	NST (°C) at 12:00	RH (%)	MWA (°C)	PLANCK (°C)	RTE (°C)	SCA (°C)
Mechanic Settlement	45°41′37.040″N	65°09′54.040″W	403	−8.1	35	−14.46	−14.31	−14.04	−15.12
Fundy Park (Alma) Cs	45°36′00.000″N	64°57′00.000″W	42.7	−6.9	50	−10.17	−10.15	−9.77	−10.98
Buctouche Cda Cs	46°25′49.006″N	64°46′05.009″W	35.9	−8.9	45	−5.40	−5.53	−5.03	−6.39
Nappan Auto	45°45′34.400″N	64°14′29.200″W	19.8	−7.9	68	N/A	N/A	N/A	N/A
Parrsboro	45°24′48.000″N	64°20′49.000″W	30.9	−5.9	49	−7.41	−7.48	−7.03	−8.33
Moncton lnternational Airport	46°06′44.000″N	64°40′43.000″W	70.7	−8	44	−4.08	−4.25	−3.72	−5.12
Gagetown A	45°50′00.000″N	66°26′00.000″W	50.6	−8.4	47	N/A	N/A	N/A	N/A
Summerside	46°26′28.000″N	63°50′17.000″W	12.2	−7.3	55	N/A	N/A	N/A	N/A
Root Mean Square Error (RMSE)						3.68	3.58	3.61	3.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).

show the values obtained from the algorithms and the near surface temperatures for Landsat 7 ETM+.

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 Name	Latitude	Longitude	Elevation (m)	NST (°C) at 11:00 am	RH (%)	MWA (°C)	PLANCK (°C)	RTE (°C)	SCA (°C)
Mechanic Settlement	45°41′37.040″N	65°09′54.040″W	403	13.1	57	9.27	9.65	9.43	8.15
Fundy Park (Alma) Cs	45°36′00.000″N	64°57′00.000″W	42.7	12.4	56	13.42	13.36	13.51	11.33
Buctouche Cda Cs	46°25′49.006″N	64°46′05.009″W	35.9	15.1	58	16.66	16.14	16.73	14.20
Nappan Auto	45°45′34.400″N	64°14′29.200″W	19.8	15.4	63	14.84	14.40	15.01	13.18
Parrsboro	45°24′48.000″N	64°20′49.000″W	30.9	15.9	42	11.70	12.10	11.69	9.17
Summerside	46°26′28.000″N	63°50′17.000″W	12.2	13.9	62	13.44	13.55	13.45	10.82
Moncton International Airport	46°06′44.000″N	64°40′43.000″W	70.7	16.9	44	17.74	17.55	17.84	14.25
Root Mean Square Error (RMSE)						2.30	2.07	2.28	3.65

^NST Near Surface Temperature; ^PLANCK Planck Equation; ^RTE Radiative Transfer Equation; ^SCA Single Channel Algorithm; ^RH Relative Humidity.

and

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 Name	Latitude	Longitude	Elevation (Meters)	NST (°C) at 11:00 am	RH (%)	MWA (°C)	PLANCK (°C)	RTE (°C)	SCA (°C)
Mechanic Settlement	45°41′37.040″N	65°09′54.040″W	403	17	27	22.24	21.80	23.07	18.85
Fundy Park (Alma) Cs	45°36′00.000″N	64°57′00.000″W	42.7	18.1	38	18.41	18.41	19.31	15.53
Buctouche Cda Cs	46°25′49.006″N	64°46′05.009″W	35.9	20.1	35	23.07	22.28	23.98	20.62
Nappan Auto	45°45′34.400″N	64°14′29.200″W	19.8	19.3	35	21.66	21.06	22.59	19.29
Parrsboro	45°24′48.000″N	64°20′49.000″W	30.9	16.2	49	22.04	21.60	22.88	18.77
Moncton International Airport	46°06′44.000″N	64°40′43.000″W	70.7	20.2	41	29.74	28.78	30.24	23.70
Gagetown A	45°50′00.000″N	66°26′00.000″W	50.6	20	31	15.36	17.58	16.34	13.32
Gagetown Awos A	45°50′20.000″N	66°26′59.000″W	50.6	20	31	22.02	21.74	22.79	18.06
Summerside	46°26′28.000″N	63°50′17.000″W	12.2	16.1	50	20.26	19.67	21.26	18.65
Root Mean Square Error (RMSE)						4.83	4.33	5.35	3.06

^NST Near Surface Temperature; ^PLANCK Planck Equation; ^RTE Radiative Transfer Equation; ^SCA Single Channel Algorithm; ^RH Relative Humidity.

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 m² 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.