Assessing the Impact of Land Use and Land Cover Changes on Surface Temperature Dynamics Using Google Earth Engine: A Case Study of Tlemcen Municipality, Northwestern Algeria (1989–2019)

Abstract

Changes in land use and land cover (LULC) have a significant impact on urban planning and environmental dynamics, especially in regions experiencing rapid urbanization. In this context, by leveraging the Google Earth Engine (GEE), this study evaluates the effects of land use and land cover modifications on surface temperature in a semi-arid zone of northwestern Algeria between 1989 and 2019. Through the analysis of Landsat images on GEE, indices such as normalized difference vegetation index (NDVI), normalized difference built-up index (NDBI), and normalized difference latent heat index (NDLI) were extracted, and the random forest and split window algorithms were used for supervised classification and surface temperature estimation. The multi-index approach combining the Normalized Difference Tillage Index (NDTI), NDBI, and NDVI resulted in kappa coefficients ranging from 0.96 to 0.98. The spatial and temporal analysis of surface temperature revealed an increase of 4 to 6 degrees across the four classes (urban, barren land, vegetation, and forest). The Google Earth Engine approach facilitated detailed spatial and temporal analysis, aiding in understanding surface temperature evolution at various scales. This ability to conduct large-scale and long-term analysis is essential for understanding trends and impacts of land use changes at regional and global levels.

1. Introduction

From the Industrial Revolution of the mid-19th century until today, the process of urbanization has continuously accelerated. Currently, around 68% of the world’s population, equivalent to 3 billion people, reside in urban areas. This figure is predicted to increase to 5 billion by 2030 [1]. The dual processes of industrialization and urbanization are expected to gradually improve the living conditions of the population but also induce adverse environmental effects, such as air, water, and soil pollution, as well as a considerable increase in temperatures in cities, owing to the urban heat island effect (UHI) [2]. According to the specialized literature, the urban heat island phenomenon may be due to the artificialization of green surfaces and the replacement of permeable natural land with mineral materials such as concrete and asphalt because of their high capacity for storage and re-emission of solar radiation [3]. This results in a decrease in green surfaces and their associated benefits. Surface properties and composition also have a direct impact on environmental air temperature. Furthermore, Oke [2] has shown that urban geometry and the albedo of reflective building surfaces have a significant role in reflecting solar radiation and contributing to the urban heat island effect.

Concentrations of air pollutants are usually 5–25 times higher in cities than in the countryside [4]. Therefore, the temperature increase due to the urban heat island has a severe impact on human health and results in higher electrical energy use [5,6]. According to Johnson and Wilson [7], these heat islands correspond to “death islands” in high-density neighborhoods of vulnerable populations in several cities, particularly in the United States. This observation has led scientists to focus on urban characteristics and surface heat islands, because the surface temperature is the main biophysical parameter used in urban health studies [8,9] and can be quantified by remote sensing [10,11].

High-resolution satellite imagery and advanced image-processing techniques have revolutionized the monitoring and mapping of LULC. Remote sensing has become indispensable for studying landscape dynamics, which are crucial in the face of rapid urbanization and environmental changes. It provides detailed and current data on various scales, assisting researchers, urban planners, and policymakers in sustainably managing natural resources, planning cities, and conserving biodiversity. However, the precise extraction of different LULC classes remains a major challenge, particularly in semi-arid areas. Complex landscapes, seasonal variations, and spectral similarities between certain LULC classes, such as urban areas and barren land, complicate the classification process [12]. Semi-arid environments pose additional challenges owing to the predominance of bare soil and low vegetation cover, making it even more difficult to differentiate between LULC classes.

Several methods have been developed to overcome these obstacles, ranging from spectral unmixing [13] (classification) to artificial neural networks [14], as well as object-oriented and knowledge-based classification methods [15]. Among these approaches, the use of combinations of multiple spectral indices has yielded promising results.

By subtracting the NDVI from the NDBI, introduced in [16], built-up and bare areas are distinguished with positive values, while other land covers have values from 0 to −1, achieving an accuracy of 92.6% for reliable urban area mapping in Nanjing, China. Rasul et al. [17], in a study conducted in Iraq, proposed the Dry Built-up Index (DBI) and the Dry Bare Soil Index (DBSI) to map built-up and bare areas in arid and semi-arid climates using Landsat 8. The results showed an overall classification accuracy of 93% and 92%, respectively. Another study [12] identifies the most effective spectral index for differentiating urban areas from bare land in the semi-arid environments of Djelfa, Messaad, and Ain Oussera using a multi-index approach with BUI, BRBA, NDTI, and DBSI. The study shows that the combination of indices (DBSI/NDTI/BUI) offers a mapping performance with an accuracy of up to 98.7%.

In 2023, a study [18] conducted in a province of Pakistan proved the effectiveness of combining the spectral indices NDVI, MNDWI, and NDBI to classify vegetation, built-up areas, bare land, and water bodies. The performance of the models was evaluated, with RF showing the best results in terms of accuracy, surpassing a support vector machine (SVM) and Classification and Regression Trees (CARTs). The Kappa coefficients for the models were 94% for CARTs, 95% for the SVM, and 97% for RF.

Additionally, the NDTI, an index proposed in [19] to differentiate built-up areas from bare land, has also been applied to soil management and agricultural practices [20,21]. The results demonstrated that the NDTI provided a better distinction between the two classes. Furthermore, the use of the SVM algorithm for classification based on the multi-index NDTI significantly improved the accuracy and precision of mapping heterogeneous urban areas.

These methodological advancements are crucial for improving the accuracy of LULC classification, particularly in complex environments such as semi-arid areas. However, beyond land classification, it is essential to understand the environmental impacts of LULC changes, particularly on land surface temperature (LST). Faisal Mumtaz [22] asserted that urbanization influences LULC and the thermal environment. Another equally important parameter influencing LULC changes is unplanned urban expansion, which affects biodiversity, ecosystems, and surface temperature. The latter is considered an indicator for assessing the urban thermal environment [23].

