LST is calculated using Landsat TM 5 band-6 and Landsat 8 OLI images bands 10 and 11. The NDVI is calculated using red band 3 for TM, and 4 for OLI, and near-infrared band 4 for TM, band 5 for OLI. The gathered Landsat data of two path/row series were mosaic using ERDAS Imagine 14, and Arc GIS 10.3 were used for others analysis and layout. To execute remote sensing-based indices and provide the required output, the pictures were georeferenced to the UTM 45° N WGS-84 datum. Calculating spatial indices and assessing and defining the link between LST and other land use indices are the two stages of the current study's strategy. The methodological procedure for the study has been summarized in the diagram below (Fig. 2).

## The Derivation of Land Surface Temperature (LST)

Landsat TM 5 images for the years 2001 and 2011 and Landsat 8 OLI (Operational Land Imager) imagery for the year 2021 with 0% cloud cover were used to calculate the land surface temperature (LST) of the Raiganj urban area. Landsat TM 5 imagery gives 7 bands thermal bands while Landsat 8 OLI imagery provides 11 bands, with band 10 and band 11 being thermal bands (Equations 1 to 6). The following are the procedures involved in turning digital values into spectral radiance when calculating temperature from TM and OLI images:

Step 1: Landsat TM 5 band 6 digital values to spectral radiances conversion.

The following formula was used to convert band 6 digital values into radiance values (Lλ) (Landsat Project Science Office 2002) \(\text{L}{\lambda }= \frac{\text{L}\text{M}\text{A}\text{X}{\lambda }-\text{L}\text{M}\text{I}\text{N}{\lambda }}{\text{Q}\text{C}\text{A}\text{L}\text{M}\text{A}\text{X}-\text{Q}\text{C}\text{A}\text{L}\text{M}\text{I}\text{N} }\times \begin{array}{cc}\left(\text{Q}\text{C}\text{A}\text{L}- \text{Q}\text{C}\text{A}\text{L}\text{M}\text{I}\text{N}\right) + & \text{L}\text{M}\text{I}\text{N}{\lambda }\end{array}\) (1)

Here, Lλ is the atmospherically corrected cell value as the radiance, QCAL is the digital image value, LMINλ is the spectral radiance scaled to QCALMIN, LMAXλ is the spectral radiance scaled to QCALMAX, and QCALMIN is the minimum quantization calibration. The radiance pixel value (usually 1), and QCALMAX are the maximum values of quantized calibrated pixels (usually 255).Landsat8 OLI spectral radiances were converted from band 10 and 11 digital data.

The following formula was used to convert the Landsat OLI band 10 and 11 digital data.

$$\text{L}{\lambda }=ML\times QCAL+AL$$

2

Where Lλ is the spectral radiance at the top of the atmosphere, ML denotes the radiance multi-band X, AL denotes the radiance add band X, QCAL denotes the quantized and calibrated standard product pixel value, and X denotes the band number.The band-specific multiplicative rescaling factor ML and the band-specific additive rescaling factor AL are obtained from the metadata file (MTL file).

Step 2: At-satellite brightness spectral radiance temperatures conversion, emissivity modifications were added to radiant temperature based on the land cover nature following (Artis and Carnahan ,1982).

$$T=\frac{{K}_{2}}{Ln\left(\frac{{\kappa }_{1}}{L{\lambda }}+1\right)}-273\cdot 15$$

3

