2.2. Method
2.2.1. Data sources and preprocessing
In this study, Landsat TM5 satellite images (1989) and Landsat 8 OLI/TIRS ones (2019) were extracted by using GEE system to prepare LST and LULC maps (Table 1). In addition, this system (Shiflett et al., 2017) was applied to calculate the indices of land surface cover such as NDVI, NDWI, NDBI, LSM, and albedo.
The correlation between the intended parameters with LST was analyzed in QGIS software through employing the GWR. Further, the relationship between land use map and LST was obtained. Figure 2 displays a flowchart of stages in the present study.
Table 1
Specifications of the utilized satellite images for 1989 and 2019
Satellite | Path/ Row | Type bands | Dataset | E.C.T* | Images Date |
Landsat TM | 157/38 157/39 | Spectral bands | C01/T1_SR | 5.55 (16-day) | 1989-01-16, 1989-03-05, 1989-04-06, 1989-04-22, 1989-05-08, 1989-05-24, 1989-06-01, 1989-06-25, 1989-08-12, 1989-08-28, 1989-09-13, 1989-10-15 |
Thermal bands | C01/T1_TOA |
Landsat OLI/TIRS | 157/38 157/39 | Spectral bands | C01/T1_SR | 6.25 (16-day) | 1989-01-16, 1989-03-05, 1989-04-06, 1989-04-22, 1989-05-08, 1989-05-24, 1989-06-01, 1989-06-25, 1989-08-12, 1989-08-28, 1989-09-13, 1989-10-15 |
Thermal bands | C01/T1_TOA |
*Equatorial Crossing Time |
2.2.2. Selection of appropriate environmental indices
Primary indices were selected according to the previous studies, some of which were eliminated by considering factors like the conditions of the intended area to reduce the number of the indices and simplify process. For example, most parts of the area are plain and possess low elevation. Accordingly, no significant difference was found between their aspect, slope, and elevation, leading to the removal of these variables. To make better decision on the primary parameters, the correlation between different indices was obtained through using GWR based on the recommendation of the various studies (Alibakhshi et al., 2020; Kashki et al., 2021; Xu et al., 2021; Liu et al., 2022) (Table 2). Then, the factors with high correlation were eliminated so that NDMI was removed due to its high correlation with NDWI and NDVI. Finally, NDVI, NDWI, NDBI, LSM, and albedo were selected to evaluate the relationship between LST and environmental parameters.
Table 2
Correlation between the indices from the previous studies
Factors | 1989-Annual Average | 2019-Annual Average |
Dependent Factor | Independent Factor | AICc | R2 | Ajusted R2 | sigma | AICc | R2 | Ajusted R2 | sigma |
Aspect | Albedo | 7514.90 | 0.56 | 0.45 | 4.10 | 7568.94 | 0.58 | 0.47 | 4.47 |
Elevation | 3220.91 | 0.78 | 0.72 | 5.73 | 3220.91 | 0.78 | 0.72 | 5.73 |
LSM | 7798.40 | 0.49 | 0.35 | 5.33 | 7444.00 | 0.43 | 0.28 | 3.96 |
NDBI | 1390.25 | 0.55 | 0.43 | 0.05 | 1555.67 | 0.39 | 0.28 | 0.04 |
NDVI | 1265.71 | 0.57 | 0.45 | 0.06 | 1090.35 | 0.54 | 0.42 | 0.07 |
NDWI | 1053.25 | 0.57 | 0.46 | 0.07 | 784.25 | 0.57 | 0.45 | 0.10 |
NMDI | 1469.02 | 0.63 | 0.53 | 0.05 | 1393.25 | 0.36 | 0.27 | 0.05 |
Slope | 1976.53 | 0.06 | 0.04 | 1.74 | 1976.53 | 0.06 | 0.04 | 1.74 |
Elevation | Albedo | 3191.74 | 0.76 | 0.72 | 5.69 | 3130.14 | 0.81 | 0.76 | 5.26 |
LSM | 3174.35 | 0.79 | 0.74 | 5.51 | 3121.48 | 0.81 | 0.77 | 5.22 |
NDBI | 3122.15 | 0.81 | 0.76 | 5.23 | 3104.47 | 0.82 | 0.77 | 5.13 |
NDMI | 3163.47 | 0.79 | 0.75 | 5.45 | 3240.00 | 0.69 | 0.68 | 6.15 |
NDVI | 3194.84 | 0.78 | 0.73 | 5.62 | 3138.25 | 0.80 | 0.76 | 5.32 |
NDWI | 3177.83 | 0.79 | 0.74 | 5.54 | 3136.77 | 0.79 | 0.75 | 5.39 |
Slope | 3186.57 | 0.79 | 0.74 | 5.55 | 3186.57 | 0.79 | 0.74 | 5.55 |
Slope | Albedo | 7510.76 | 0.57 | 0.45 | 4.22 | 7567.76 | 0.56 | 0.462 | 4.56 |
LSM | 7811.78 | 0.47 | 0.33 | 5.71 | 7418.45 | 0.44 | 0.31 | 3.87 |
NDBI | 1341.24 | 0.49 | 0.37 | 0.05 | 1569.26 | 0.37 | 0.29 | 0.04 |
NMDI | 1460.29 | 0.61 | 0.52 | 0.05 | 1405.55 | 0.35 | 0.28 | 0.05 |
NDVI | 1253.47 | 0.55 | 0.44 | 0.06 | 1103.34 | 0.54 | 0.43 | 0.07 |
NDWI | 1033.02 | 0.55 | 0.43 | 0.08 | 824.35 | 0.59 | 0.49 | 0.09 |
Albedo | LSM | 6804.41 | 0.89 | 0.87 | 2.06 | 6834.03 | 0.90 | 0.88 | 2.12 |
NDBI | 7218.89 | 0.75 | 0.70 | 3.12 | 7234.01 | 0.78 | 0.73 | 3.17 |
NMDI | 7208.40 | 0.76 | 0.70 | 3.09 | 7491.89 | 0.53 | 0.51 | 4.29 |
NDVI | 7298.86 | 0.71 | 0.65 | 3.38 | 7229.70 | 0.78 | 0.73 | 3.165 |
NDWI | 7257.90 | 0.73 | 0.67 | 3.26 | 7197.31 | 0.78 | 0.74 | 3.10 |
LSM | NDBI | 7180.49 | 0.85 | 0.81 | 3.01 | 6745.44 | 0.86 | 0.82 | 1.94 |
NDMI | 7555.79 | 0.68 | 0.61 | 4.06 | 6940.30 | 0.73 | 0.72 | 2.24 |
NDVI | 7459.32 | 0.73 | 0.67 | 3.97 | 6873.38 | 0.81 | 0.77 | 1.25 |
NDWI | 7449.53 | 0.74 | 0.68 | 3.98 | 6888.44 | 0.79 | 0.76 | 2.28 |
NDMI | 1823.25 | 0.80 | 0.75 | 0.03 | 2791.85 | 0.93 | 0.93 | 0.01 |
NDBI | NDVI | 2023.58 | 0.86 | 0.83 | 0.03 | 2384.36 | 0.89 | 0.86 | 0.02 |
NDWI | 2016.61 | 0.86 | 0.83 | 0.03 | 2217.73 | 0.83 | 0.80 | 0.86 |
NDVI | NDMI | 1912.36 | 0.87 | 0.84 | 0.03 | 1351.54 | 0.63 | 0.62 | 0.06 |
NDWI | 3051.47 | 0.98 | 0.98 | 0.01 | 2791.33 | 0.98 | 0.97 | 0.01 |
NDMI | NDWI | 1914.45 | 0.92 | 0.90 | 0.03 | 927.84 | 0.56 | 0.54 | 0.09 |
2.2.3. LST calculation
The mean annual LST (MALST) was computed through converting spectral radiance, recovering brightness temperature, and calculating ratio vegetation index and surface emissivity (Eq. 6).
The digital values are converted to spectral radiance based on the spectrum radiation reference as follows (Eq. 1).
$${L}_{\lambda }=Gain\times DN+Bias$$
1
Where \({L}_{\lambda }\) indicates spectral radiance and Gain demonstrates band specific multiplicative rescaling factor. Additionally, DN reflects digital numbers for each pixel and Bias reveals add-on rescaling factor (Yue et al., 2019).
Planck equation is used to convert the amount of spectral radiance to brightness temperature (Eq. 2), in which Tsensor represents the brightness temperature of sensor (K) (Wukelic et al., 1989). \({K}_{1}\) and \({\text{K}}_{2}\) indicate the calibration constants for TM, ETM+/ OLI sensor thermal band. Regarding TM, ETM+, and OLI, K1 equals 607.76, 666.09, and 774.89 mWcm-2, while K2 is 1260.56, 1282.71, and 1321.08 K, respectively.
$${\text{T}}_{\text{s}\text{e}\text{n}\text{s}\text{o}\text{r}}= \frac{{\text{K}}_{2}}{ln\left({K}_{1}/{L}_{\lambda }+1\right)}$$
2
Further, NDVI is used to obtain ratio vegetation index by using Eq. 3 (Artis & Carnahan, 1982; Griend & Owe, 1993).
(3)\({P}_{v} ={\left.\left( \frac{NDVI-{NDVI}_{min}}{{NDVI}_{max}-{NDVI}_{min}}\right.\right)}^{2}\)
where NDVImin and NDVImax are respectively considered as the least NDVI in bare soil pixel and highest value in the vegetation-containing pixel.
Eq. 4 is applied for the conditional estimation of surface emissivity.
$$\epsilon ={E}_{v\lambda }{P}_{v}+{E}_{s\lambda }\left(1-{P}_{v}\right)+{C}_{\lambda }$$
4
in this respect, C demonstrates surface roughness factor, and Ev and Es are the emissivity of vegetation cover and soil, respectively (Artis & Carnahan, 1982).
Eq. 5 reveals the conditions for emissivity determination by using NDVI (Griend & Owe, 1993).
$$\epsilon =\left\{\begin{array}{c}0.995 NDVI\le 0 \\ 0.970 0<NDVI\le 0.157 \\ 1.0094+0.047\text{ln}NDVI 0.157<NDVI\le 0.727\\ 0.986 NDVI>0.727 \\ \end{array}\right.$$
5
where \(\lambda\) shows emitted wavelength (11.435 µm for Landsats 5 and 7; 10.9 µm for the band 10 of Landsat 8). Furthermore, \(\rho\) is equal to 1.438 × 10− 2 mK (constant), and \(\epsilon\) reflects land surface emission (Artis & Carnahan, 1982).
$$LST= \frac{{T}_{sensor}}{ln\epsilon \left(\lambda \times {T}_{sensor}/\rho \right)+1}$$
6
Given that the climate changes caused by temporal variations only affect absolute temperature, the normalization method utilized by Xu and Chen (2004) was employed to better evaluate the changing pattern of LST distribution. Accordingly, this approach allowed assess LST under the different climate conditions induced by temporal alterations through changing LST scale between zero and one (Eq. 7).
$${NLST}_{i}=\frac{{LST}_{i}-{LST}_{min}}{{LST}_{max}-{LST}_{min}}$$
7
where LSTi and NLSTi refer to LST and its normalized value in the pixel i, respectively. The subscripts of min and max represent the minimum and maximum LST in each image, respectively. To examine the trend of LST changes in the period under study, normalized land surface temperature was classified according to the mean and standard deviation. As outlined in Table 3, LSTmean and LSTSTD are the mean and standard deviation of LST in a normalized land surface temperature image (Firozjaei et al., 2018).
Table 3
Classification of LST map (Firozjaei et al., 2018)
LST | Class range |
Very low | LST ≤ LST mean − 1.5LSTSTD |
low | LST mean − 1.5LSTSTD < LST ≤ LST mean - LSTSTD |
Moderate | LST mean - LSTSTD < LST ≤ LST mean + LSTSTD |
High | LST mean + LSTSTD < LST ≤ LST mean + 1.5LSTSTD |
Very High | LST > LST mean + 1.5LSTSTD |
2.2.4. Environmental index calculation
The five environmental indices influencing LST were utilized in this study.
NDVI:
This index reflects the trend of vegetation cover variations. Given the dimensionlessness of this index, its values range − 1 to + 1. Practically, the amounts below 0.1 are related to land without water, while the higher levels are associated with agricultural activities and forest regions (Eq. 8).
NDVI = \(\frac{\text{N}\text{I}\text{R}-\text{R}}{\text{N}\text{I}\text{R}+\text{R}}\) (8)
where red and near-infrared bands are denoted as R and NIR, respectively (Rouse et al., 1974).
NDWI:
It can be applied as a strong normalized difference spectral index to specify the correlation with LST, which is determined by green and near-infrared (NIR) bands. As shown in Table 3, the second and fourth bands are respectively considered as green and NIR bands in both TM and ETM, while bands 3 and 5 are respectively used in this regard for OLI/TIRS data. This parameter varies between − 1 and 1 so that the negative values reveal built-up and bare lands without water, while the positive NDWI indicates water bodies and vegetation cover (Eq. 9) (McFeeters, 2013).
NDWI = \(\frac{\text{G}\text{r}\text{e}\text{e}\text{n}-\text{N}\text{I}\text{R}}{\text{G}\text{r}\text{e}\text{e}\text{n}+\text{N}\text{I}\text{R}}\) (9)
LSM:
This variable is among the most important environmental indices. Soil moisture at 1–2 m above land surface has been widely identified as a key factor in the numerous environmental studies such as those in the field of meteorology, hydrology, agriculture, and climate changes. Thus, it should be estimated in the climate, hydrological, and agricultural research, the monitoring of which acts a crucial role in supervising and forecasting flood, drought, and other climatic phenomena (Hu & Xu, 2018; Sun & Pinker, 2004).
$${\text{L}\text{S}\text{M}}_{\text{T}\text{M}}=0.0315{{\rho }}_{\text{B}\text{l}\text{u}\text{e}}+0.2021{{\rho }}_{\text{G}\text{r}\text{e}\text{e}\text{n}}+0.3102{{\rho }}_{\text{R}\text{e}\text{d}}+0.1594{{\rho }}_{\text{N}\text{I}\text{R}}-0.6806{{\rho }}_{\text{S}\text{W}\text{I}\text{R}1}-0.6109{{\rho }}_{\text{S}\text{W}\text{I}\text{R}2}$$
10
$${\text{L}\text{S}\text{M}}_{\text{O}\text{L}\text{I}}=0.1511{{\rho }}_{\text{B}\text{l}\text{u}\text{e}}+0.1972{{\rho }}_{\text{G}\text{r}\text{e}\text{e}\text{n}}+0.3283{{\rho }}_{\text{R}\text{e}\text{d}}+0.3407{{\rho }}_{\text{N}\text{I}\text{R}}-0.7117{{\rho }}_{\text{S}\text{W}\text{I}\text{R}1}-0.4559{{\rho }}_{\text{S}\text{W}\text{I}\text{R}2}$$
11
In this regard, Blue, Green, Red, NIR, SWIR1, and SWIR2 illustrate blue, green, red, and near-, middle-, and far-infrared bands, respectively ( Baig et al., 2014).
NDBI:
It is one of the indices for obtaining human land use intensity, human effect, and built-up uses (Eq. 12) (Hu & Xu, 2018), in which near- and middle-infrared bands are demonstrated as NIR and SWIR1, respectively.
NDBI = \(\frac{SWIR1-\text{N}\text{I}\text{R}}{\text{S}\text{W}\text{I}\text{R}1+\text{N}\text{I}\text{R}}\) (12)
Albedo:
Surface albedo is defined as the ratio of the reflected light to the total energy incident on a surface in a hemisphere space. It plays a fundamental role in global climate changes, which influences land energy level. The variations in the physicochemical properties and spatial structure of local objects in turn alter surface albedo (Alibakhshi et al., 2020). Currently, this index is extensively applied in modeling for radiation-energy balance, numerical weather prediction, atmospheric circulation, and land surface processes (Eq. 13) (Bonafoni et al., 2017; Trlica et al., 2017). Regarding the factor, blue, green, red, and near-, middle-, and far-infrared bands are respectively denoted as Blue, Green, Red, NIR, SWIR1, and SWIR2.
$${\alpha }_{toa}=0.356{\rho }_{Blue}+0.130{\rho }_{Red}+0.373{\rho }_{NIR}+0.085{\rho }_{SWIR1}+0.072{\rho }_{SWIR2}-0.0018$$
13
2.2.5. Classification of LULC map
Yousefi et al. (2015) reported non-parametric Support Vector Machine (SVM) method as the best approach to produce land use map in the arid regions. This method was introduced by Vapnik and Chervonenkis (1971) as a linear classifier for the first time. In the present study, it was used to prepare the LULC map related to the lands of Zabol Basin in the GEE system in the five categories of irrigated and rainfed agriculture, bare land, water body, and built-up areas for 1989 and 2019. The accuracy of land use classification in the basin under study was evaluated by using Overall Accuracy (OA) and Kappa Coefficient (KC) as assessment criteria based on the widely applicable approach of error matrix calculation. The mathematical model of OA and KC can be expressed as follows (Eqs. 14 & 15).
OA = \(\frac{\sum _{i=1}^{q}nii}{n}\) * 100 (14)
Kappa = \(\frac{\sum _{i=1}^{q}nii-\sum _{i=1}^{q}ni+n+i}{n2- \sum _{i=1}^{q}ni+n+i}\) (15)
where q reflects the number of classes, n refers to the total number of intended pixels, and \(\sum nii\) indicates the sum of main diagonal elements in error matrix. Further, n + reveals the margin of rows, and n + i demonstrates the marginal sum of the columns in the matrix (Rousta et al., 2018).
2.2.6. Correlation between LST and environmental indices
During the recent years, an extensive range of various regression techniques has been utilized to predict and model environmental indices (Taghadosi et al., 2019). The GWR, proposed by Brunsdon et al. (2010), is known as a potent method for modeling the location-dependent data at local level. Furthermore, this approach locally acts in the basin, in which each observation point is weighted by its distance from reference point and the coefficients are not considered in the different fixed or same locations (Wheeler & Páez, 2010). Regarding the present study, GWR method was applied to assess the relationship between LST and intended indices. The regression is determined by using Eq. 16.
$$y={\beta }_{0}+\sum _{j=1}^{n}{\beta }_{j}{x}_{i}+\epsilon$$
16
where n represents the number of independent variables, and y and \({x}_{i}\) illustrate dependent and independent parameters, respectively. Additionally, \({\beta }_{0}\), \({\beta }_{1}\), and \(\epsilon\) denote intercept, coefficient, and error, respectively. Eq. 17 can be used for the local estimation of this regression.
$${y}_{i}={\beta }_{0}\left({u}_{i},{v}_{i}\right)+\sum _{j=1}^{n}{\beta }_{j}({u}_{i},{v}_{i}){x}_{ij}+{\epsilon }_{i}$$
17
in which \(\left({u}_{i}{,v}_{i}\right)\) specifies the location of data i, \({\beta }_{0}\left({u}_{j},{v}_{j}\right)\) is intercept, and \({y}_{i}\) demonstrates independent variable. The \({\beta }_{j}\left({u}_{i},{v}_{i}\right)\) indicates the value of parameter j and \({\epsilon }_{i}\) is random in the ith location.
Further, the data weight is obtained according to the distance from the location i. The parameter estimation matrix for i is as Eq. 18.
$$\beta \left({u}_{i},{v}_{i}\right)={⌈{X}^{T}W\left({u}_{i},{v}_{i}\right)X⌉}^{-1}{X}^{T}W\left({u}_{i},{v}_{i}\right)Y$$
18
where \(W\left({u}_{i},{v}_{i}\right)\) represents a spatial weight matrix, and X and Y reflect dependent and independent variables, respectively.
Gaussian function is applied to compute weighted operator based on the Eq. 19, in which \({W}_{ij}\) demonstrates the data weight observed in the location j for determining dependent parameter in the ith location, and h is considered as bandwidth. In this respect, the weight of location j declines by distancing from location i (Mirchooli et al., 2020).
$${W}_{ij}={e}^{{-0.5\left({d}_{ij}/h\right)}^{2}}$$
19
Furthermore, the validity or efficiency of multivariate regression models was compared by using the coefficient of determination (R2), as well as Akaike information criterion (AIC). The R2 indicates the percentage of the dependent variable variance explained by an independent one. The numerical amount of this coefficient ranges between 0–1 so that zero reflects that an independent parameter has no role in estimating the dependent variable, while one reveals the estimation of 100% of dependent parameter variance by the independent one. The AIC is used to measure the relative efficiency of model, balancing the accuracy and complexity of the model. A small amount of this criterion exhibits the closeness of the level estimated by the model to the observation value or ground truth (Fotheringham et al., 2002).
2.2.7. Spatial autocorrelation
In this study, Moran's index (Moran's I) was determined in Geoda software to assess spatial autocorrelation. This index was introduced by Moran in 1948, which is divided into global and local spatial autocorrelation. The global Moran's I (Eq. 20) captures the spatial pattern according to the location of phenomena, and possesses clustered, random, and dispersed patterns. Another one, local Moran's I, (Eq. 21) is utilized for presenting the Local Indicator of Spatial Analysis (LISA).
$${\text{I}}_{\text{g}}=\frac{N{\sum }_{i}{\sum }_{j}{w}_{ij}\left({x}_{i}-\mu \right)\left({x}_{j}-\mu \right)}{\left({\sum }_{\text{i}}\sum _{\text{j}}{\text{w}}_{\text{i}\text{j}}\right)\sum _{\text{i}}{\left({\text{x}}_{\text{i}}-{\mu }\right)}^{2}}$$
20
where xi indicates the amount of feature in location i and n is the total number of grades in the area under study. In addition, Wij denotes matrix weight (i = 1, 2, 3,… and j = 1, 2, 3,…), which equals zero if i and j are adjacent. Moran's I is equal to + 1 and − 1 in the cases of positive and negative spatial autocorrelation, respectively. However, zero shows the lack of spatial autocorrelation.
LISA can reflect the continuity of local spatial relations by analyzing Moran's I value in each spatial unit (Anselin, 2010).
$${I}_{\text{l}}=\frac{{\text{x}}_{\text{i}}-{\mu }}{{\sum }_{\text{i}}{\left({\text{x}}_{\text{i}}-{\mu }\right)}^{2}}\sum _{\text{j}}{\text{w}}_{\text{i}\text{j}}\left({\text{x}}_{\text{j}}-{\mu }\right)$$
21
Further, local Moran's I is applied to identify hot and cold spots based on comparing with neighbor samples, the negative levels of which refer to data dispersion. The positive local Moran's I exhibits spatial autocorrelation and clustered pattern in the data distribution. In this respect, high-high and low-low clusters demonstrate a high and a low value in a neighborhood of high (hot spots) and low values (cold points), respectively. Furthermore, high-low clusters represent the non-clustered distribution of data so that a high value is surrounded by low values, while low-high ones are related to a low value surrounded by high value feature (Wang et al., 2020).