Other studies have demonstrated that the spatial variation in LST influences LULC, indicating that the increase in LST is caused by changes in LULC. Several studies have examined the impact of LULC dynamics on surface urban heat islands in different cities worldwide [10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25]. The specific objectives were to evaluate the spatial variability of LST with LULC, study the relationship between different LULC indices and LST, and explore the application of remote sensing to calculate LST.

Another study [26] addressed the impact of LULC changes on the LST in the megacity of Dhaka. The objective was to quantify the impact of LULC changes on LST and map changes in the spatial and temporal distribution of LST. Landsat TIR data from 1990 to 2011 were used to retrieve LST, and transition matrices were utilized to determine LULC changes over several periods. Statistical analyses were conducted to establish the relationship between the LST and biophysical parameters. The results revealed that the expansion of urban built-up areas at the expense of natural cover (floodplains and agricultural lands) significantly affected the spatial and temporal distribution of surface temperature. The urban built-up areas consistently displayed the highest ambient radiant temperatures during the study period. A decrease in vegetation cover and an increase in urban cover were associated with an increase in LST, indicating amplification of the UHI effect over time.

One study [27] focused on the impact of LULC changes on LST and the intensity of the UHI in the city center of Sivas, Turkey, from 1989 to 2015. Landsat-4 TM, Landsat-7 ETM+, and Landsat-8 OLI satellite images were used to classify LULC into five categories based on the Maximum Likelihood algorithm: agriculture, vegetation, urban/built-up areas, water, and bare land. The classification was validated with auxiliary data and GPS reference points, achieving Kappa coefficients of 0.79 in 1989, 0.87 in 1999, and 0.92 in 2015. LST was calculated from the thermal bands of the Landsat images, and the UHI intensity was determined by comparing the urban and rural temperatures. The results show an increase in built-up urban areas and agricultural land, along with a decrease in bare land. Urban areas and bare land had the highest surface temperatures, whereas areas with dense vegetation had lower temperatures. UHI intensity increased overall by 2.25 °C during this period, indicating a rising trend in urban temperature compared to rural areas.

Another study [28] conducted in the Oued Fekan sub-basin in northwest Algeria analyzed the relationship between LST, NDVI, and NDBI using remote sensing data from Landsat-8 OLI-TIRS images for four different seasons. The results revealed that urban expansion has a significant impact on LST and that NDBI is a better indicator of urban heat island effects than NDVI. The study suggests measures, such as increasing green spaces and better urban planning, to mitigate these effects and promote sustainable development.

In an analysis [29] of the spatiotemporal dynamics of LULC changes and their effect on the creation of UHI in Tehran from 1988 to 2018, the authors used Landsat 5 TM and Landsat 8 OLI/TIRS satellite data for supervised classification based on the Maximum Likelihood algorithm to delineate five LULC classes: impervious land (IL), vegetative land (VL), water bodies (WBs), farmland (FL), and open land (OL). The results revealed an increase in impermeable areas and a decrease in vegetative land, with higher temperatures in impermeable areas and lower temperatures in vegetative land. A negative correlation was found between NDVI and LST, indicating that an increase in vegetation reduces LST, whereas a positive correlation between NDBI and LST shows that an increase in impermeable surfaces increases LST. The Kappa coefficient, reflecting classification accuracy, showed a consistent improvement from 0.841 in 1988 to 0.909 in 2018.

These studies share the common essential task of selecting and sorting appropriate satellite images to meet their specific objectives. These images must be processed and corrected to eliminate distortions and atmospheric error. Subsequently, LST and LULC maps were compiled using remote sensing and geographic information system (GIS) software. This process is lengthy and complex, increasing the risk of errors at each step from image selection to processing and data interpretation. This was also highlighted by [29] in the section “Limitations and Future Scope of the Study.”

To overcome the limitations of the traditional methods mentioned earlier and take advantage of the multi-index approach, the GEE platform has proven to be a powerful and efficient tool. With its cloud processing capability, GEE allows the analysis of large amounts of geospatial data on a global scale. Supervised classification algorithms, such as RF, available in GEE, enable the automation and refinement of LULC detection and classification with high accuracy. To fully leverage the benefits offered by the GEE platform and overcome the limitations of the traditional methods, we set the following objectives:

• Extract LST using Landsat 5 and 8 OLI (operational land imager) and TIRS (thermal infrared sensor) data from GEE and build spatial distribution maps of LST for the region of interest, Tlemcen municipality;
• Improve the accuracy of LULC classification: utilize the multi-index approach by combining the indices NDBI, NDVI, and NDTI with the GEE to achieve a Kappa coefficient greater than 0.90;
• Extract the LULC indicators NDVI, NDBI, and NDLI via the GEE platform and determine the correlations of LST with these indices.

The outcomes of this study could aid urban planners and decision makers in mitigating the adverse impacts of urban expansion on surface temperature in semi-arid regions. This study’s findings offer crucial insights into the trajectories of land use and land cover changes, their effects on LST, and the potential for implementing measures to mitigate the urban heat island effect, which could be of great importance for urban planners and municipal experts.

2. Materials and Methods

2.1. Study Area

Our study area was Tlemcen municipality, which is situated in the northwest region of Algeria, and is part of the wilaya of Tlemcen. The study area spans an approximate area of 42.1 square kilometers, ranging in latitude from 34°56′ 4.5882″ to 34°51′28.6128″ N and in longitude from −1°22′8.1006″ to −1°15′26.3154″ O (Figure 1). The average elevation of the study area is 60 m, with some points reaching up to 1000 m. The region features a semi-arid climate with a dry summer and a rainy winter, with the dry period sometimes stretching for up to 6 months. The temperature varies from 6 °C to 33 °C throughout the year, rarely dropping below 2 °C or exceeding 38 °C. On average, the annual rainfall is approximately 370 mm; February has the highest number of rainy days, with an average of 5.7 days.

The municipality of Tlemcen was selected based on several important criteria, including the fact that it is a medium-sized city with no previous studies conducted on the region, as well as the following:

-

Tlemcen is a pivotal city that connects the Sahara to northern Algeria.

-

It is particularly renowned for its significant agricultural lands, primarily consisting of wheat and similar crops.

-

There is a gap in the understanding of environmental dynamics and the potential impacts of LULC changes on LST, as no research has been conducted in this region to date.

