Revealing Land-Use Dynamics on Thermal Environment of Riverine Cities Under Climate Variability Using Remote Sensing and Geospatial Techniques

Abstract

Urbanized riverine cities in southern Asian developing countries face significant challenges in understanding the spatiotemporal thermal impacts of land use/land cover (LULC) changes driven by rapid urbanization and climatic variability. While previous studies have investigated factors influencing land surface temperature (LST) variations, gaps persist in integrating Landsat imagery (7 and 8), meteorological data, and Geographic Information System (GIS) tools to evaluate the thermal effects of specific LULC types, including cooling and warming transitions, and their influence on air temperature under variable precipitation patterns. This study investigates LST variations in Islamabad, Pakistan, from 2000 to 2020 using quantile classification at three intervals (2000, 2010, 2020). The thermal contributions of each LULC type across the LST-based temperature classes were analyzed using the Land Contribution Index (LCI). Finally, Warming and Cooling Transition (WCT) maps were generated by intersecting LST classes with 2000 as the baseline. Results indicated a rise in LST from 32.39 °C in 2000 to 45.63 °C in 2020. The negative LCI values revealed that vegetation and water bodies in lower temperature zones (Ltc_1 to Ltc_3) contributed to cooling effects, while positive LCI values in built-up and bare land areas in higher temperature zones (Ltc_5–Ltc_7) exhibited warming effects. The WCT map showed a general warming trend (cold-to-hot type) from 2000 to 2020, particularly in newly urbanized areas due to a 49.63% population increase, while cooling effects (hot-to-cold type) emerged in the newly developed agricultural lands with a 46.46% rise in vegetation. The mean annual air temperature gap with LST narrowed from 11.55 °C in 2000 to 2.28 °C in 2020, reflecting increased precipitation due to increasing yearly rainfall from 982.88 mm in 2000 to 1365.47 mm in 2020. This change also coincided with an expansion of water bodies from 2.82 km² in 2000 to 6.35 km² in 2020, impacting the local climate and hydrology. These findings highlight the importance of green spaces and water management to mitigate urban heat and improve ecological health.

1. Introduction

Changing climatic conditions, compounded by anthropogenic activities, have disrupted the energy balance between urban landscapes and local climates in riverine regions globally [1,2,3,4]. Urbanization has accelerated, with 55% of the population now living in urban areas, and this percentage is projected to rise to 68% by 2050 [5]. Unplanned urban expansion has transformed natural landscapes, i.e., soil, water, and vegetation, into impervious surfaces, industrial infrastructure, and high-rise buildings [6,7], altering surface energy balance, increasing land surface temperatures (LST), and degrading air quality [8,9].

The urban heat island (UHI) effect, linked to land use/land cover (LULC) changes and LST variations, is particularly pronounced in developing regions [10,11,12,13,14,15]. Urban populations are increasingly vulnerable to extreme climate events linked to UHIs, including heatwaves, droughts, storms, and heavy rainfall [10,12,16,17,18]. More critically, rising temperatures exacerbate energy demands for electricity, water, and cooling systems [19] while also heightening health risks, such as respiratory and cardiovascular diseases [20]. Global surface temperatures are also expected to rise by 0.3–0.7 °C by 2035, depending on greenhouse gas emission trends [21], further straining riverine ecosystems [9,22]. Numerous studies have been conducted to understand the relationship between LULC and UHIs in several areas around the globe, such as in China [15,23,24], India [14], and Pakistan [4,25,26].

Monitoring LST is crucial for evaluating UHIs intensity, especially given the significant regional variability in rapidly growing cities [4]. Traditional methods relying solely on in situ meteorological data fail to capture the complex interplay between LULC dynamics and LST. Remote Sensing (RS), integrated with Geographical Information System (GIS) frameworks, offers a robust alternative to LST analysis [11,14,27,28]. Advances in multispectral sensors, including Landsat Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+), Operational Land Imager (OLI), Sentinel-2, and ASTER, have provided essential datasets for assessing urban land cover (LC) types [29,30]. Notably, the Landsat series, with freely accessible multispectral and Thermal Infrared Sensor (TIRS) imagery, has enabled the comprehensive spatiotemporal monitoring of LULC changes and UHI effects over the past 51 years [31,32,33,34,35]. Numerous studies have examined the spatiotemporal variability of LST and its key drivers, including urbanization levels [14], urban landscape characteristics [36], canopy density [37], and LC types [11]. Recent research has applied various spectral LC indices to explore their relationship with LST [14,33,38,39]. Automated GIS-based indices effectively identify built-up areas [40], bare lands [41], vegetation [42], water bodies [43,44], and wetlands [45]. It is widely agreed that these underlying LC types are crucial in the formation of UHIs [18,31]. For example, natural water bodies can substantially reduce the UHI effect [46,47], whereas urban areas tend to absorb more heat due to the presence of tall buildings, which store, retain, and generate heat [9,10].

While GIS offers an efficient platform for integrating and analyzing multitemporal LULC changes and associated LST in rapidly urbanizing areas [48], offering critical insights for UHIs effects, and quantifying LULC impacts on UHI intensity [14,15,49]. However, the assessment of specific LC changes on LST variations in rapidly urbanizing cities facing shifting climate patterns over extended periods remains insufficiently explored. This research presents a novel approach, offering a comprehensive analysis of the thermal impacts of different LC types on LST fluctuations, highlighting both cooling and warming effects, and their influence on local air temperatures in a densely populated riverine city.

