## 1 Introduction
### Reference Parameters and Datasets
The spectral radiance was calculated as follows:  
\[L_{\gamma}=\frac{(L_{max}-L_{min})*Q_{cal}}{(Q_{calmax}-Q_{calmin})}+L_{min}-O_ {i} \tag{2}\]  
\[L_{\gamma}=M_{L}*Q_{cal}+A_{L} \tag{3}\]  
where \(L_{max}\) is the maximum radiance (Wm\({}^{-2}\)sr\({}^{-1}\)um\({}^{-1}\)), \(L_{min}\) is the minimum radiance (Wm\({}^{-2}\)sr\({}^{-1}\)um\({}^{-1}\)), \(Q_{cal}\) is the value of the pixel, \(Q_{cal}\) is the maximum value of pixels, \(Q_{cal}\) is the minimum value of pixels, and \(O_{i}\) is the correction value for Band 10. These values can be obtained from the metadata file of each raster. The metadata values were the same for all years. Another method is to determine the top of the atmosphere, where \(M_{L}\) is the band-specific multiplicative rescaling factor from metadata, and \(A_{L}\) is the band-specific additive rescaling factor from metadata (Table 1). The brightness temperature conversion (From \(Ly\) to BT) was then calculated as follows:  
\[\text{BT}=\frac{K_{2}}{\ln[\left(\frac{K_{1}}{Dy}\right)+1]}-273.15 \tag{4}\] \begin{table}
\begin{tabular}{c c c} \hline \hline
**Variables** & **Description** & **Values** \\ \hline \(K_{1}\) & Thermal constants Band 10 & 774.8853 \\ \(K_{2}\) & 1321.0789 \\ \(M_{L}\) & Band-specific multiplicative rescaling factor & \(3.34\times 10^{-4}\) \\ \(A_{L}\) & Band-specific additive rescaling factor & 0.1 \\ \(L_{\text{max}}\) & Maximum value or radiance, Band 10 & 22.0018 \\ \(L_{\text{min}}\) & Minimum value or radiance, Band 10 & 0.10033 \\ \(Q_{cal}\)max & Maximum values of quantized calibration, Band 10 & 65,535 \\ \(Q_{cal}\)min & Minimum value of quantized calibration, Band 10 & 1 \\ \hline \hline \end{tabular}
\end{table} 
Table 1: Metadata of the satellite images from 2013–2022 for Landsat 8–9 OLI-TRIS. where \(K_{1}\) and \(K_{2}\) are the thermal constants of Band 10. These values can be obtained from the metadata file attached to the raster. \(L\) is the result of Equation (1). To obtain the results in Celsius, we must then subtract 273.15. To calculate NDVI, we used two methods: (1) MOD13Q1 (16 days) at a 250 m spatial resolution and a valid range from \(-2\) to 10, in which we applied a 0.0001 scale factor as recommended by LP DAAC [23], and (2) using Landsat 8-9 OLI with Bands 6 and 4:  
\[\text{NDVI}=\text{MOD13Q1-NDVI}^{*}\;0.0001 \tag{5}\]  
\[\text{NDVI}=\frac{(\text{Band6}-\text{Band4})}{(\text{Band6}+\text{Band4})} {=\frac{NIR-RED}{NIR+RED}} \tag{6}\]  
where \(NIR\) is the near-infrared band value of a pixel (Band 6), and \(RED\) is the red band value of the same pixel (Band 4). It was noticed that NDVI estimates from Landsat were lower than NDVI values obtained from MODIS. Proportional vegetation (Pv) was computed as follows:  
\[\text{Pv}=\frac{\text{NDVI}-\text{NDVI}_{min}}{\text{NDVI}_{max}-\text{NDVI}_ {min}} \tag{7}\]  
where \(\text{NDVI}_{min}\) is the minimum value and \(\text{NDVI}_{max}\) is the maximum value of NDVI obtained by either Equation (5) or Equation (6). On a global scale, it is proposed to use values of \(\text{NVDI}_{min}\) = 0.2 and \(\text{NVDI}_{max}\) = 0.5, but in some cases, these values tend to underestimate \(\text{NVDI}\)\(min\) and \(max\) values. Due to the variation between global regions, we used the NDVI range for each year. Land surface emissivity (LSE) was then calculated:  
\[\text{LST}=\frac{\text{BT}}{\left(1+\left[\left(\frac{\lambda\text{BT}}{ \rho}\right)*\text{Ln}(\text{E})\right]\right.} \tag{8}\]  
where \(\lambda\) is the average wavelength of Band 10 equal to 0.00115; and \(\rho\) is equal to 1.4388 \(\times\) 10\({}^{-2}\) mK, which derives from \(\left(\rho=\hbar*c/\sigma\right)\), where \(\hbar\) is [PERSON]'s constant, \(c\) is the velocity of light, and \(\sigma\) is the [PERSON] constant. BT and E are the results of steps 2 and 5, respectively. The result of LST (Equation (8)) is given in Celsius.