2.2. Presentation of the Data Used

2.2.1. Data for LULC

This study utilized the GEE platform to prepare data; we acquired satellite data over the years 1989, 1999, 2009, and 2019. Choosing 10-year intervals enabled the detection and analysis of significant changes in land use, land cover, and the evolution of surface temperature over an extended period. This approach allows for a precise and targeted analysis of long-term environmental changes. Landsat 5 Landsat thematic Mapper (Landsat TM) images were used to obtain data for 1989, 1999, and 2009, whereas Landsat 8 (operational land imager/thermal infrared sensor) OLI/TRIS images were used for the year 2019. The spectral bands used in the combination of the three indices (NDBI, NDVI, NDTI) are detailed in

Table 1. Landsat band characteristics for NDVI, NDBI, NDLI, and NDTI.

Satellite	Bands	Pixel Size	Wavelength	Description
Landsat 5 TM	ρGREEN: B2	30 m	0.53–0.61 µm	Green
Landsat 5 TM	ρRed: B3	30 m	0.63–0.69 µm	Red-band reflectance
Landsat 5 TM	ρNIR: B4	30 m	0.76–0.90 µm	NIR-band reflectance
Landsat 5 TM	ρMIR: B5	30 m	1.55–1.75 µm	midinfrared
Landsat 5 TM	ρSWIR: B6	60 m	10.4–12.5 µm	Thermal Infrared
Landsat 5 TM	ρSWIR: B7	30 m	2.09–2.35 µm	Short-wave Infrared
Landsat 8 OLI	ρGREEN: B3	30 m	0.52–0.60 µm	Green
Landsat 8 OLI	ρRed: B4	30 m	0.64–0.67 µm	Red-band reflectance
Landsat 8 OLI	ρNIR: B5	30 m	0.85–0.88 µm	NIR-band reflectance
Landsat 8 OLI	ρMIR: B6	30 m	1.56–1.66 µm	midinfrared
Landsat 8 OLI	ρSWIR: B7	30 m	2.10–2.30 µm	Short-wave Infrared

.

2.2.2. Data for Retrieval of LST, NDVI, NDBI, NDLI, NDTI

During the study, four Landsat images, all from August, were used from 1989 to 2019; these images were acquired from GEE. The spatial resolution of all these images was 30 m (

Table 2. Details of Landsat datasets (1989–2019).

Satellite Type/Sensor	Acquisition Date	Time (GMT)	Path/ Row	Cloud Cover	Thermal	Conversion Constants
					K1	K2
Landsat 5 LANDSAT/LT05/C01/T1_SR	11 August 1989	10:05:25	198/36	6	607.76	1260.56
Landsat 5 LANDSAT/LT05/C01/T1_SR	7 August 1999	10:16:19	198/36	0	607.76	1260.56
Landsat 5 LANDSAT/LT05/C01/T1_SR	2 August 2009	10:27:45	198/36	0	607.76	1260.56
Landsat 8 OLI LANDSAT/LC08/C01/T1_SR	5 August 2019	10:44:29	199/35	6.81	774.855	1321.078

). The images were filtered using the desired date and a cloud cover of less than 10%. We used the split window algorithm, which is strongly recommended by the scientific literature [30]. Furthermore, all images used for surface temperature extraction included previously atmospherically corrected SR from the Landsat 8 OLI/TIRS sensors.