South Asian developing countries face significant challenges in achieving sustainable urban development due to increasing poverty, greenhouse gas emissions, and energy demands [50,51,52]. With the regional population expected to reach 2.1 billion by 2030, those earning less than USD 2 per day remain particularly vulnerable to climate change impacts [21,50]. Pakistan, with 36.4% of its population residing in urban areas, has the highest urbanization rate in South Asia [25]. This rapid urban growth amplifies the UHI effect in its arid and semi-arid regions (Figure 1a), leading to higher energy demands [53]. In 2021, the World Bank identified Pakistan as highly vulnerable to rising temperatures with projections of up to a 2.5 °C increase by 2050 [54]. Studies show that rapid urbanization in Pakistan has increased LST due to the expansion of built-up areas and impervious surfaces. Research on Islamabad (1993–2018) indicates significant urban growth, higher LST, and reduced thermal comfort [11,55,56]. Similarly, urbanization in Lahore over three decades has elevated LST, highlighting the need for green infrastructure to mitigate thermal stress [57]. Karachi has also experienced heightened UHI effects due to land-use changes, emphasizing the importance of monitoring urban growth and its ecological impacts [58]. These findings align with global research, linking urbanization and land-use changes to thermal environmental alterations, underscoring the need for sustainable urban planning and green infrastructure. Islamabad, the capital city of Pakistan in northern Punjab (Figure 1b,c), was chosen as the case study, and it experienced a 3.5 °C temperature rise between 1960 and 2010 [59,60]. The city has undergone significant changes in LC, which are driven by rapid urbanization, population growth, industrialization, and real estate development [26,49,61,62,63].

While some studies have exclusively focused on the dynamic changes in LULC and their effects on UHIs in Islamabad, this research shed insights on the thermal influence of various LC types on LST-based temperature classes, encompassing both cooling and warming effects, and their impact on local air temperature under anthropic interventions and erratic rainfall patterns. The study integrates Landsat imagery (7 and 8), ground-based meteorological data, and GIS technologies to (1) investigate the spatiotemporal changes in LULC from 2000 to 2020, identifying key driving factors and their influence across different LST ranges and (2) evaluate the contributions of different LC types to warming and cooling effects in the period 2000–2020, along with the resultant impacts on localized air temperature variations under variable precipitation patterns, as derived from ground observations. This research explores the impact of driving forces, such as human activities and population growth, on UHIs and related environmental degradation in Islamabad. The findings aim to improve understanding of how urban planners can mitigate anthropogenic stressors, aiding decision-makers in formulating effective climate adaptation strategies.

2. Materials and Methods

2.1. Study Area Characteristics

The study area spans approximately 906 km², encompassing Islamabad, which was strategically developed as the capital city. Located at 33°28′ N and 72°48′ E, between the Potohar Plateau and Rawalpindi district, it lies at an elevation of 540 m (Figure 1). Attracting migrants due to better living standards, education, and job opportunities [18], the city experiences an urban growth rate of 3.48% [64]. The city’s climate is classified as humid subtropical [49]. June and July are the hottest months, with average temperatures reaching 36 °C, while December and January are the coldest with average monthly temperatures ranging from 2 to 3 °C. The monsoon season occurs in July and August, which is characterized by heavy rainfall.

2.2. Research Methodology

Figure 2 illustrates the comprehensive methodological framework of this research. This study employed Landsat (7 and 8) satellite imagery, meteorological data (temperature and rainfall), and population census data spanning two decades (2000–2020). Detailed explanations are provided in the subsequent sections.

2.2.1. Data Collection

Medium spatial resolution Landsat 7 (L7) and 8 (L8) imagery were freely obtained from the United States Geological Survey (USGS) Earth Explorer platform for the years 2000, 2010, and 2020. The study employed satellite sensors from L7 and L8, specifically the ETM+ and the OLI, along with Thermal Infrared Sensors (TIRSs), to examine the correlation between LULC and LST. Cloud-free (less than 5%) scenes were acquired during the summer season in June. L7 imagery was utilized for 2000 and 2010, while Landsat 8 imagery was employed for 2020, and their main qualifications are outlined in

Table 1. Description of dataset characteristics and related sources.