Where T is the at-satellite brightness temperature in Kelvin (K), Lλ is the at-satellite radiance in W/(m2srµm), and \({K}_{1}\) and \({K}_{2 }\)are the thermal calibration constants in W/(m2sr µm), respectively. (For Landsat-5 TM,\({K}_{1}\)= 607.76, \({K}_{2}\) = 1260.56 for band 6, and for Landsat 8 OLI,\({K}_{1}\) for band 10 and 11 is 774.8853 and 480.8883 respectively, and, \({K}_{2}\) for band 10 and 11 is 1321.0789 and 1201.1442 respectively. The metadata file provided the values for \({K}_{2}\)and\({K}_{1}\). For better understanding, the thermal constant values for Landsat TM and Landsat OLI are converted from Kelvin (K) to degrees Celsius (°C) using the Eq. 0°C = 273.15K.

Step 3: Emissions from the ground surface are measured (E)

The temperature values derived above are compared to a black body. As a result, spectral emissivity (E) adjustments are required. These can be done according to the land cover type (Snyder et al. 1998) or by calculating the emissivity values for each pixel from the proportion of vegetation (Pv) data.

Where the proportion of vegetation (PV) can be calculated as:

$${P}_{V}= {\left\{\frac{\left({NDVI}_{max}-{NDVI}_{min}\right)}{\left({NDVI}_{max}-{NDVI}_{\text{m}\text{i}\text{n}}\right)}\right\}}^{2}$$

5

Step 4: Land surface temperature (LST). LST is calculated using the equation below.

$$\frac{BT}{1}+W*\left(\frac{BT}{P}\right)*Ln\left(E\right)$$

6

Where BT is the brightness temperature at the satellite, W is the wavelength of emitted radiance, P = h*c/s (1.438 10 − 2 m K), h is the Planck's constant (6.626 10–34 Js), s is the Boltzmann constant (1.38 10–23 J/K), and is the light velocity (2.998 10 8 m/s).

## Retrieval of LULC Indices

The following method was used to determine the relationship between LST and several spatial indices such as NDVI (Normalized Difference Vegetation Index), NDWI (Normalized Difference Water Index), and NDBI (Normalized Difference Built-Up Index) of the studied area (Equations7 to 9).

Calculation of Normalized Difference Vegetation Index (NDVI)

The NDVI value was extracted using the (Townshend and Justice, 1986) approach.

$$NDVI=\frac{(\text{N}\text{I}\text{R} \text{b}\text{a}\text{n}\text{d}-\text{R} \text{b}\text{a}\text{n}\text{d})}{(\text{N}\text{I}\text{R} \text{b}\text{a}\text{n}\text{d}+\text{R} \text{b}\text{a}\text{n}\text{d})}$$

7

Where NIR is the near-infrared band's DN (digital number) value and R is the red band's DN value. The NDVI value is a number that runs from 0 to 1. Low vegetation cover is indicated by values close to 0, and high vegetation density is indicated by values close to 1.

Calculation ofNormalized Difference Water Index (NDWI)

$$NDWI=\frac{(Greenband-NIRband)}{(\text{G}\text{r}\text{e}\text{e}\text{n} \text{b}\text{a}\text{n}\text{d}-\text{N}\text{I}\text{R} \text{b}\text{a}\text{n}\text{d})}$$

8

To avoid the problem of the built-up area being included in the NIR band, NDWI is used (Gu et al. 2008), ∑where green refers to the green band and NIR refers to the near-infrared band.

Calculation of Normalized Difference Built-up Index (NDBI)

The following formula (Zha et al. 2003) is used to calculate NDBI, with a value closer to 1 indicating a high density of built-up land.

$$NDBI=\frac{(\text{M}\text{I}\text{R} \text{b}\text{a}\text{n}\text{d}-\text{N}\text{I}\text{R} \text{b}\text{a}\text{n}\text{d})}{(\text{M}\text{I}\text{R} \text{b}\text{a}\text{n}\text{d}+\text{N}\text{I}\text{R} \text{b}\text{a}\text{n}\text{d})}$$

9

Where MIR is the DN from the mid-infrared band and NIR is the near-infrared band.

Pearson’s product-moment correlation coefficient (r)and Linear Regression method used ((Patra and Gavsker, 2021))

$$r=\frac{n(\sum xy)-(\sum x)(\sum y)}{\sqrt{[n\sum }{x}^{2}-\left({\sum x)}^{2}\right][n\sum {y}^{2}-(\sum {y}^{)2}]}$$

10

$${y}_{i}=\beta ˳+{\beta }^{1}x{ᵢ}^{1}+{\beta }^{2}x{ᵢ}^{2}+\dots +\beta ₚxᵢₚ+ℇ$$

11

where, for i = n observations,

yi = dependent variable,

xi = expanatory variables,

β0 = y-intercept (constant term),

βp = slope coefficients for each explanatoryvariable,

ϵ = the model’s error term (also known as the

Output maps were created after the values of several spatial indices were calculated to depict the geographic distribution of vegetation cover, water coverage, and built-up areas throughout time.