The images contained five visible and near-infrared (NIR) bands and two short-wave infrared bands, processed to orthorectified surface reflectance (SR), and two thermal infrared bands processed to orthorectified brightness temperature (BT). These data were atmospherically corrected using Land Surface Reflectance Code (LaSRC) [31] and include a cloud, shadow, water, and snow mask produced using the C version of Function of Mask (CFMASK) and a per-pixel saturation mask (https://gee.stac.cloud/B69zEhWZA8UDdv95Jq9wpQUFv9mXwakFcBX7TsndhevENXNKqRVMheoKCQ2tWvJGp1N5eXUrsaU1KErZLfDRiiy8a3tEeDrk4baT3TYq/2FFvzA2zeqoVZ5Nv2es9fLhN3K3iBAvr84aKpTytQ3v5rrR6KJ1qJdc1Y7nmnFmdVwP5KsDC45XctKuHj1hPSkP6NhTeRU3zECHa4iVpd1xdXua4HWnF5hLF9s2eR), accessed on 1 June 2024. For Landsat 5, SRs were derived using the Landsat Ecosystem Disturbance Adaptive Processing System algorithm, which calculates the radiative transfer for atmospheric data from the moderate resolution imaging spectroradiometer (MODIS) and National Centers for Environmental Prediction (NCEP).

The following section describes the step-by-step process of mapping LULC, generating LST, NDBI, NDVI, and NDLI and quantifying the correlation between LST and indices. Figure 2 summarizes the conceptual research framework of this study.

2.3. LULC Classification

To extract LULC information, we used the random forest model with supervised classification, as it is considered the most effective algorithm for this purpose [30]. Subsequently, the land use of the study area was classified into four distinct categories, namely urban, vegetation, bare land, and forest. The urban category encompasses residential, commercial, industrial, and road areas, as well as all impermeable surfaces. The vegetation category primarily includes agricultural lands, wheat fields, or similar crop- and other grass-covered areas, whereas bare land refers to exposed soil and surfaces lacking vegetation or construction.

This section explains the code script computing in GEE for extracting the LULC of Tlemcen municipality. Figure 2, step 1, summarizes the LULC extraction in the GEE.

To select the date range for Landsat, we defined the start and end dates. Next, we loaded Landsat images and filtered them using the specified geographic zone and date range. Additionally, we applied a filter to select images with less than 10% cloud cover, which was necessary to reduce the impact of cloud cover and seasonal variations on the classification results.

We printed a collection of filtered images to the console. Then, we calculated the median of the filtered image collection and clipped it to the specified zone. To display the image on the map, we added it as a layer with specific visualization parameters (min and max values, and the bands used for RGB composition). We centered the map view on the specified zone with a zoom level of 9 and selected specific bands (Table 1) from the image for further processing. We then used a multi-index approach by combining the NDVI, NDBI, and NDTI to provide essential information on vegetation, agricultural lands, and built-up areas.

We created training data by merging different land cover points (urban, bare soil, vegetation, and forest (

Table 3. Description of LULC classes in the case study: Tlemcen.

No.	LULC Class	Description of Class
1	Urban	Residential, commercial, industrial, roads, and all impervious surfaces
2	Barren land	Exposed soil and surfaces without any vegetation or construction
3	Vegetation	Agricultural land, wheat fields, or similar crop- and other grass-covered areas
4	Forest	Area with forest features

)) into a single collection called training points. Next, we created training regions by ‘overlaying’ these points on the image stack, selecting all the bands and associating them with the land cover property. The scale was set at 30 m.

We then randomized the collection to ensure a good mix of training and validation data and split the dataset into approximately 70% for training and 30% for validation.

Next, we classified the data using a random forest model with supervised classification. This determines the number of trees from a collection of the sample training data. Previous research has shown that RF algorithms effectively handle data overfitting and possess a high processing power [30]. Furthermore, RF methods tend to offer higher accuracy than other classification approaches, such as support vector machines, Maximum Likelihood, and single decision trees [30]. This approach involves training the model with labeled data to ensure high classification accuracy. The next phase is necessary to validate the previously performed classification.

Accuracy Assessment of Supervised Classification

We used an accuracy assessment process to verify the accuracy of computer-classified maps and produced descriptive statistics to compare the classification results with actual field data [32,33]. The GEE platform enabled us to obtain the medians of bands 2, 3, 4, 5, 6, and 7 for each scene throughout the study year. Approximately 450 samples were generated for each year to ensure the dependability of the results. To evaluate the performance of the random forest model, the samples were divided into training and validation sets, with 70% of the samples being employed as training samples and the remaining 30% used for validation.

The accuracy of the retained classified imagery was assessed using 150 reference pixels for each class identified by a sample selected from the displayed Landsat imagery. For each selected year, about 600 samples were created to ensure the reliability of the results. We then analyzed the categorized image accuracy using the validation samples from 4 years.

The confusion matrix is frequently used to obtain analytical and descriptive data for classification accuracy. It consists of values shown in columns and rows, indicating the contrast between different sample points (polygons, pixels, or pixel clusters) assigned to a particular land cover class and the ground conditions of the class as a whole [34]. The matrix’s assessment indices for overall accuracy are producer and user accuracy and the Kappa coefficient (KC) [35]. The producer accuracy is the ratio of all classified pixels in the error matrix diagonals to all classified pixels in that category in the error matrix column.

The KC, based on [29], was calculated using Equation (1):

$$KappaCoefficient\left(KC\right)=\frac{N{\sum }_{i=1}^{r}xii-{\sum }_{i=1}^r\left(x_{i+}\times x_{+i}\right)}{N^2-{\sum }_{i=1}^r\left(x_{i+}\times x_{+i}\right)}$$

(1)

where N is the sum of pixels in the error matrix; r is the sum of columns/rows; xii is the number of correctly classified pixels in the ith column and row; x + i is the sum of pixels in the i the column; xi+ is the sum of pixels in the i the row; and N² is the square of the total number of pixels.

2.4. LST Retrieval

This second step involves calculating the LST on the GEE platform. The main factors considered in the LST calculation were NDVI, brightness temperature (BT), and land surface emissivity (LSE). The thermal bands from Landsat 5 and Landsat 8 were also considered in the LST calculation process. Figure 2, step 2, summarizes the LST retrieval process.

First, the spectral radiance of the thermal band was converted to BT using Equation (2) in the tool’s algorithm [36,37]:

$$T_B=\frac{k_2}{\mathit{ln}\left(\frac{K_1}{L_{\lambda }}+1\right)}$$

(2)

where $T_B$ is the at-satellite BT (C°), $L_{\lambda }$ is the top-of-atmosphere (TOA) spectral radiance (W/(m² srad µm)), $K_1$ is the band-specific thermal conversion constant from the metadata (W/(m² srad µm)), $k_2$ is the band-specific thermal conversion constant from the metadata (K), and ln is the natural logarithm. We then derived the emissivity values of the corrected LST (C°) with the aid of the at-satellite $T_B$ [38,39] using Equation (3):

$$LST\left(K°\right)=\frac{T_B}{[1+\lambda \left(\frac{T_B}{E}\right)Ln\left(\epsilon \right)]}$$

(3)

where λ is the wavelength of the emitted radiance (i.e., 11.5 µm in band 6 for Landsat 4/5/7 and 10.8 µm in band 10 for Landsat 8), E is (h × v)/s (1.4388 × 10⁻² mK), h is Planck’s constant (6.626 × 10⁻³⁴ mK), v is the velocity of light (2.998 × 108 m/s), s is the Boltzmann constant (1.38 × 10⁻²³ JK), and ε is the emissivity of the land surface. Then, the emissivity of the land surface (ε) was calculated using Equation (4) based on [33,40], as follows:

$$\left(\epsilon \right)=N\left(Pv\right)+\mathrm{n}$$

(4)

where N is 0.004; n is 0.986; and Pv is the vegetation proportion expressed by Equation (5) as described in [40,41]:

$$Pv={\left(\frac{NDVI-NDVImin}{NDVImax-NDVImin}\right)}^2$$

(5)

where NDVI is the value of the Digital Number (DN) obtained from the NDVI image, and NDVImax and NDVImin are the highest and lowest DN values obtained from the NDVI image. Last, we converted the LST value (in °K) into degrees Celsius (°C) using Equation (6) as described in [41]:

$$LST\left(°\mathrm{C}\right)=LST\left(°\mathrm{K}\right)-273.15$$

(6)

2.5. Calculation of Spectral Indices

The NDVI is one of the most commonly used indices in remote sensing [42]. This index uses the ratio of red to NIR bands (Table 1) to identify live green vegetation and takes values ranging between −1 and +1, where negative values indicate non-vegetated areas and positive values represent areas covered in vegetation [43].

We used the NDVI to investigate the vegetation distribution and extract emissivity values within the study area for the four years of interest: 1989, 1999, 2009, and 2019. The NDVI was calculated using Equation (7):

$$NDVI=\frac{\left(\rho NIR-\rho Red\right)}{\left(\rho NIR+\rho Red\right)}$$

(7)

where ρNIR is the NIR-band reflectance and ρRed is the red-band reflectance.

The NDBI is another indicator used to understand urban climates [44]. This index uses the ratio of the mid-infrared (MIR) band to the near-infrared (NIR) band (Table 1) to identify built-up areas and impervious surfaces [45]. Its values range from −1 to +1, where a significant positive value denotes a built-up area, a small positive value indicates bare soils, and a negative value represents water bodies and vegetation [43]. This index allowed us to map and analyze land cover and obtain spatial information for built-up and impervious areas for 1989, 1999, 2009, and 2019. The NDBI was calculated using Equation (8):

$$NDBI=\frac{\left(\rho MIR-\rho NIR\right)}{\left(\rho MIR+\rho NIR\right)}$$

(8)

where ρMIR and ρNIR represent the reflectance values for the mid-infrared and near-infrared bands of Landsat images, respectively.

The NDLI has been shown to be a reliable indicator for estimating potential surface evapotranspiration, incorporating the three regularly utilized satellite mission channels of green, red, and shortwave-infrared (Table 1). The NDLI has the unique capacity to optimize spectrum sensitivity on biophysical land surface parameters [46]. This index was calculated using Equation (9):

$$NDLI=\frac{\left(\rho GREEN-\rho RED\right)}{\left(\rho GREEN+\rho RED+\rho MIR\right)}$$

(9)

where ρGREEN, ρRED, and ρMIR represent the reflectance values for the green, red, and mid- infrared - bands of Landsat images, respectively.

2.6. Correlation Analysis

We assessed the correlations of LST with the NDVI, NDBI, and NDLI of the study area (reflecting the surface and LST changes in the Tlemcen municipality) for each year using the correlation coefficient (R). Through the collection of independent pixels selected at random from the entire study region, these indices were used to develop quantitative connections with LST.

$$r=\frac{\sum _{\mathrm{i}=1}^{\mathrm{n}}\left({\mathrm{x}}_{\mathrm{i}}-\overline{\mathrm{x}}\right)\times {(\mathrm{y}}_{\mathrm{i}}-\overline{\mathrm{y}})}{\sqrt{{\sum _{\mathrm{i}=1}^{\mathrm{n}}\left({\mathrm{x}}_{\mathrm{i}}-\overline{\mathrm{x}}\right)}^2}\times \sqrt{\sum _{\mathrm{i}=1}^{\mathrm{n}}{({\mathrm{y}}_{\mathrm{i}}-\overline{y})}^2}}$$

(10)

where r represents Pearson’s correlation coefficient; x represents the independent variables measuring the value of x_i; y represents the dependent variable measuring value of y_i; x_i and y_i represent the individual sample points indexed (i); and x and y represent the mean of the sample.

3. Results

This section is divided into four parts: the first part presents the results of the accuracy assessment of LULC classification, the second part covers the results of LULC, the third part presents the results of LST, and the fourth part discusses the relationship between LST, LULC, and indices.

3.1. Accuracy Assessment of LULC Classification

The Kappa coefficient was employed to assess the LULC classification’s correctness in more detail. The total accuracy throughout the four periods was above 90%, suggesting a robust land cover classification and good agreement between the referenced and classified maps [47].

Previous studies [48] indicate that a Kappa coefficient greater than 75 means that the classification and reference data are sufficient to assess accuracy. If the accuracy of the classification is within acceptable limits, the requirements are met.

Table 4. Producer accuracy and Kappa coefficient for LULC.

Years	Kappa Coefficient	Producer Accuracy
2019	0.98	0.98
2009	0.97	0.98
1999	0.96	0.97
1989	0.96	0.96

presents the detailed accuracy of LULC.

3.2. LULC of the Study Area

Tlemcen’s classified land cover maps for the periods of interest (i.e., 1989, 1999, 2009, and 2019) are presented in Figure 3 and quantified in

Table 5. LULC distribution in 1989, 1999, 2009, and 2019.

Years	Urban (ha)	%	Barren Land (ha)	%	Vegetation (ha)	%	Forest (ha)	%	Total
1989	1079	25.65	1149	27.3	1821	43.27	159	3.78	100%
1999	1343	31.91	1382	32.84	1315	31.25	168	4.00	100%
2009	1349	32.06	1416	33.66	1266	30.08	177	4.20	100%
2019	1786	42.44	499	11.85	1726	41.01	197	4.70	100%

and Figure 4. The LULC was classified into four broad classes, comprising urban areas, vegetation, barren land, and forest. The findings of the LULC classification study for 1989 indicated that vegetation was widespread, particularly in the central and southeastern regions of the study area. Urban zones were mostly concentrated in the central and western parts of the municipality (old town and residential neighborhoods). On the other hand, the northern region of the municipality consisted of bare and rocky soils, while the southern region was characterized by dense forest cover. The results of the LULC classification study conducted in 2019 revealed a significant evolution in the different land cover classes, with a notable change in their spatial distribution. Urban expansion was observed towards the north, mainly due to geographical constraints in the south of Tlemcen. This urban growth resulted in a reduction in bare soil in the north and a decrease in vegetation, primarily consisting of agricultural land, in the east.

The classification revealed that built-up areas increased by 264 ha (3%) during the first decade (1989–1999) and 440 ha (10.39%) during the last decade (2009–2019). Conversely, built-up areas increased slightly (by 6 ha; 0.14%) between 1999 and 2009, whereas the increase was considerable between 2009 and 2019. The amount of vegetation cover decreased by 12.02% (506 ha) during the first 10 years (1989–1998). However, it increased by 460 ha (10.93%) between 2009 and 2019. Overall, the results indicated a decrease of 95 ha (−2.26%) in vegetation cover. Barren land surface increased by 267 ha (6.35%) during the first three decades (1989–2009). Conversely, it decreased sharply (by 917 ha; 21.8%) in the last decade. This decrease was due to use of barren land to construct new extensions totaling 551 ha. The barren land class was the most dominant between 1999 and 2009. However, it considerably decreased in 2019 (21.84%).

Finally, the results showed a significant shift in urban classes, contributing significantly to the formation of urban heat islands and other environmental issues. These results are consistent with previous findings [29,49,50].

3.3. LST of the Study Area

The LST was analyzed for August, which is considered the hottest month in our study area, over the 4 years of study. Figure 5 depicts the spatial distribution of surface temperatures for 1989, 1999, 2009, and 2019.

Table 6. Summary statistics of LST for Tlemcen during the decades of interest.

Year	Minimum	Maximum	Mean	Standard Deviation
1989	24.95	37.26	32.23	2.07
1999	24.16	38.06	31.31	2.42
2009	27.95	44.46	37.40	2.98
2019	27.30	41.03	35.12	2.42

summarizes the statistical temperature results. The surface temperatures in Tlemcen municipality were in the ranges of 27.30–41.03 °C, 27.95–44.46 °C, 24.15–38.06°, and 24.95–34.26 °C during the four distinct periods (i.e., 2019, 2009, 1999, and 1989, respectively).

Furthermore, the LST analysis revealed that the minimum and maximum surface temperature values in 1989 were 24.95 °C and 37.26 °C, respectively. In 1999, the minimum was 24.16 °C, and the maximum reached 38.06 °C. In 2009, the minimum and maximum both increased, reaching values of 27.96 °C and 44.46 °C, respectively. Finally, in 2019, the minimum surface temperature was 27.30 °C, and the maximum was 41.03 °C. The spatialization results for the surface temperature revealed an increase of 2.89 °C between 1989 and 2019. Consequently, the year 2009 is marked by the highest average temperature.

The industrial area located in the northeast of the town had the warmest area over the 4 years of this research; this was owed to the materials covering the roofs of the buildings (asbestos and reinforced concrete). The two northeastern and northwestern parts also had high-temperature values since 1999, as these areas have clusters of new dwellings (collective housing, dense individual housing, and informal housing). Notably, the spatial patterns of surface temperature showed that the urban core had cooler temperatures compared with the neighboring districts. This distinction arises from the lush environment of the core, which has a significant abundance of trees.

3.4. Relationship between LULC and LST

For a thorough analysis of the relationship between LST and LULC, it is essential to study the thermal patterns of the different LULC classes [51].

Table 7. Mean LST (°C) over each LULC class in Tlemcen from 1989 to 2019.

Year	Urban	Barren Land	Vegetation	Forest
1989	30.28	31.064	30.54	26.67
1999	30.7	33.46	30.41	26.4
2009	36.04	39.73	36.69	31.36
2019	35.01	37.67	35.27	30.02

displays the mean surface temperature values by LULC class (urban, barren land, vegetation, and forest) during 1989, 1999, 2009, and 2019. The barren land class recorded the highest temperatures of 34.21 °C in 1989 and 33.34 °C in 1999, increasing to 39.54 °C in 2009 and 37.27 °C in 2019, an average increase of 4.62°. The average temperatures of the urban class were 32.02 °C, 30.87 °C, 36.72 °C, and 35.05 °C in 1989, 1999, 2009, and 2019, respectively. The LST of the urban class increased by 5 °C between the first two and the last two decades. Ranking third, the vegetation class also experienced an increase in temperature between 1989 and 2009, with an average increase of 5.03 °C (Table 7). Thus, the barren land and urban areas had the highest temperatures, compared with the vegetation and forest classes. This difference can be explained by the high evaporation rates of vegetation, reducing the temperatures at the soil’s surface [52,53].

Although the average temperature values for each class increased, 2009 was the hottest year (Table 7). The spatial distribution of temperatures showed elevated temperatures in the northeast and west of Tlemcen. These areas include the new urban extensions. We also found that temperatures were highest in the northeast industrial area during the four years studied (i.e., 1989, 1999, 2009, and 2019). According to the remote sensing data, Tlemcen experienced intense urbanization during the last decade of this study. In fact, the northern section of the municipality experienced soil artificialization. In addition, the center of the municipality experienced a transition from agricultural surfaces to mineral surfaces. However, this change was not consequence-free, as it harmed the natural environment and contributed considerably to the increase in surface temperatures of the region.

3.5. Mapping of NDVI and Relationship with LST

NDVI maps of Tlemcen are presented in Figure 6 for 1989, 1999, 2009, and 2019. Statistical information for NDVI is shown in

Table 8. Summary statistics of NDVI for Tlemcen during the decades of this study.

Year	Minimum	Maximum	Mean	Standard Deviation
2019	0.038	0.73	0.26	0.13
2009	0.054	0.68	0.29	0.10
1999	0.042	0.43	0.18	0.07
1989	0.06	0.40	0.18	0.05

. The mean NDVI observed in 1989 and 1999 was 0.18. However, after this first decade, the value of NDVI increased in 2009 to about 0.29, then decreased in 2019 to 0.26. Plausible reasons for the increase in the average value of NDVI are the building of new parks in urban areas, such as Park El Hartoun and Park Lalla Setti in the south, and the increase in forest area, as shown in Figure 4. All these activities were initiated to organize the “Tlemcen Islamic Culture Capital” event in 2011 (chosen by ISESCO: the Islamic World Educational, Scientific, and Cultural Organization). The highest NDVI values were observed between 1989 and 2009, in the south (Les Petits Perdreaux Forest) and the center (presence of agricultural land) of the region. However, the northeast and northwest of the municipality had low NDVI values, which can be explained by the prevailing rocky soil and the presence of abundant impervious surfaces (Figure 6).

The correlation between NDVI and LST was negative during the four years of interest, as shown in

Table 9. Correlations between NDVI and LST in 1989, 1999, 2009, and 2019.

	1989	1999	2009	2019
Correlation	−0.131 1	−0.472 1	−0.382 1	−0.408 1

¹ The correlation is meaningful at the level of 0.01 (bilateral).

. This correlation indicates that LST increases when NDVI decreases, and decreases when the NDVI increases. Les Petits Perdreaux Forest, located in the south of the municipality, recorded the highest NDVI and the lowest LST values. The outcomes of this study are in agreement with those of previous investigations [10,29,54,55].

3.6. Mapping of NDBI and Its Relationship with LST

NDBI maps of Tlemcen are presented in Figure 7 for 1989, 1999, 2009, and 2019. Statistical information for NDBI is provided in

Table 10. Summary statistics for NDBI in Tlemcen during the decades of interest.

Year	Minimum	Maximum	Mean	Standard Deviation
1989	−0.13	0.20	0.065	0.07
1999	−0.17	0.19	0.060	0.06
2009	−0.36	0.66	0.056	0.09
2019	−0.31	0.20	0.033	0.09

. The mean NDBI values observed in 1989 and 1999 were 0.065 and 0.060, respectively. Following the first decade, the value of NDBI decreased in 2009, whereas it rose to values of 0.056 and 0.033 in 2019. This was clearly the result of a decrease in barren land surface areas.

As shown in Figure 7, the highest NDBI values in 1989, 1999, and 2009 were recorded in the northern area of the region. The latter is characterized by barren land that underwent unplanned urbanization in the northwest during the 1990s, whereas rapid urbanization took place (in terms of collective and individual dwellings) in the northeast area, coupled with the absence of vegetation. The average NDBI values were spread over Tlemcen municipality, indicating significant impermeability of surfaces along with vegetation surfaces (Figure 7a).

There were positive correlations between NDBI and LST during the last four years, as shown in

Table 11. Correlations between NDBI and LST in 1989, 1999, 2009, and 2019.

	1989	1999	2009	2019
Correlation	0.319 1	0.677 1	0.626 1	0.659 1

¹ The correlation is meaningful at the level of 0.01 (bilateral).

, with the R value increasing between 1989 and 2019. This phenomenon is represented by a maximal value of NDBI covering the whole municipality. These results are in accordance with those of previous studies [52,56].

3.7. Mapping of NDLI and Its Relationship with LST

NDLI maps of Tlemcen are presented in Figure 8 for the four years of interest (i.e., 1989, 1999, 2009, and 2019). Statistical information for NDLI is provided in

Table 12. Summary statistics for NDLI in Tlemcen during the decades of interest.

Year	Minimum	Maximum	Mean	Standard Deviation
1989	−0.10	0.003	−0.062	0.015
1999	−0.11	0.0022	−0.058	0.017
2009	−0.13	0.013	−0.052	0.020
2019	−0.12	0.02	−0.062	0.020

. The mean NDLI observed in 1989 was −0.062, increasing to −0.058 and −0.052 in 1999 and 2009, respectively.

As shown in Figure 8, the maximum values of NDLI in 1989 show a peak, mainly in the Les Petits Perdreaux Forest surface, located on the road to the villages of Safsaf and Sidi Othmane. The latter is known for its fertile agricultural land. Over the years considered, the maximal values dissipated. However, these values were centered at Les Petits Perdreaux Forest owing to the loss of agricultural surface and mineralization of the vegetation surfaces in the center of the municipality. Furthermore, there was a gap between the NDLI values, based on the shift of the standard deviation from 0.015 in 1989 to 0.020 in 2019 (Table 12).

There were negative correlations between NDLI and LST, with r = −0.742, −0.647, and −0.579 for 1999, 2009, and 2019, respectively (

Table 13. Correlations between NDLI and LST in 1989, 1999, 2009, and 2019.

	1989	1999	2009	2019
Correlation	−0.260 1	−0.742 1	−0.647 1	−0.579 1

¹ The correlation is meaningful at the level of 0.01 (bilateral).

), which can be attributed to high-water-content surfaces, as the amount of water transferred to the atmosphere through evaporation reduces surface temperature. Consequently, the results align with those of previous studies [46,57].

4. Discussion

4.1. Urban Transformation and Thermal Changes in Tlemcen: 30 Years of Evolution

Over the past three decades, the evolution of urban areas has been marked by irregular fluctuations. The increase was very slight between 1999 and 2009 due to urban policies implemented during that period, when the country was emerging from a period of widespread insecurity, which slowed down the construction of major housing projects. During the last decade, from 2009 to 2019, the city experienced a boom in the construction of major housing projects, such as the National Agency for Housing Improvement and Development (AADL) and participatory social housing (LSP), primarily in the northeastern part of the study area. These results align with those obtained in [58] for the intercommunal grouping of Guelma, where urbanization saw a slight increase between 2000 and 2010.

The study of urban expansion, particularly in the northeast and northwest of Tlemcen, provides crucial information on temperature fluctuations in the city. Notably, neighborhoods built after 1990 exhibited higher temperatures than those built in the old city. Conversely, the historic center within the city walls displays significantly cooler temperatures. This observation may be attributed to the age of this area of the city, which encompasses the confluence of the colonial fabric, constructed according to colonial principles, and the medieval urban fabric, constructed with specific materials characterized by high thermal inertia, developed according to Arab-Muslim urbanization principles. This situation underscores the importance of preserving the unique characteristics of Tlemcen’s historic center, such as the plane trees and horse chestnuts, as well as the traditional architecture, which contributes to urban sustainability and the reduction in the urban heat island effect.

This study revealed a significant increase in temperature across all categories, with the barren land class being the hottest each year. This situation can be explained by its low thermal inertia and low evapotranspiration due to the absence of vegetation cover. These results are consistent with the conclusions of [59], and the data from the correlation matrix presented in the

Table 14. Correlation matrix between the indices in 2019.

2019	Lst	Ndvi	Ndbi	Ndli
Lst	1	−0.401	0.659 1	−0.59 1
Ndvi	−0.40 1	1	−0.74 1	0.73 1
Ndbi	0.65 1	−0.74 1	1	−0.80 1
Ndli	−0.59 1	0.73 1	−0.80	1

¹ The correlation is meaningful at the level of 0.01 (bilateral).

,

Table 15. Correlation matrix between the indices in 2009.

2009	Lst	Ndvi	Ndbi	Ndli
Lst	1	−0.38 1	0.62 1	−0.64 1
Ndvi	−0.38 1	1	−0.66 1	0.64 1
Ndbi	0.62 1	−0.66 1	1	−0.66 1
Ndli	−0.64 1	0.64 1	−0.66	1

¹ The correlation is meaningful at the level of 0.01 (bilateral).

,

Table 16. Correlation matrix between the indices in 1999.

1999	Lst	Ndvi	Ndbi	Ndli
Lst	1	−0.47 1	0.67 1	−0.74 1
Ndvi	−0.47 1	1	−0.76 1	0.69 1
Ndbi	0.67 1	−0.76 1	1	−0.82 1
Ndli	0.67 1	0.69 1	−0.82 1	1

¹ The correlation is meaningful at the level of 0.01 (bilateral).

and

Table 17. Correlation matrix between the indices in 1989.

1989	Lst	Ndvi	Ndbi	Ndli
Lst	1	−0.13 1	0.31 1	−0.26 1
Ndvi	−0.13 1	1	−0.76 1	0.71 1
Ndbi	0.31 1	−0.76 1	1	−0.87 1
Ndli	−0.26 1	0.71 1	−0.87 1	1

¹ The correlation is meaningful at the level of 0.01 (bilateral).

show a strong negative correlation between NDBI and NDLI over the four years studied.

The results of this study showed that the average temperatures of urban and vegetated areas were almost similar. However, the correlations of these two types of areas with surface temperature were completely opposite. This divergence can be explained by the fact that the vegetated areas identified in the images were mainly agricultural lands, particularly dedicated to wheat cultivation and other similar crops, characterized by very sparse vegetation. In August, when the surface temperatures were measured, these areas had a cover of dehydrated wheat stalks after harvest or were devoid of cover for other crops. According to [27], these areas are considered bare lands, and the temperatures of this type of class vary depending on the soil type (clayey, sandy, peaty) and exogenous parameters, such as air temperature. Indeed, Mallik, J. et al. [60] demonstrated that the average air temperature corresponds to that measured using remote sensing on Landsat images. These results corroborate the conclusions of [27,29,58].

In accordance with the results obtained by other researchers [27,33,61], the correlation between vegetation and surface temperature was negative, which can be attributed to the decrease in vegetation cover and its transformation into mineral surfaces. However, it is important to note the variability in the standard deviation, which reached its maximum in 2019 with a value of 0.13, compared with a value of 0.05 in 1989. This reveals an increased heterogeneity in vegetation cover in 2019 compared to 1989, when the vegetation cover was more homogeneous, leading to a weaker correlation between vegetation and surface temperature.

The mapping of evapotranspiration over the three-decade period covered in this study presents a clear visualization of the thermal dynamics at play. However, an intriguing aspect emerges, with a significant decrease in evapotranspiration at the historical center. This reduction can be justified by the progressive degradation of the landscape framework including the construction of new buildings with low-albedo and low-inertia materials, at this specific location in the city over time.

These results underscore the impact of urban development on local thermal conditions, thus emphasizing the need for integrated and sustainable urban management practices. It is essential to adopt sustainable approaches in urban planning to minimize adverse effects on the microclimate of the city. This involves promoting environmentally friendly architecture, developing green spaces, preserving historical areas, and implementing urban policies that aim to balance growth with the conservation of natural resources. Sustainable urban management is crucial for preserving the quality of life of residents while maintaining the unique and historical character of the city of Tlemcen.

4.2. GEE: Advancing LULC and LST Research

Conventional approaches to evaluating the effects of alterations in land usage on surface temperature usually entail the manual compilation of data from diverse sources, including meteorological stations, field assessments, and satellite information. Such methods can be laborious, requiring substantial investment of time and effort. Researchers are required to carry out on-site investigations, establish meteorological stations, and accumulate data over designated periods to determine the impact of land use changes on surface temperature. However, the development of the GEE platform has had great significance in this field, yielding notable advantages in the exploration of alterations in land use and their effects on surface temperature. Various pivotal factors warrant careful consideration within this comparative evaluation, which is summarized in

Table 18. Contrasting features of Google Earth Engine (GEE) and traditional approaches in geospatial data analysis.

Feature	GEE	Conventional Method
Data scale	Access to global data with high spatial and temporal resolution	Limited to local or regional scales
Processing time	Rapid processing of large amounts of data using cloud computing	Requires manual pre-processing of data before analysis
Data selection	Offers a variety of pre-processed data, reducing manual data collection	Requires effort to select, filter, and pre-process data
Advanced analysis	Provides advanced analysis tools such as spectral imaging, time series analysis, and image classification	Traditional analysis approach, less equipped with advanced tools
Accessibility	Accessible online to researchers worldwide	May be limited by availability of local resources

below.

5. Conclusions

This study evaluated the impact of LULC changes on LST in Tlemcen municipality, from 1989 to 2019, leveraging the power of GEE. Using a multi-index approach by combining the three spectral indices NDVI, NDTI, and NDBI, we achieved results with Kappa coefficients greater than 0.95, thus meeting our set objective for Kappa coefficients. The analysis included four LULC classes, urban, vegetation, barren land, and forest, showing significant urban growth and a corresponding decrease in barren and agricultural lands. Surface temperatures increased by an average of 2.89 °C, with an increase of 5 °C in urban and vegetative areas, 3.35 °C in forests, and 6.61 °C in barren lands. Negative correlations were found between surface temperature (LST) and both NDVI and NDLI, indicating lower temperatures with denser vegetation. Positive correlations with NDBI were associated with higher temperatures due to increased impermeable surfaces. This study highlights the effectiveness of GEE in analyzing the impacts of LULC changes and suggests its use for better urban planning and land management. For future studies, it will be essential to combine other spectral indices to more precisely distinguish between different LULC classes. The use of other satellites, such as Sentinel, could improve the quality and resolution of the data obtained. Additionally, varying the studies across different seasons would be beneficial to capture the seasonal variations of LULC classes and better understand their dynamics.

Finally, remote sensing techniques and the new GEE platform have dramatically changed how we access and apply geographic data. The extensive geodata library, advanced algorithms, and cloud computing capabilities of GEE, including the classification model, allow for easier, faster, and more accessible data analysis, allowing us to better study the impact of LULC changes on the temperature in a semi-arid climate. The results presented here could serve as a theoretical basis for urban planning and decision-making for better future land management to prevent excessive and anarchic land use.