Data Type	Datasets	Date	Specification	Source
Ancillary Data	Meteorological Station	2000, 2010, and 2020	Monthly air temperature and rainfall data for the years 2000, 2010, and 2020.	Pakistan Meteorological Department (PMD).
	Land cover	2000, 2010, and 2020	Land cover types for the LULC classification for assessment.	Google Earth imagery
	Population Census	2000–2017	Populace number in Islamabad.	Pakistan Bureau of Statistics (https://www.pbs.gov.pk/content/final-results-census-2017-0, accessed on 12 December 2022)
Remote sensingdata	Landsat 7 ETM+ (L7)	14 June 2000	One scene (Path/Row = 150/37). Band 1—Blue Band 2—Green Band 3—Red Band 4—near infrared (NIR) Band 5—Short infrared (SWIR-1) Band 7—SWIR-2 Band 6—Thermal (60 m) Spatial Resolution 30 m	https://earthexplorer.usgs.gov/ (accessed on 12 December 2022)
		26 June 2010
	Landsat-8 OLI (L8)	26 June 2020	One scene (Path/Row = 150/37). Band 2—Blue (B) Band 3—Green (G) Band 4—Red (R) Band 5—near infrared (NIR) Band 6—Short infrared (SWIR-1) Band 7—SWIR-2 Band 10—Thermal Infrared (TIR1) (60 m) Spatial Resolution 30 m

. Scan-Line Errors were employed for all visible, infrared, and thermal bands of L7 images using the Landsat Toolbox in ArcMap version 10.5. All bands from the L7 and L8 scenes were subsequently calibrated to reflectance using radiometric calibration and atmospheric correction tools within the ERDAS software version 16.

The ancillary geospatial data included ground meteorological monthly records (Appendix A, Figure A1), Google Earth imagery, and population data spanning the years 2000 to 2020. Particularly, monthly air temperature and rainfall data for Islamabad for the years 2000, 2010, and 2020 were obtained from the Pakistan Meteorological Department (PMD). These measurements were recorded at the weather station located near Zero Point in Islamabad. Population data were obtained from the official website of the Pakistan Bureau of Statistics for the census years 1998, 2017, and 2020 [65,66]. All data were integrated into a GIS-based relational database implemented in ArcGIS to streamline subsequent image classification, accuracy assessment, and analysis.

2.2.2. Land Use/Land Cover (LULC) Retrieval and Classification

Notably, only the non-thermal bands—specifically bands 1, 2, 3, 4, 5, and 7 from L7 and 2, 3, 4, 5, 6, 7 from L8—were selected and subset to the investigated area. A contrast stretching technique was then applied, which was followed by the generation of false-color composite (FCC) images to enhance feature visibility (Figure A2). Supervised classification was employed to monitor LULC classes using ERDAS Imagine 2015 and ArcMap 10.5. The Maximum Likelihood Classifier (MLC) was applied to generate LULC maps for the three time periods (2000, 2010, and 2020). Therefore, five geospatial land-use indices were derived from the corrected Landsat bands to retrieve specific LC types (

Table A1. Spectral indices used for Landsat images.

Name and Abbreviation	Equation	Reference
Normalized Difference Vegetation Index (NDVI)	$\mathrm{NDVI}={\displaystyle \frac{\left(\mathrm{NIR}-\mathrm{RED}\right)}{\left(\mathrm{NIR}+\mathrm{RED}\right)}}$	[42]
Normalized Difference Water Index (NDWI)	$\mathrm{NDWI}={\displaystyle \frac{\left(\mathrm{Green}-\mathrm{NIR}\right)}{\left(\mathrm{Green}+\mathrm{NIR}\right)}}$	[61]
Normalized Difference Build-up Index (NDBI)	$\mathrm{NDBI}={\displaystyle \frac{\left(\mathrm{SWIR}1-\mathrm{NIR}\right)}{\left(\mathrm{SWIR}1+\mathrm{NIR}\right)}}$	[62]
Normalized Difference Bareness Index (NDBaI)	$\mathrm{NDBaI}={\displaystyle \frac{\left(\mathrm{SWIR}1-\mathrm{TIR}1\right)}{\left(\mathrm{SWIR}1+\mathrm{TIR}1\right)}}$	[63]

), including vegetation, water bodies, soil, and built-up areas, as well as to assess urban heat dynamics in Islamabad. Spectral indices were independently classified using predefined threshold values from the relevant literature to enhance accuracy in identifying regions of interest (ROIs). A random sample comprising 80% of the ROIs database was selected to train the classification of the LULC maps for each year. The reliability of the classified images was assessed by comparing 20% of randomly selected ROIs with the corresponding classified pixels using a confusion matrix (

Table A2. Accuracy assessments of the classified images of the years 2000, 2010, and 2020.

Years		Reference Data					Classified Accuracy %
	Classes	Vegetation	Build-Up Areas	Water Bodies	Bare Land	Producer’s Accuracy	User’s Accuracy
2000	Vegetation	39	0	5	2	97.50	85
	Built-Up Areas	1	40	4	0	100	89
	Water Bodies	0	0	31	0	77.50	100
	Bare Land	0	0	0	38	95	100
	Overall Accuracy (%)	92.5
	Kappa Coefficient	0.90
2010	Vegetation	40	3	4	2	100.00	81.63
	Built-Up Areas	0	36	1	1	90	95
	Water Bodies	0	0	35	0	87.5	100
	Bare Land	0	1	0	37	92.5	97.40
	Overall Accuracy (%)	92.5
	Kappa Coefficient	0.90
2020	Vegetation	40	2	1	4	100	89
	Built-Up Areas	0	36	5	1	90	86
	Water Bodies	0	0	34	0	85.5	100
	Bare Land	0	2	0	35	87.5	95
	Overall Accuracy (%)	90.6
	Kappa Coefficient	0.88

). Subsequently, the study area was classified into four categories: vegetation, water bodies, bare land, and built-up areas. Lastly, a post-classification change detection method was employed to pinpoint pixel-level differences between the independently classified images from 2000 and 2020, utilizing the intersect tool in ArcMap [67].

To this end, the Normalized Difference Vegetation Index (NDVI), which leverages the spectral contrast between the red and near-infrared (NIR) bands, was used to assess the greenness of healthy and natural vegetation across the study area [42]. The Normalized Difference Water Index (NDWI) was developed to assess the spatial variability of surface wetness and water bodies using near-infrared (NIR) and green spectral bands [68]. The NDBI is commonly employed to map and monitor changes in built-up areas [38,69] and is calculated using NIR and SWIR bands. The Normalized Difference Bareness Index (NDBaI) is used to classify bare lands and open barren spaces [70]. It is determined by the difference between SWIR and Thermal Infrared (TIR1) bands.

2.2.3. Retrieving Land Surface Temperature (LST)

The LST was obtained through a three-step process: calculating brightness temperature, estimating vegetation proportion, and adjusting emissivity using the NDVI. Initially, the digital numbers (DNs) in the thermal band of the satellite imagery were converted into spectral radiance using Equation (1).

$$L_{\lambda }=\left({\displaystyle \frac{L_{Max \lambda }–L_{Min \lambda }}{Q_{calMax}–Q_{calMin}}}\right)\left(Q_{cal} –Q_{calMin}\right)+L_{Min \lambda }$$

(1)

where $L_{\lambda }$ = spectral radiance, $L_{Max \lambda }$ = spectral radiance of satellite sensor scaled to $Q_{calMax}$, $L_{Min \lambda }$ = spectral radiance of satellite sensor scaled to $Q_{calMin}$.

Second, the spectral radiance was determined using the DN of thermal band 6 from ETM+ and band 10 from OLI, as outlined in Equation (2).

$$L_{\lambda }=M_L .Q_{cal}+A_L$$

(2)

where $L_{\lambda }$ = atmospheric spectral radiance, $M_L$ = band-specific multiplicative rescaling factor, $A_L$ = band-specific additive rescaling factor, $Q_{cal}$ = quantized and calibrated standard product pixel values (DN).

After that, the brightness temperature (T B) was obtained from the top of the atmosphere (TOA) reflectance using a calibration constant (

Table A3. Value of calibration constant ($K_1$ and $K_2$) for Landsat datasets.

S. No.	Constant	Landsat 7 (ETM+)	Landsat 8 (OLI/TIRS)
			Band 10	Band 11
1	$K_1$	666.09	774.89	480.89
2	$K_2$	1282.71	1321.08	1201.14

), as presented in Equation (3).

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

(3)

where $T_B$ = brightness temperature in Kelvin, $L_{\lambda }$ = spectral radiance, and $K_1$ and $K_2$ = calibration constant given in Table 1 from the metadata.

The emissivity correction ($P_{\upsilon }$) was completed by applying the NDVI threshold method using Equation (4).

$$P_{\upsilon }=\left({\displaystyle \frac{\mathrm{NDVI}-{\mathrm{NDVI}}_{\min }}{{\mathrm{NDVI}}_{\max }-{\mathrm{NDVI}}_{\min }}}\right)$$

(4)

where the emissivity correction (Ԑ), based on the proportion of vegetation ($P_{\upsilon }$), was computed using Equation (5). Subsequently, LST was determined using Equation (6).

$$\epsilon =0.004P_{\upsilon }+0.986$$

(5)

$$T={\displaystyle \frac{T_B}{1+\left(\lambda \frac{T_B}{P}\right)\ln \epsilon }}$$

(6)

where $T_B$ = brightness Temperature, $\lambda$ = wavelength of radiance emitted, and $P=h\times c/s$, where h is Planck’s constant (6.624 × 10⁻³⁴ J s), c is the velocity of the light (2.998 $\times$ 108 m/s), and s is the Boltzmann constant (1.38 $\times$ 10⁻²³ J/K).

The resulting LST values are in Kelvin (T) (${\mathrm{LST}}_K$), which can be converted to Celsius (${\mathrm{LST}}_{°\mathrm{C}}$) using the relation where 0 degrees Celsius equals 273.15 K; this conversion is performed using Equation (7).

$${\mathrm{LST}}_K – 273.15={\mathrm{LST}}_{°\mathrm{C}}$$

(7)

2.2.4. Categorization of LST-Based Temperature Classes

To assess UHIs, a phenomenon associated with LST, each LST raster image was categorized into seven distinct classes using the quantile method in ArcMap [48]. The quantile method ranks all pixels in a specific order and divides them based on percentage intervals [71]. This approach not only reduces comparison errors between different dates but also tends to prevent the overclassification of areas as moderate, avoiding an overrepresentation of smaller regions as hot zones [48]. In this study, the LST values were classified into seven categories: extremely low, low, sub-low, medium, sub-high, high, and extremely high temperatures. The classes were denoted as Ltc_1 through Ltc_7 with Ltc_1 corresponding to the extremely low-temperature category and Ltc_7 representing the extremely high-temperature category. Subsequently, the proportions of underlying LC of different years were determined in each temperature-based LST class using the intersect tool in ArcMap.

2.2.5. Assessment of Land Cover Contributions to Temperature Patterns

The Land Contribution Index (LCI) was employed as a quantitative measure to evaluate the thermal influence of dominant land cover types on temperature variations within a specified region, as proposed by Huang et al. [72]. The LCI was utilized to quantify the thermal contribution of each LC type across different UHI-based LST classes in Islamabad. It counts the temperature difference of individual LC and its proportion to the area into consideration using Equation (8):

$$\mathrm{LCI}= \left({\mathrm{LST}}_i-{\mathrm{LST}}_m\right)\times A_i$$

(8)

where ${\mathrm{LST}}_i$ denotes the average temperature of the $i$-th land cover, ${\mathrm{LST}}_m$ represents the overall average LST of the study area, $i$ refers to the four LC types, and $A_i$ indicates the proportion of the i-th area relative to the total area. The mean temperature for each overlapping region between LC types and temperature categories across multiple years was determined using the “Zonal Statistics as Table” function within the Spatial Analyst tool in ArcMap.

Using the year 2000 as a reference point, the overlap between each LST class from 2000 and 2020 was calculated through the intersect tool to produce the Warming and Cooling Transition (WCT) map. WCT map were categorized into three groups: no change (same class), hot-to-cold, and cold-to-hot transitions. The cold-to-hot group represents areas that were in a lower temperature class in 2000 but converted to a higher temperature class by 2020, while the hot-to-cold group shows the reverse pattern. Additionally, the proportion of each LC type within each overlapping area was calculated to investigate the causes of these differences.

2.2.6. Assessment of Air Temperature and Rainfall Patterns

Rainfall significantly contributes to cooling processes through evaporation and perspiration, making it essential to consider the interaction between temperature and humidity, especially during extreme heat events. In this work, monthly minimum and maximum air temperatures, along with the annual average air temperature, were calculated for the years 2000, 2010, and 2020 at the weather station situated near Zero Point in Islamabad (Table 1).

3. Results

3.1. Spatiotemporal LULC Dynamics

Spectral indices were independently classified using predefined threshold values from to identify ROIs. For example, the vegetation class, representing both healthy and sparse plants, exhibited high positive NDVI values (>0) in the eastern parts and the northern Margalla hills (Figure 3). Notably, NDWI values were higher (NDWI > 0) in areas with surface water coverage, particularly in the central and western regions. The NDBaI index proved highly effective in delineating bare lands throughout the 2000–2020 period (Figure 3). Additionally, the analysis of NDBaI maps revealed a significant negative correlation with NDVI in vegetated areas. The NDBI index identified built-up (roads and buildings) areas with values greater than 0.1 across the investigated periods. Subsequently, the study area was classified into four categories: vegetation, water bodies, bare land, and built-up areas (Figure 4 and

Table 2. The area and percentages of change in different LULC categories from 2000 to 2020 in Islamabad.

Land Cover Classification	2000		2010		2020
	(km2)	(%)	(km2)	(%)	(km2)	(%)
Bare lands	504.06	55.70	399.68	44.16	94.33	10.42
Built-up areas	114.47	12.65	156.34	17.28	280.52	31.00
Vegetation	280.05	30.94	345.36	38.16	520.16	57.48
Water bodies	6.44	0.71	3.62	0.40	9.98	1.10

). The generated LULC maps were found to be reliable, achieving an overall accuracy exceeding 85% with no user or producer accuracy for any class falling below 70%. The Kappa coefficient was close to 1, indicating a high level of agreement [73]. Table A2 provides the overall accuracy and Kappa coefficient for the classified maps corresponding to the years 2000, 2010, and 2020. The results demonstrated both the negative effects of urban sprawl and the positive effects of land reclamation on most LC changes (Figure 4 and Figure 5). The area and proportions of the LULC classes were calculated as shown in Table 2, while

Table 3. Gains and losses analysis of different intervals for land cover classes.

Classes	Area Change
	2000–2010		2010–2020		2000–2020
	km2	%	km2	%	km2	%
Bare lands	−104.38	−20.71	−305.35	−76.40	−409.73	−81.29
Built-up areas	41.87	36.58	124.18	79.43	166.05	145.06
Vegetation	65.31	23.32	174.81	50.62	240.12	85.74
Water bodies	−2.82	−43.77	6.35	175.49	3.54	54.97

presents the net change for each LULC class across different periods, with proportions calculated relative to the area of the preceding date (Figure 5).

Bare lands play a crucial role in urban planning, especially in land use and reclamation endeavors. From 2000 to 2020, bare lands areas diminished from 504.06 km² to 94.33 km², representing a 81.29% decrease. This decrease was primarily driven by the conversion of 190.45 km², 226.06 km², and 3.02 km² of bare lands into built-up areas, vegetation, and water bodies, respectively, predominantly in the western, northwestern, and southwestern regions (Figure 5c and Table 3 and

Table 4. LULC transitions matrix (km²) of Islamabad from 2000 to 2020.

Classes		2020
		Bare Lands	Built-Up	Vegetation	Water Bodies
2000	Bare lands	84.33	190.45	226.06	3.02
	Build-Up	2.64	69.44	41.70	0.65
	Vegetation	7.02	20.24	251.27	1.45
	Water bodies	0.30	0.30	0.98	4.86

Values in bold refer to the unchanged area for each class.

). For example, urban expansion efforts were successfully implemented over bare lands in areas between Tarnol, Islamabad, Margalla Hills, and Barab Koh districts from 2000 to 2010 (Figure 5a). Similar activities, coupled with intensive agricultural reclamation, were also observed in the southern Nilore and Ali Pur districts between 2010 and 2020 (Figure 5b).

Unplanned urban expansion, particularly in the northeastern regions like Islamabad, Margalla Hills, and Barab Koh, led to the degradation of nearby vegetated areas from 2000 to 2010 (Figure 5a). The vegetation cover increased significantly from 280.05 km² in 2000 to 520.16 km² in 2020, reflecting a 85.74% (240.12 km²) growth over two decades. This expansion primarily resulted from the conversion of 226.06 km² of bare lands, 41.7 km² of built-up areas, and 0.98 km² of water bodies into vegetated land (Figure 5c; Table 3 and Table 4). Meanwhile, water bodies decreased by 43.77% (−2.82 km²) between 2000 and 2010, then grew substantially by 175.49% (6.35 km²) between 2010 and 2020, resulting in an overall increase of 54.97% (3.54 km²) with the most significant growth observed in the central areas around Rawal Lake (Figure 5c).

3.2. Spatial Distribution of LST-Based Temperature Classes

Descriptive statistics for LST in Islamabad, including minimum, maximum, mean, and standard deviation, were computed using the data analysis add-in in Microsoft Excel and are summarized in

Table 5. Descriptive statistics of LST for the study area.

Indices	Year	Min	Max	Mean	Standard Deviation
LST	2000	21.32	32.39	27.71	2.70
	2010	26.20	43.21	37.22	3.58
	2020	25.64	45.63	38.12	4.36

. A notable rise in LST is evident, from 32.39 °C in 2000 to 45.63 °C in 2020. The cold temperature classes (Ltc_1 to Ltc_3) were primarily situated in the northern Margalla Hills, Rawal Lake, and other water bodies in the western portion of the study area throughout this period (Figure 6). While these cold zones in Margalla maintained relatively lower temperatures, higher heat temperature classes (Ltc_5 to Ltc_7) were predominantly found in substantial areas of the city in 2000, leading to increased thermal stress (Figure 6a). By 2010, a significant temperature rise was observed in the western and southeastern regions (Figure 6b). In 2020, the spatial distribution of higher heat temperature classes intensified in the western and southern regions of the study area. Conversely, colder temperature classes were predominantly observed in the eastern regions, encompassing the Barab Koh, Nilore, and Ali Pur districts (Figure 6c).

3.3. Land Cover Contributions on LST Patterns

Figure 7 illustrates the distribution of various LC types and their proportions across different temperature-based LST classes. As temperatures rise, the proportion of vegetation decreases progressively within each temperature class, while the proportions of bare lands and built-up areas increase (Figure 7). In the lowest temperature class, vegetation constitutes the largest proportion, exceeding 95%, followed by urban areas, which account for approximately 8–18%. Bare lands contribute the smallest proportion, remaining below 5% in both 2000 and 2010. Notably, water bodies are present only in Ltc_1 to Ltc_4 classes with the proportion of water bodies in Ltc_1 being higher in 2020 than in the other years. In the highest temperature class, the proportion of built-up areas is lowest, ranging from 10 to 55%, while bare lands exceed 70%, and vegetation comprises about 10% (Figure 7). In this context, spectral indices, including NDVI, NDWI, NDBI, and NDBaI, were used to assess the spatiotemporal relationship between LULC changes and UHI intensity. Figure 8 presents scatter plots showing correlations between these indices and LST, revealing a negative correlation between NDVI and NDWI with LST, and a positive correlation between NDBI and NDBaI with LST. These findings highlight the significant impact of urban expansion and development on LST and UHI intensity.

Figure 9 represents the LCI of each LULC class against each temperature class in the entire area from 2000 to 2020. An LCI value of ≥0 indicates that the corresponding LULC class contributes to a temperature increase in the area, whereas an LCI value of <0 suggests that the LULC class contributes to a temperature decrease in the region. In the Ltc_1 class, the LCI for vegetation, water bodies, bare lands, and built-up areas exhibited negative values in 2000 and 2010, initially contributing to cooling. Notably, vegetation demonstrated the least negative LCI among these categories, reaching up to −0.45 for vegetation cover. By 2020, only vegetation and water bodies retained negative LCI values, suggesting their continued positive impact on temperature reduction over time. Between 2000 and 2010, the LCI of vegetation, water bodies, bare lands, and built-up areas showed negative LCI values from the Ltc_2 to Ltc_3 classes. However, these values became positive in Ltc_4, slightly contributing to a temperature increase in the class. The vegetation and built-up areas played a significant role in mitigating heat in the Ltc_2 to Ltc_4 regions. Meanwhile, the LCI for water bodies does not exist from Ltc_2 to Ltc_7 classes. The LCI for bare land and built-up areas is notably higher in the Ltc_5 and Ltc_6 classes compared to the lower temperature classes. Within the Ltc_7 class, bare lands exhibited the most pronounced positive temperature in 2000 and 2010, while built-up areas recorded the greatest positive values in 2020. Overall, the LCI of vegetation and water bodies consistently exert a cooling effect in lower temperature zones (Ltc_1 to Ltc_3), while bare land and built-up areas show a weaker cooling impact during the study period. In higher temperature zones (Ltc_5 to Ltc_7), bare land and built-up areas significantly contribute to regional warming. Over the years, the highest LCI in 2000 was observed for Ltc_1 and Ltc_5 to Ltc_7, with Ltc_4 peaking in 2010 and Ltc_2 and Ltc_3 in 2020.

3.4. Spatiotemporal Land-Cover Warming and Cooling Conditions

Figure 10 illustrates the WCT map, highlighting spatial changes in regions where LST classes overlap from 2000 to 2020. These overlapping areas are concentrated within 42 adjacent categories, which are evenly distributed between cold-to-hot and hot-to-cold temperature transitions as defined in this study. The spatial distribution of temperature increased over the study period (Figure 10), particularly in newly developed urban areas such as western Tarnol, Islamabad, Barar Koh, Ali Pur, and Sihala districts. Consequently, regions with overlapping cold and hot areas experience a lower warming degree compared to other regions (Figure 6a). This has led to an upgrade in their LST classification by 2020 (Figure 6b). As depicted in Figure 9, the presence of built-up areas and bare land contributes significantly to surface heating. In contrast, the northern regions, characterized by high terrains and lakes, such as Margalla Hills and Rawal Lake, have maintained relatively stable temperature conditions. On the contrary, areas with overlapping hot and cold characteristics exhibited a relatively high warming degree in 2000, particularly in the eastern and middle parts, resulting in a downgrade of their temperature classification (Figure 10).

Figure 11 shows the proportion of each land cover type in the overlapping areas of the cold-to-hot and hot-to-cold groups over the studied period. Reasonably, cold-to-hot transitions occur more frequently in the lower temperature range, whereas hot-to-cold transitions are more prevalent in the higher temperature range (Figure 11). In the hot–cold overlap category, the proportion of vegetation was notably higher than that of other LC types during transitions from higher to lower temperature classes except for transitions between Ltc_6 to Ltc_5, Ltc_7 to Ltc_5, and Ltc_7 to Ltc_6 classes (Figure 11a). In the cold–hot overlap category, the proportion of bare land and built-up areas significantly increased as the transitions moved from lower to higher temperature classes (Figure 11b). The proportion of vegetation is greater in the Ltc_2 to Ltc_1 class compared to other transition classes.

These overlapping groups reflect overall warming and cooling trends from 2000 to 2020. The cooling effect in the hot–cold overlap area, which results in a reduction in LST classification over the study period, is higher than in other regions. Conversely, the cold–hot overlap area exhibits a stronger warming trend, leading to an increase in its LST classification. When examining the proportion of each land cover type in the overlapping cooling and warming categories (

Table 6. The percentage of intersected land cover with different cold–hot transition types.

	Cold–Hot (%)	Hot–Cold (%)
Bare lands	14.72	6.97
Built-up areas	48.89	19.00
Vegetation	35.80	73.45
Water bodies	0.59	0.58

), it is evident that built-up areas constitute the largest proportion of the cold–hot category, accounting for approximately 48.89%. This proportion is significantly greater than the 19% observed in the hot–cold category, suggesting a stronger warming effect associated with built-up areas. Additionally, bare lands and water bodies exhibit a higher proportion in the cold–hot category compared to the hot–cold category, indicating their greater contribution to warming. In contrast, vegetation shows a notably higher proportion in the hot–cold category, indicating a more pronounced cooling effect.

3.5. The Influence of Rainfall Variability on LST and Air Temperature

The LST was compared with air temperature data from the PMD in Islamabad for the years 2000, 2010, and 2020 (Figure A1). Meteorological data indicated that the highest average temperatures occurred in May and June in 2000 and 2010, respectively, while July recorded the highest temperatures in 2020. This shift indicates a transition of the hottest months from May and June to an extension into July. Interestingly, the estimated LST values closely aligned with the ground measurements from PMD for 2010 and 2020, ranging from 35 to 39 °C (Figure 12). Figure 12 presents the comparison between LST derived from satellite data and air temperature obtained from ground measurements. The maximum temperature difference between the mean annual land-based air temperature and LST was 11.55 °C in 2000, decreasing to 2.28 °C in 2020, indicating a significant reduction in monthly air temperatures over the observed years. To further explore the reasons behind this decrease in air temperature, monthly rainfall patterns for 2000, 2010, and 2020 were analyzed (Figure 13). The monsoon season in Islamabad typically spans from June to September, with the highest rainfall recorded at the end of June, which was followed by July and August in 2010. In contrast, a more consistent rainfall pattern was observed throughout the year in 2020, with total annual rainfall reaching 1365.47 mm, which was considerably higher than the 982.88 mm recorded in 2000. Notably, August 2000 and March 2020 experienced the highest rainfall in Islamabad.

4. Discussion

Islamabad serves as an exemplary case of a rapidly expanding urban environment that has undergone significant LC transformations driven by urbanization, population growth, industrialization, and real estate development [26,49,61,62]. While previous studies have highlighted that differences in LST primarily arise from variations in LULC dynamics [11,18,49,60,61,62,74], this study further explores the relationship between population growth and vegetation expansion. The findings of the study identified the drastic effects of rapid population growth and urbanization on changes in land use in Islamabad. As indicated in Table 3, the expansion of built-up areas has resulted in a significant decrease in bare land cover, which is a trend probably driven by Islamabad’s increasing population. The population of Islamabad rose from 568,689 to 804,337 between 2000 and 2010 [65], and it subsequently approached 1,129,198 by 2020 [66]. These data clearly illustrate substantial growth in both population and built-up areas. Over the two decades, the population nearly doubled, exhibiting a remarkable increase of 49.63%. Similarly, the built-up area expanded significantly, from 114.47 km² in 2000 to 156.34 km² by 2010, and further to 280.52 km² in 2020, resulting in an overall increase of 145.06% (166.05 km²) during the study period (Table 3). These changes in land use, along with the construction of multistory high-rise buildings, have profoundly transformed the city’s linear layout toward western, northeastern, eastern, and southeastern regions (Figure 4c). This development has extended not only into suburban areas along major highways but also into established urban sectors. The city has experienced a significant transformation, transitioning from predominantly low-density, horizontal expansion to a denser, vertically oriented development pattern [75]. Such trends corroborate findings from prior research [76,77], which documented the replacement of low- and mid-rise structures with high-rise commercial and residential complexes, [77], reshaping Islamabad’s landscape [74].

Analysis of the distribution of LC types across temperature-based LST classes revealed that as temperatures rise each year, vegetation cover decreases, which is accompanied by an increase in bare land and built-up areas (Figure 7). For instance, bare lands areas consistently show higher LST values primarily due to extensive urban development (Figure 9), which involves paving the land for construction projects, particularly along the city’s periphery [30]. Correlations between normalized indices confirmed these observations: negative relationships for NDVI and NDWI and positive relationships for NDBI and NDBaI with LST indicate the significant role of urban expansion in intensifying LST and UHI intensity. Persistent high temperatures in bare land areas underscore the necessity of low-impact development practices to mitigate urbanization-induced warming [18,77]. Although vegetation has a relatively lesser impact on LST within highly vegetated areas (Ltc_7 class), deforestation and vegetation loss have exacerbated local climatic conditions and intensified UHI effects [11,14,27,28]. Comparable patterns are observed in other rapidly urbanizing cities, where reduced green spaces have amplified UHI impacts [78]. Vegetation, through shading and evapotranspiration, reduces temperatures by 2–4 °C in highly vegetated areas such as parks and greenbelts in Islamabad [79]. Similarly, water bodies like Rawal Lake contribute to localized cooling through evaporation, functioning as natural heat sinks [80].

During the study period, regions dominated by built-up areas and bare land exhibited the highest temperature increases, while vegetation and water bodies provided cooling effects, highlighting their critical role in mitigating UHI impacts through urban planning interventions (Figure 10 and Figure 11 and Table 6). The slight increase in water bodies observed during the study period can be attributed to the heavy rainfall recorded in 2020 (Figure 13). For example, the highest monsoon rainfall was documented in June, July, and August 2010, which led to severe flooding in subsequent years [81]. In 2020, consistent rainfall was observed throughout the year with a total of 1365.47 mm—substantially higher than the 982.88 mm recorded in 2000. This increased rainfall corresponded with notable growth in vegetation and water bodies by 2020 (Figure 5c), contributing to a decline in air temperature in Islamabad (Figure 12). These findings highlight the importance of water bodies, particularly in mitigating lower temperature classes [46,47], and emphasize the potential for urban water management to reduce UHI intensity. Such features can play a pivotal role in minimizing UHI effects in urban regions [47]. Despite this, the current extent of water bodies remains insufficient to satisfy the city’s growing demands [26]. Consequently, rapid urban expansion and projected population growth are expected to intensify pressures on agricultural lands and groundwater resources. These findings corroborate earlier studies on the climatic impacts of urbanization and LULC changes. For instance, [82] explored the impact of urbanization-induced LULC changes, along with adjacent river morphological alterations, on urban microclimates. The study concluded that river migration increased heat-emitting lowlands, resulting in reduced annual rainfall and wind speeds, as well as elevated relative humidity. Similarly, Ullah et al. [83] conducted a detailed investigation of land use trends and the impacts of climate change on the hydrological systems of Pakistan’s Ravi River basin. Their results highlighted a significant increase in urbanized areas from 1990 to 2020 along with projected rises in rainfall and temperature under future climate scenarios. Rao et al. [84] analyzed three decades of data from the Emilia–Romagna region in Italy and identified buildings and vegetation as the primary land cover types influencing surface temperatures, reaffirming their significant role in the formation of urban heat islands.

These findings highlight the intricate relationship between population growth, human activities, and land-use changes during the study period, providing valuable insights for sustainable resource management in Islamabad, significantly identifying vulnerable areas and evaluating the ecological and environmental health of the city. Transforming these susceptible locations into green spaces, such as parks, green roofs, and landscaped areas around buildings, as well as incorporating water bodies, could improve ecological conditions [36,85]. Future studies should quantify the long-term effects of green and blue infrastructure on reducing LST in Islamabad. For instance, the WCT zoning map should be used to prioritize the preservation and enhancement of green spaces while promoting sustainable urban expansion. Recent initiatives, such as urban afforestation policies and the development of greenbelts using Miyawaki techniques, have enhanced ecological services in Islamabad, including carbon sequestration and temperature regulation [86]. This relatively effective integration of urban growth with environmental conservation offers a promising framework for sustainable urban development that is aligned with international best practices [87,88]. Additionally, the potential inclusion of urban wetlands and rainwater retention ponds should be further explored as sustainable solutions to mitigate UHIs effects and manage water resources in the city.

5. Limitations and Future Research

Despite these advancements, several limitations must be acknowledged. First, the study selected the years 2000, 2010, and 2020 as representative of summer seasons. However, this approach may introduce uncertainty due to the potential influence of unpredictable weather events, such as those observed in 2020. Second, land cover classification was conducted using medium-resolution Landsat imagery with a spatial resolution of 30 m to align with the LST data. While this resolution is appropriate for broader analyses, it may not adequately capture the finer details of specific underlying surface characteristics. Third, the analysis primarily focused on the overall impact of vegetation types on LST variations. However, other critical variables, such as vegetation conditions (e.g., sparse vegetation, agricultural fields, forested areas) and biomass volume across different seasons, were not accounted for in this study. Shiflett et al. [89] argued that overlooking micro-advection processes within the near-surface air can lead to erroneous associations between air temperature and LST due to the complex relationship between these two factors and the generation of warming and cooling air fluxes across adjacent regions. To address these limitations and enhance the robustness of future studies, subsequent research will aim to characterize distinct land cover types, such as sparse vegetation, agricultural zones, forested areas, water bodies, impervious surfaces, and their surrounding environments. This will involve using monthly LULC data in conjunction with LST imagery processed through the Google Earth Engine platform, providing a more detailed and seasonally comprehensive assessment of land cover influences on LST variations within the study area.

6. Conclusions

Riverine cities in southern Asian developing countries encounter substantial challenges in understanding the spatiotemporal thermal effects of specific LULC types, which are driven by rapid urbanization under climate changes. This study successfully employed an integrated approach utilizing Landsat imagery (5, 8, and 9), ground-based meteorological data, and GIS to comprehensively analyze the thermal impact of various LC types on LST variations in Islamabad over the past two decades. The analysis revealed both cooling and warming effects on local air temperature amidst evolving precipitation patterns. Key findings include the following:

• LST rose from 32.39 °C in 2000 to 45.63 °C in 2020. A negative correlation was observed between NDVI and NDWI with LST and a positive correlation was observed between NDBI and NDBaI with LST.
• Negative LCI values indicated that vegetation and water bodies in lower-temperature zones (Ltc_1 to Ltc_3) contributed to cooling effects, whereas positive LCI values for built-up and bare land areas in higher-temperature zones (Ltc_5–Ltc_7) demonstrated warming effects.
• The WCT map showed a general warming trend (cold-to-hot transition) from 2000 to 2020 especially in newly urbanized areas due to a 49.63% population increase. In contrast, Margalla Hills and Rawal Lake exhibited stable temperatures, while newly developed agricultural lands showed cooling effects with a 46.46% rise in vegetation cover.
• Built-up areas accounted for 48.89% of the cold–hot category, greatly exceeding the 19% in the hot–cold category, indicating a stronger warming effect. Vegetation areas increased by 85.74%, accounting for 35.8% of the cold–hot transition compared to 73.45% in the hot–cold transition, reflecting their significant cooling effects.
• The mean annual air temperature and LST gap narrowed from 11.55 °C in 2000 to 2.28 °C in 2020. This reduction coincided with increased annual rainfall from 982.88 mm in 2000 to 1365.47 mm in 2020 and the expansion of water bodies from 2.82 km² in 2000 to 6.35 km² in 2020, influencing the local climate and hydrology.

The findings provide critical insights for identifying vulnerable zones and promoting urban ecological health. Expanding green spaces and water bodies while preserving agricultural lands can mitigate greenhouse gas impacts and urban heat, offering a sustainable planning model for similar riverine cities worldwide.