Assessment of Changes in Land Use/Land Cover and Land Surface Temperatures and Their Impact on Surface Urban Heat Island Phenomena in the Kathmandu Valley (1988–2018)

Abstract

More than half of the world’s populations now live in rapidly expanding urban and its surrounding areas. The consequences for Land Use/Land Cover (LULC) dynamics and Surface Urban Heat Island (SUHI) phenomena are poorly understood for many new cities. We explore this issue and their inter-relationship in the Kathmandu Valley, an area of roughly 694 km², at decadal intervals using April (summer) Landsat images of 1988, 1998, 2008, and 2018. LULC assessment was made using the Support Vector Machine algorithm. In the Kathmandu Valley, most land is either natural vegetation or agricultural land but in the study period there was a rapid expansion of impervious surfaces in urban areas. Impervious surfaces (IL) grew by 113.44 km² (16.34% of total area), natural vegetation (VL) by 6.07 km² (0.87% of total area), resulting in the loss of 118.29 km² area from agricultural land (17.03% of total area) during 1988–2018. At the same time, the average land surface temperature (LST) increased by nearly 5–7 °C in the city and nearly 3–5 °C at the city boundary. For different LULC classes, the highest mean LST increase during 1988–2018 was 7.11 °C for IL with the lowest being 3.18 °C for VL although there were some fluctuations during this time period. While open land only occupies a small proportion of the landscape, it usually had higher mean LST than all other LULC classes. There was a negative relationship both between LST and Normal Difference Vegetation Index (NDVI) and LST and Normal Difference Moisture Index (NDMI), respectively, and a positive relationship between LST and Normal Difference Built-up Index (NDBI). The result of an urban–rural gradient analysis showed there was sharp decrease of mean LST from the city center outwards to about 15 kms because the NDVI also sharply increased, especially in 2008 and 2018, which clearly shows a surface urban heat island effect. Further from the city center, around 20–25 kms, mean LST increased due to increased agriculture activity. The population of Kathmandu Valley was 2.88 million in 2016 and if the growth trend continues then it is predicted to reach 3.85 million by 2035. Consequently, to avoid the critical effects of increasing SUHI in Kathmandu it is essential to improve urban planning including the implementation of green city technologies.

1. Introduction

People increasingly prefer to live in urban and surrounding areas since these are perceived to offer improved opportunities for employment, good healthcare, and higher level of education [1,2,3]. As a result, the global urban population has increased to 55% by 2018, and is predicted to increase further to 68% by 2050 [4]. Such growth in the urban population has resulted in increased spread of built up urban areas which have directly or indirectly altered the dynamics of land use/land cover (LULC) at large scales [1,5,6]. Globally, land surface air temperature has increased by 1.53 °C and the global mean surface (ocean and land) air temperature by 0.87 °C since the pre-industrial period (1850–1900) [7]. The dramatic rise in urbanization has negative environmental consequences especially increases in Surface Urban Heat Island (SUHI) effects which severely impact human health, air quality, and climate change, all of which are major topics of global concern [8,9]. The observed higher land surface temperatures (LST) over city centers compared to surrounding areas, the so-called SUHI, results largely from the conversion of natural spaces into built-up areas with largely impervious surfaces [8].

Understand LULC has become increasingly important to help understand aspects of urban dynamics such as geography, morphology, ecology, and sustainability, with the goals of strengthening concepts of land use patterns, urban intensity, urban diversity, UHI phenomena, among other factors [9]. LULC of urban improvement and urban expansion have been explored widely in both developed and developing nations with spatial expansion appearing to be more complicated in developing nations [10]. Assessment of LULC in some cities has been particularly influential in aiding future urban planning, transport network development, economic prosperity growth, policy development, and environmental enhancement [10,11,12].

The surfaces of city landscapes have properties that lead to a wide variety of behaviors with respect to phenomena such as electromagnetic radiation absorbance or evaporation, longwave radiation, and prevailing winds, resulting in increased heat discharge [13]. The physical surface of a city comprises different materials such as asphalts, gravels, and stones as well as other building material surfaces, which enhance the sensitivity and lower evapotranspiration in the city with consequent effects on city climate [14,15]. This leads to higher absorption of sun radiation, higher retention of infrared radiation in street canyons, higher uptake, and delayed release of heat from buildings, higher proportion of absorbed sun radiation than latent heat forms, and higher release of latent heat from combustion fuels for industrial processing, urban transport, and domestic space heating [16]. As a result, core urban spaces experience warmer weather than their surrounding areas, which leads to the creation of SUHI. This SUHI phenomenon drastically transforms the local climate with adverse effects on plants as well as animals including humans [10,17]. SUHI has rapidly increased due to the large-scale increase in areas of impervious surfaces in the process of urbanization [18,19]. As a result, SUHI is of increasing concern for urban planners, health authorities, urban investors, policy makers, ecologists, and others due to its adverse effects on various aspects of the environment such as rainfall, temperature, air quality, energy balance, and carbon storage [20,21].

To help mitigate the effects of SUHI, the spatial dynamics of city landscapes need to be assessed [13]. Plans for improving cities should include improved design of building structures, including the use of environmentally friendly materials to minimize thermal energy storage, excessive impermeability, and solar radiance absorbance. Developing sound mitigation strategies will help ‘green’ cities, such as improving the energy balance through the use of cool and green roofs, increased urban tree plantings, and the development of gardens [22]. To aid the development of these strategies there is a need for spatial quantification of LULC and SUHI dynamics using Remote Sensing since these tools provide cost-effective and reliable data [23,24,25,26,27]. Elsewhere, assessment of LULC dynamics has provided land–resource information which has aided sustainable development [28,29,30,31]. Recently, the SUHI phenomena [32,33], the interrelationship between LST vs. Normal Difference Vegetation Index (NDVI) [34,35,36], LST vs. Normal Difference Built-up Index (NDBI) [37,38] and LST vs. Normal Difference Moisture Index (NDMI) [39] have been assessed in several city landscapes using Pearson’s correlation coefficients [40,41]. These studies confirm that impervious areas increase SUHI and increasing vegetation reduces this effect.

Studies of the growing urban history of Kathmandu in Nepal date back to the late 1950s [42,43] and showed that urban areas have expanded rapidly since this time. While the processes of urbanization and LULC change in Kathmandu has been previously explored [43,44,45,46,47], the role of LULC in LST intensification is not understood. Here we address this knowledge gap assessing the spatiotemporal dynamics of LULC and their influence on LST intensification, creation of SUHI, and the relationship of NDVI and NDBI with LST using Pearson’s correlation coefficient. We specifically examine these issues in the geopolitical area of Kathmandu Valley, which includes the districts of Bhaktapur, Kathmandu, and Lalitpur [42]. We used Landsat datasets (1988–2018) at decadal intervals selecting dry, summer season (April) data [48]. This study area was selected because of its location in the Himalayas, its heterogeneous physiographic diversity, and because it is one of the most rapidly expanding cities in South Asia [43,45]. Considerable LULC change has taken place as urban areas have expanded from the city center altering the dynamics of land–resource prospects and this analysis is essential for future sustainable development in the Kathmandu Valley.

We explore the following: (a) the spatiotemporal dynamics of LULC and their influence on intensification of SUHI for the Kathmandu Valley; (b) for each distinct LULC class we assess the mean LST along with their influence on the city landscape; (c) the relationship between LST and indicators of LULC (NDBI, NDVI, and NDMI) to aid assessment of the relative importance of impervious space and vegetation space with respect to LST in the Kathmandu Valley; and d) LST, NDVI, and NDBI on an urban–rural gradient from the city center to the urban periphery.

2. Materials and Methods

2.1. Study Area

Kathmandu Valley is situated in central Nepal and comprises 694.27 km², which includes three government districts: Bhaktapur, Lalitpur, and Kathmandu [47], and is located between 27°31′ to 27°50′ North, and between 85°11′ to 85°34′ East in the Himalaya mountains (Figure 1). It had a population of 1.11 million in 1991 which more than doubled to 2.52 million by 2016 [49]. Elevation ranges from 837 m to 2723 m, with the central part of the valley ranging from 1200 m to 1500 m. It has a dry-winter humid subtropical climate (Cwa) [50]. The valley has been created by the Bagmati River [43] which flows from north-east to south-west. Mean annual temperature and precipitation are 18.1 °C and 1407 mm, respectively [47].

2.2. Data

For our study, LULC data were based on images for individual days in April (early summer) at decadal intervals using Landsat TM and Landsat 8 OLI/TIRS at four different time-points (1988, 1998, 2008, 2018) which were retrieved from the USGS website (https://earthexplorer.usgs.gov) (

Table 1. Description of the Landsat datasets used for the study area, Kathmandu Valley.

Sensor	Path/Row	Resolution	Acquisition Date	Time (GMT)	Constants of Thermal Conversion		Source
					${\mathit{K}}_1$	${\mathit{K}}_2$
Landsat-5 TM	141/41	30 m	3 April 1988	04:18:24	607.76 (Band 6)	1260.56 (Band 6)	United States Geological Survey (USGS) web portal (https://earthexplorer.usgs.gov/)
			15 April 1998	04:25:19	607.76 (Band 6)	1260.56 (Band 6)
			26 April 2008	04:37:22	607.76 (Band 6)	1260.56 (Band 6)
Landsat-8 OLI/TIRS			22 April 2018	04:47:35	774.8853 (Band 10)	1321.0789 (Band 10)
ASTER	-		17 October 2011	-	-	-

). These data were also used to create LST based on their thermal bands (Band-6 of Landsat 5-TM, and Band-10 of Landsat 8-TIRS) and NDVI, NDBI, and NDMI based on optical and near-infrared bands (Band-3, 4, and 5 of Landsat 5-TM, and Band-4, 5, and 6 of Landsat 8-OLI), respectively. For the validation of weather conditions, we have explored NASA’s POWER project data on selected sample location for our study area and we have found that weather conditions of our selected time points are satisfactory (please refer to supplementary excel file).

Topographical data at the scale of 1:25,000 [51] and Google Earth images of the study area also were used. To overcome the atmospheric error in satellite images, they were pre-processed, including atmospheric and radiometric correction. Decadal intervals were chosen to determine growth in urban expansion and its influence in the transformation of LST, as it has been suggested as an optimal time period to measure LULC change, as well as SUHI development [10]. The Global Digital Elevation Model (GDEM) of ASTER was used to depict the topography of Kathmandu valley and to help understand its physiographic dynamics. We recognize the fact that just selecting single days to represent a decadal point in time does not address the problem of variation in LST between days or even between years but we consider our data to be a starting point for further analysis that may later address such issues.

2.3. Retrieval of LULC

We used freely available terrain corrected Landsat (Level 1T) datasets with best cloud-free data (less than 10% cloud cover) of UTM zone 45N. Our area of interest (AOI) had zero numbers of spiked digital numbers (DN) due to the availability of cloud free imageries and hence, we did not need to do masking and exclusion for correction of co-registration, cloud, cloud shadow, and gap-filling. Radiometric and geometric corrections are part of the pre-processing function of satellite images. In this study, we used ENVI software to carry out the atmospheric image correction process. Following this the DN values of images were converted into radiance values. All images were further verified for their accuracy and the root mean square (RMS) of the geometric rectification of less than 15 m (0.5) pixels was accepted. The Flash Line-of-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) model was used for atmospheric correction and LULC data were extracted using ENVI software [12,46].

There are multiple image classification algorithms, such as support vector machine (SVM), maximum likelihood classifier (MLC), and decision tree (DT) [52]. In this study we chose to use SVM for the classification of images as it is flexible, uses a non-parametric approach and is widely used for the extraction of LULC data bases [53]. Further, according to some authors, SVM has a higher level of accuracy than MLC [54]. Similarly, others have shown SVM performed better in the assessment of land cover changes and urban development than MLC and DT [55]. SVM can be grouped into four kernels function, such as linear, polynomial, radial function, and sigmoid [56]. The Radial Basic Function kernel was used at the time of extraction of different LULC classes as it usually provides better results than other machine learning [57]. Detailed field visits were conducted to ground-proof our analyses and five main LULC classes were identified based on the Anderson Classification Scheme at Level-I [58]: Impervious land (IL), Agriculture land (AL), Vegetation Land (VL), Open Land (BL), and Water Body (WB) (Table 1).

Accuracy assessment is essential for land cover data developed from remote sensing technology [59,60]. Overall accuracy (OA), User’s accuracy (UA), and producer’s accuracy (PA) assessment were assessed based on the field reference information. A topographical map of the study area for the scale of 1:25000 was used as developed by the Survey Department of Nepal in 1995 [51]. Similarly, Google Earth images for multiple dates were used in the assessment. OA represents what proportion of references sites were mapped correctly whereas UA is defined as the accuracy from the knowledge of a map user, not the map maker. PA is defined as the accuracy of a map from the knowledge of the map maker (the producer) [60].

Accuracy assessment reports were generated for each class of LULC using 200 stratified random samples points (total 1000 random points) and an error matrix for each time point was created to allow accuracy assessments [61]. UA, PA, OA, and Kappa coefficients were calculated based on the error matrix for each time point (1988–2018).UA, PA, OA, and Kappa coefficients were calculated based on error matrix of respective time points (1988–2018). To optimize classification accuracy: User, Producer, Overall accuracy, and Kappa coefficients were estimated using Equations (1)–(4) [10,61].

$$UserAccuracy=\left\{\frac{{{\displaystyle \sum }}^{\text{}}\in }{{{\displaystyle \sum }}^{\text{}}\vartheta \left(\omega \right)}\times 100\right\}$$

(1)

$$ProducerAccuracy=\left\{\frac{{{\displaystyle \sum }}^{\text{}}\in }{{{\displaystyle \sum }}^{\text{}}\vartheta \left(\phi \right)}\times 100\right\}$$

(2)

$$OverallAccuracy=\left\{\frac{{{\displaystyle \sum }}^{\text{}}\beta }{{{\displaystyle \sum }}^{\text{}}\gamma }\times 100\right\}$$

(3)

$$KappaCoefficients=\frac{N{{\displaystyle \sum }}_{i=1}^r{X}_{ii}-{{\displaystyle \sum }}_{i=1}^r(X_{i+}\times X_{+i})}{N^2-{{\displaystyle \sum }}_{i=1}^r(X_{i+}\times X_{+i})}$$

(4)

where, $\in$ defines corrected classified pixels (CCP) (category); $\vartheta \left(\omega \right)$ defines classified pixels (CP) in that category (row total (RT)); $\vartheta \left(\phi \right)$ defines CP in that category (column total (CT)); $\beta$ defines CCP (diagonal); $\gamma$ defines classified reference pixels in that category; $N$ defines total samples;$r$ defines number of rows error matrix (EM); $X_{ii}$ defines total corrected samples in $i^{th}$ row and column:; $X_{+i}$ defines RT: $X_{i+}$ defines CT.

2.4. Land Surface Temperature (LST) Retrieval

2.4.1. Land Surface Emissivity

Land surface emissivity (LSE) is an important parameter to calculate LST [62,63]. NDVI threshold (NDVI^THR) method was used to calculate LSE because it can differentiate pixels of vegetation, water, and soil significantly [64]. It was calculated using Equation (5).

$${\epsilon }_{sv}={\epsilon }_v{P}_V+{\epsilon }_s\left(1-P_V\right)+C$$

(5)

where, ${\epsilon }_{sv}$ is emissivity of soil and vegetation; ${\epsilon }_v$ is emissivity of vegetation; ${\epsilon }_s$ is emissivity of soil; $P_V$ is proportion of vegetation (using Equation (6)); $C$ is defining constant for surface characteristics (using Equation (7)).

$$P_V={\left[\frac{NDVI-NDV{I}_s}{NDV{I}_v-NDV{I}_s}\right]}^2$$

(6)

where, $NDVI$ is estimated using Equation (19) in Section 2.5; $NDV{I}_s$ is NDVI of pure soil; $NDV{I}_v$ is NDVI of pure vegetation.

$$C=\left(1-{\epsilon }_s\right){\epsilon }_{v}F\left(1-P_V\right)$$

(7)

where, $F$ is geometric factor (it depends on surface geometry, commonly considered as 0.55 [64,65]).

Skokovic et al. [66], Sorbino et al. [64], Sekertekin and Bonafoni [65] calculated LSE using NDVI^THR values considering three distinct cases shown in Equation (8). The first case had NDVI values < 0.2, for bare soils; the second case had NDVI values, 0.2 ≤ NDVI ≥ 0.5 for mixed bare soils and vegetation; and the third case had NDVI values > 0.5 for fully vegetated areas.

$$\epsilon =\{\begin{array}{c}{a}_i{\rho }_R+b_i\\ {\epsilon }_v+{\epsilon }_s\left(1-P_V\right)+C\\ {\epsilon }_v+C\end{array}$$

(8)

where, $\epsilon$ is LSE; ${\rho }_R$ is the reflectance value of the red band, a_i and b_i are calculated using the empirical relationship for reflectance of red band and Moderate Resolution Imaging Spectroradiometer (MODIS) emissivity library [66].

We used the defined constant values of $\epsilon$ from Sekertekin and Bonafoni [65] for Landsat 5 (TM) and Landsat 8 (OLI/TIRS) in given Equation (9) and Equation (10), respectively.

$$\epsilon =\{\begin{array}{c}0.979+0.035{\rho }_R\to NDVI<0.2\\ 0.004P_V+0.986{\rho }_R\to 0.2\le NDVI\ge 0.5\\ 0.99\to NDVI>0.5\end{array}$$

(9)

$$\epsilon =\{\begin{array}{c}0.979+0.046{\rho }_R\to NDVI<0.2\\ 0.989P_V+0.977{\rho }_R\to 0.2\le NDVI\ge 0.5\\ 0.987+C\to NDVI>0.5\end{array}$$

(10)

2.4.2. LST

To estimate LST, we used the Radiative Transfer Equation (RTE) method [65] as follows:

$$B_i\left(T_i\right)=\frac{L_{max}-L_{min}}{QCA{L}_{max}-QCA{L}_{min}}\times \left(QCAL-QCA{L}_{min}\right)-L_{min}$$

(11)

where, $B_i\left(T_i\right)$ is the spectral radiance of the top of atmosphere (TOA) (watts/(m² $\times$ sr $\times$ µm)), $QCAL$ is DN, $QCA{L}_{min}$ and $QCA{L}_{max}$ defines the minimum and maximum DN values of the images, respectively; $L_{min}$ and $L_{max}$ are spectral radiance of TIR band at $QCA{L}_{min}$ and $QCA{L}_{max}$ respectively: these rescaling factor values can be found in the metadata of Landsat images. Equation (11) was used for Landsat 5 (TM) for Band-6.

$$B_i\left(T_i\right)=M_L\times Q_{cal}+A_L$$

(12)

where, $M_L$ is a multiplicative rescaling factor in the specific band from the metadata, $Q_{cal}$ is the quantized and calibrated DN values of standard product, and $A_L$ is additive rescaling factor in the specific band from the metadata. Equation (12) was used for Landsat 8 (OLI/TIRS) for Band-10.

We incorporated the RTE method for estimating LST as in Equation (13).

$$B_i\left(T_i\right)=\left[\epsilon B_{\lambda }\left(T_s\right)+\left(1-\epsilon \right)L_d^{\downarrow }\right]\tau +L_u^{\uparrow }$$

(13)

where, $B_i\left(T_i\right)$ is the spectral radiance of the top of atmosphere (TOA) (watts/(m² $\times$ sr $\times$ µm)) for band i in which have $T_i$ i.e., at-satellite brightness temperature; $\tau$ is atmospheric transmittance; $L_d^{\downarrow }$ is downwelling radiance; $L_u^{\uparrow }$ is upwelling radiance; $T_s$ is LST; $\epsilon$ is emissivity of band i. We obtained atmospheric values such as $\tau$ $L_d^{\downarrow }$,and $L_u^{\uparrow }$ using an online calculator tool, called, ‘Atmospheric Correction Parameter Calculator (ACPC)’ based on the given radiative transfer code of MODTRAN (http://atmcorr.gsfc.nasa.gov). $B_{\lambda }$ is Blackbody radiance at a temperature of $T_s$ as calculated in Equation (14) using the inversion of Equation (13).

$$B_{\lambda }\left(T_s\right)=\frac{B_i\left(T_i\right)-L_u^{\uparrow }-\tau \left(1-\epsilon \right)L_d^{\downarrow }}{\tau \epsilon }$$

(14)

$T_s$ in Kelvin was calculated using Equation (15):

$$T_s=\frac{K_2}{\ln \left(\frac{K_1}{\frac{B_i\left(T_i\right)-L_u^{\uparrow }-\tau \left(1-\epsilon \right)L_d^{\downarrow }}{\tau \epsilon }}+1\right)}$$

(15)

where, K₁ and K₂ were obtained from metadata file of respective time point images mentioned in Table 1.

2.4.3. LST Calculation (before Year 2000)

Because of the limitation of Atmospheric parameters values before the year 2000 in ACPC, we used a conventional method to extract LST values for Landsat 5 (TM) before the year 2000 only. Brightness temperature at sensor value was estimated using Equation (16) [10,13,32,41,67]:

$$B_i{T}_b=\left[\frac{K_2}{\ln \left(\frac{K_1}{B_i{T}_i}+1\right)}\right]$$

(16)

where, $B_i{T}_b$ is brightness temperature (At sensor) in Kelvin; $K_1$ and $K_2$ are the thermal conversion constants from the metadata (Landsat 5 TM (Band 6), and Landsat 8 OLI (Band 10)) (Table 1).

Using Equation (17), the derived LST (in Kelvin (K)) through correction of emissivity was assessed from the brightness temperature [10,32,41]:

$$T_s=\left[\raisebox{1ex}{$B_i{T}_b$}\!\left/ \!\raisebox{-1ex}{$1+w\left(\frac{B_i{T}_b}{p}\right)\ln \left(\epsilon \right)$}\right.\right]$$

(17)

where, $T_s$ is temperature (At sensor) in Kelvin; $\mathrm{w}$ is wavelength of emitted radiance (10.8 µm for Band 10 in Landsat 8 OLI and 11.5 µm for Band 6 in Landsat 5 TM); $\mathrm{p}=\mathrm{h}\times \mathrm{c}/\mathrm{s}\left(1.438\times {10}^{-2}\mathrm{mK}\right)$, $\mathrm{h}$ is Plank’s constant ($6.626\times {10}^{-34}\mathrm{Js}$), $s$ is the Boltzmann Constant ($1.38\times {10}^{-23}\mathrm{J}/\mathrm{K}$), and $\mathrm{c}$ is the velocity of light ($2.988\times {10}^8\mathrm{m}/\mathrm{s}$); $\epsilon$ is LSE.

Finally, we converted the $T_s$ (LST (Kelvin)) value into Celsius (°C) using Equation (18) [10,13,32]:

$$LST\left(°C\right)=T_s-273.15$$

(18)

2.5. NDVI

An important indicator of urban climate is the Normal Difference Vegetation Index (NDVI) [10,68]. It varies between −1 and +1, where large negative (and adjacent to negative) values, positive values, and low positive values indicate water bodies, vegetation, built-up areas, or bare soils, respectively [10]. It describes information about the amount of vegetation, and its phenology and health [38]. NDVI was estimated using Equation (19) using data from the red and near infrared (NIR) bands [19,20]:

$$NDVI=\left[\frac{NI{R}_{Band}-Re{d}_{Band}}{NI{R}_{Band}+Re{d}_{Band}}\right]$$

(19)

where, in Landsat 5 TM: NIR band denotes Band 4 (0.76–0.90 μm (wavelength)) and Red band denotes Band 3 (0.63–0.69 μm (wavelength)); while in Landsat 8 OLI: NIR band denotes Band 5 (0.85–0.88 μm (wavelength)) and Red band denotes Band 4 (0.64–0.67 μm (wavelength)).

2.6. NDBI

Another important urban climate indicator is the Normal Difference Built-up Index (NDBI) [10,68]. It varies from −1 to +1, where negative values specify water bodies and vegetation, positive value indicates the built-up area, and low positive value specifies types of bare soils [48]. It represents information about the imperviousness of the landscape [38] and was estimated using Equation (20) with NIR band and mid infrared (MIR) [19,20]:

$$NDBI=\left[\frac{MI{R}_{Band}-NI{R}_{Band}}{MI{R}_{Band}+NI{R}_{Band}}\right]$$

(20)

where, in Landsat 5 TM: MIR band denotes Band 5 (1.55–1.75 μm (wavelength)) and NIR band denotes Band 4 (0.76–0.90 μm (wavelength)); while in Landsat 8 OLI: MIR band denotes Band 6 (1.57–1.65 μm (wavelength)) and NIR band denotes Band 5 (0.85–0.88 μm (wavelength)).

2.7. NDMI

Another significant urban climate indicator is the Normal Difference Moisture Index (NDMI) [39]. It also varies from −1 to +1, where a positive value specifies water bodies and vegetation and a negative value indicates built-up areas and bare soils. It signifies information about the moisture in the landscape and was estimated using Equation (21) with NIR band and mid infrared (MIR) [39]:

$$NDMI=\left[\frac{NI{R}_{Band}-MI{R}_{Band}}{NI{R}_{Band}+MI{R}_{Band}}\right]$$

(21)

where, in Landsat 5 TM: MIR band denotes Band 5 (1.55–1.75 μm (wavelength)) and NIR band denotes Band 4 (0.76–0.90 μm (wavelength)); while in Landsat 8 OLI: MIR band denotes Band 6 (1.57–1.65 μm (wavelength)) and NIR band denotes Band 5 (0.85–0.88 μm (wavelength)).

2.8. Analysis of Urban–Rural Gradient

The gradient approach is frequently used to evaluate spatiotemporal differences in the environment [10,48]. Here we used it to assess the spatial dynamics of mean LST, NDBI, and NDVI at 1 km intervals from the center of the city to the periphery/suburban/rural area of the city up to a maximum of 26 km (see Figure 1d for the transect of 1 km of buffers up to 26 km) [10].

2.9. Statistical Analysis (Pearson’s Correlation Coefficient)

To visualize the effects of the environmental variables (NDVI/NDBI/NDMI) on LST intensification, scatter plots were made for all time-points, (1988, 1998, 2008, and 2018) using linear regression. LST, NDBI, NDVI, and NDMI pixels were transformed into point data [10,68]. Pearson’s correlation coefficient (‘r’) was used to measure the relationship among LST v/s NDVI, LST v/s NDBI, and LST v/s NDMI, where LST was the dependent variable, and NDVI/NDBI/NDMI were the independent variables. Pearson’s ‘r’ was calculated through Equation (22):

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

(22)

where, r is Pearson’s correlation coefficients, $x$ represents NDVI/NDBI/NDMI measuring value of $x_i$, $y$ represents LST measuring value of $y_i$. $x_i$ and $y_i$ are single sample indexed with $i$. $\overline{x}$ and $\overline{y}$ defining the single samples indexed of $x_i$ and $y_i$, respectively.

3. Results

3.1. Accuracy Assessment of LULC Classification

LULC classification has been performed by SVM method. Accuracy assessments has been done by random sampling method by 200 sample points for each LULC classes and then error matrix has been made for each time points (i.e., 1998, 1998, 2008, and 2018). User and producer accuracy were greater than 80%, an overall accuracy was greater than 90%. Kappa coefficients were 0.90 in 1988, 0.92 in 1998, 0.94 in 2008, and 0.96 in 2018 (Appendix A

Table A1. Accuracy assessment report of Land Use/Land Cover classification in Kathmandu (1988–2018).

	LULC Class	Year
		1988	1998	2008	2018
User Accuracy (%)	IL	88.0	89.0	91.0	96.0
	AL	95.5	96.0	96.5	96.0
	VL	84.5	90.0	95.5	96.5
	OL	93.5	94.0	95.0	96.0
	WB	97.5	98.0	96.5	98.5
Producer Accuracy (%)	IL	92.6	87.3	94.8	93.7
	AL	83.8	86.9	91.0	93.7
	VL	90.0	95.7	96.4	97.8
	OL	99.4	99.5	99.5	100
	WB	94.7	97.0	98.0	98.5
Overall Accuracy (%)		95.2	93.4	95.2	96.6
Kappa Coefficient		0.897	0.917	0.940	0.958

).

3.2. LULC Analysis

Spatiotemporal analyses of LULC of Kathmandu Valley (Figure 2 and Figure 3) showed that AL, VL, and IL were the major land use classes between 1988 and 2018, and that there was large increase in IL almost entirely at the expense of AL. Over this period, IL expanded by 113.44 km² (16.34% of total area) whereas AL reduced by 118.29 km² (17.03% of total area). Small changes in land were observed for OL and WB. Summaries of changes in LULC for each class are provided in

Table A2. The summary Land Use/Land Cover statistics of Kathmandu at 10-years interval from 1988 to 2018.

Year		LULC Classes
		IL	AL	VL	OL	WB	Total
1988	Km2	37.98	418.32	230.69	1.97	5.31	694.27
	%	5.47	60.25	33.23	0.28	0.76	100
1998	Km2	58.96	391.46	235.85	3.45	4.55	694.27
	%	8.49	56.38	33.97	0.5	0.66	100
2008	Km2	95.45	356.63	235.4	3.34	3.45	694.27
	%	13.75	51.37	33.91	0.48	0.5	100
2018	Km2	151.42	300.03	236.76	3.05	3.01	694.27
	%	21.81	43.22	34.1	0.44	0.43	100

, and gains and losses for each class in

Table A3. Summary Land Use/Land Cover change statistics for Kathmandu Valley.

Year		LULC Type
		IL	AL	VL	OL	WB
1988–1998	Km2	20.99	−26.86	5.16	1.48	−0.76
	%	3.02	−3.87	0.74	0.22	−0.1
1998–2008	Km2	36.49	−34.83	−0.47	−0.11	−1.1
	%	5.26	−5.01	−0.06	−0.02	−0.16
2008–2018	Km2	55.97	−56.6	1.38	−0.3	−0.44
	%	8.06	−8.15	0.19	−0.04	−0.07
1988–2018	Km2	113.44	−118.29	6.07	1.07	−2.3
	%	16.34	−17.03	0.87	0.16	−0.33

and Figure A1). A small increase was observed for VL class during 1988–2018, increasing from 230.69 km² in 1988 to 236.76 km² in 2018 with a total conversion of 6.07 km² (0.87% of total area) from other LULC classes.

3.3. LST Analysis

Mean LST (

Table 2. Maximum, minimum, and mean LST (°C) values in Kathmandu (1988–2018).

Date	Minimum (°C)	Maximum (°C)	Mean (°C)	Standard Deviation
3 April 1988	11.88	35.25	23.08	3.09
15 April 1998	17.93	39.15	25.46	3.20
26 April 2008	17.47	40.30	28.86	3.31
22 April 2018	16.96	39.46	28.35	3.40
The Mean LST difference in different time periods (°C)
1988–1998	1998–2008	2008–2018		1988–2018
2.38	3.4	−0.51		5.27

) increased from 23.08 °C in 1988 to 28.36 °C in 2008 and remained at that level in 2018. The spatial dynamics of changes over time in LST for the Kathmandu Valley are shown in Figure 4. The greatest transformation in LST between 1988 and 2018 was in the central, southern, and top north-eastern parts of the valley due to the massive intensification in LST of AL and IL, whereas the lowest LST was observed at the top of the northern, south-western (most bottom), and western margins of the valley because of the vast VL presence (see Figure 2). In 1988, LST values for the suburban/rural areas in Kathmandu Valley such as Katunje, Bhaktapur city, Gapali, Liwali (eastern side of the valley), Chandragiri (western side of the valley) (Figure 2 and Figure 4), and Godavari (southern side of the valley) were in the range of 15–25 °C, whereas in the center of the valley, Kathmandu City, LST was 20–28 °C (Figure 4), giving a difference of 3–13 °C between suburban/rural areas and urban areas (Figure 4). In 1998, a difference of 5–15 °C LST was experienced across the city area, with urban areas experiencing LST values of 23–33 °C, some 18–23 °C higher than that of suburban/rural areas. In 2008, a difference of approximately 7–10 °C LST was experienced with LST in the range of 25–40 °C for urban areas, indicating a difference of 18–30 °C (Figure 4). In 2018, the city center observed a range of 25–39 °C, some 7–9 °C higher than that of the suburban/rural area (range of 18–30 °C), largely because Bhaktapur and Chandragiri had developed into dense urban areas (Figure 2).

3.4. LULC Differences in LST

The greatest changes in LST were observed for IL (

Table 3. Mean LST for different of LULC classes in Kathmandu Valley during 1988–2018.

LULC Class	Mean LST (°C)				The Difference of Mean LST (°C)
	1988	1998	2008	2018	1988–1998	1998–2008	2008–2018	1988–2018
Impervious Land	23.81	23.63	30.97	30.92	−0.18	7.34	−0.05	7.11
Agriculture Land	23.69	23.80	30.38	29.67	0.11	6.58	−0.71	5.98
Vegetation Land	21.80	22.03	25.65	24.98	0.23	3.62	−0.67	3.18
Open Land	24.89	23.89	31.22	30.36	−1	7.33	−0.86	5.47
Water Body	22.62	22.46	27.55	26.81	−0.16	5.09	−0.74	4.19

). Mean LST for IL increased from 23.8 °C in 1988 to 23.6 °C in 1998 and again to 31.0 °C in 2008 and to 30.9 °C in 2018, indicating a mean increase of 7.11 °C of over the whole time period but noting this increase was not even over the time period (Table 3 and Figure 5). Similar, but slightly lower results were observed for AL, (23.7 °C, 30.4 °C, 23.8 °C, and 29.7 °C for the same time periods showing an overall increase of 6.0 °C (Table 3 and Figure 5). Changes for VL over the whole time period of 3.2 °C was much less than that for other classes, as it was 21.8 °C in 1988, 22.0 °C in 1998, 25.7 °C in 2008, 25.0 °C in 2018 (Table 3 and Figure 5). For OL and WB classes the overall change was an increase of 5.5 °C and 4.2 °C, respectively (details in Table 3).

Mean LST for IL was greater than that for VL by 2.1–5.9 °C and for WB by 1.2–3.4 °C at all time-points (Figure 6 and

Table A4. Comparisons of mean LST between Land Use/Land Cover classes from 1988 to 2018.

Comparison of LULC Class	LULC Class (Cross Cover Comparison)	Mean LST Magnitude (°C)
		1988	1998	2008	2018
IL vs. Other Class	IL-VL	2.01	1.6	5.32	5.94
	IL-WB	1.19	1.17	3.42	4.11
	IL-AL	0.12	−0.17	0.59	1.25
	IL-OL	−1.08	−0.26	−0.25	0.56
VL vs. Other Class	VL-IL	−2.01	−1.6	−5.32	−5.94
	VL-AL	−1.89	−1.77	−4.73	−4.69
	VL-OL	−3.09	−1.86	−5.57	−5.38
	VL-WB	−0.82	−0.43	−1.9	−1.83
WB vs. Other Class	WB-IL	−1.19	−1.17	−3.42	−4.11
	WB-AL	−1.07	−1.34	−2.83	−2.86
	WB-VL	0.82	0.43	1.9	1.83
	WB-OL	−2.27	−1.43	−3.67	−1.55

). However, mean LST was lower for IL than OL by 0.3–1.1 °C (except in 2018, where IL had a higher mean LST than OL by 0.6 °C). IL had a greater LST mean than AL of 0.1–1.3 °C at all consecutive time-points (except in 1998 see Table A4). It is evident that mean LST for VL is lower than that for all other classes of LULC on all dates. Similarly, it was evident that IL had higher mean LST than VL and WB because of the presence of the vegetation and water surface, respectively [13,69].

3.5. Spatiotemporal NDVI, NDBI, and NDMI Patterns and Their Influence on LST

NDVI mean was 0.37 in 1988, increasing to 0.38 in 1998, declining to 0.36 in 2008, but increasing again to 0.43 in 2018 (

Table 4. Statistics of Land Indices and their correlation with LST.

Date	Statistics of NDVI and NDBI					Correlation with LST
		Minimum	Maximum	Mean	Standard Deviation	Correlation Coefficient	Significance (p)
3 April 1988	NDVI	−0.073	0.704	0.365	0.111	−0.6665	p < 0.001
	NDBI	−0.549	0.313	−0.017	0.136	0.7780	p < 0.001
	NDMI	−0.313	0.549	0.017	0.136	−0.7780	p < 0.001
15 April 1998	NDVI	−0.074	0.704	0.367	0.111	−0.6527	p < 0.001
	NDBI	−0.671	0.416	−0.095	0.144	0.7783	p < 0.001
	NDMI	−0.416	0.671	0.095	0.144	−0.7783	p < 0.001
26 April 2008	NDVI	−0.042	0.725	0.356	0.135	−0.7678	p < 0.001
	NDBI	−0.549	0.318	−0.058	0.141	0.8414	p < 0.001
	NDMI	−0.318	0.549	0.058	0.141	−0.8414	p < 0.001
22 April 2018	NDVI	−0.025	0.787	0.430	0.168	−0.8045	p < 0.001
	NDBI	−0.519	0.378	−0.106	0.150	0.8264	p < 0.001
	NDMI	−0.378	0.519	0.106	0.150	−0.8264	p < 0.001

). The raised mean NDVI in 2018 was most likely due to the expansion of vegetation (trees and grasses) on the northern side (where Shivpuri Nagarjun National Park is located over part of Budhanilakantha, Gokarneshwor, Tokha, and Shankharapur), the western side (Nagarjun Forest Reserve), the southern side (Godavari Forest), the north-eastern side (open forest which has increased in Shankharapura and Changunarayan), and the south-western side (open forested areas increasing over Dakshinkali, Kirtipur, and Chandragiri). Consequently, the highest NDVI was observed in the same northern, southern, and western sides of the Valley for each date due to the presence of dense vegetation in national parks and forest reserves. However, by 2018 the area of vegetation also had increased in the central part of the valley in areas such as Gaucharan, Pashupatinath, Bhandarkhal Jungle, Panipokhari, Samakhusi, and Amideva Budha Park. The lowest NDVI was in some central and eastern parts of the Valley, perhaps due to the large areas of AL (bare land type) and IL, respectively (Figure 7).

The correlations between LST and NDVI were significantly negative for all time points (p < 0.001), declining from −0.67 in 1988, to −0.80 in 2018, but increasing in 1998 to −0.65 (Table 4 and Figure A2). It was very clear that reduction in vegetation cover resulted in higher LST and vice versa (Figure 4).

Mean NDBI also declined from −0.02 in 1988, to 0.11 in 2018 with an increase in 2008 to −0.06 (Table 4). The highest values for NDBI were in the central part in the Valley at each time-point because of densely built-up areas of Kathmandu City, Latlipur, Bhaktapur, and Madhyapur Thimi, whereas the lowest NDBI values were in some central and eastern parts of the Valley, such as Shivpuri Nagarjun National Park, Nagarjun Forest Reserve, Godavari Forest, and Bhandarkhal Jungle because of the high density VL (Figure 8).

A significantly positive correlation was observed between LST and NDBI (Figure A3), where the coefficient of determination was greater than 0.60 at each time-point with 0.61, 0.61, 0.71, and 0.68 in 1988, 1998, 2008, and 2018 (all p < 0.001), respectively, indicating that there is strong relationship between LST and NDBI (Table 4, Figure 4). The significant correlation between LST and NDBI (Figure A2) varied over time: 0.78, in 1988, 0.78 in 1998, and 0.84 in 2008, and 0.83 (all p < 0.001) in 2018 (Table 4 and Figure A3). The positive correlation coefficients suggested that with increased NDBI (Figure 8), LST also increased at all time-points (Figure 4). This was particularly so in 2008 when mean LST was 3.4 °C higher than 1998 (Table 3) because the distribution of high NDBI increased due to the rapid development of IL (Figure 2 and Table A4. It is clear that less built-up areas resulted in lower LST and vice versa (Figure 2 and Figure 4).

Mean NDMI also increased from 0.062 in 1988 to 0.101 in 2018 but an interesting drop of mean NDMI value of 0.059 was seen in 2008 (Table 4). Additionally, in 2008, mean LST increased and a greater mean NDMI was estimated in 2018 than the mean NDMI of 2008 and this may be the reason mean LST was lower in 2018 than 2008. The lowest values for NDMI were seen over agricultural land and impervious land in the Valley at each time-point due to the densely built-up areas of Kathmandu City, Latlipur, Bhaktapur, and Madhyapur Thimi, whereas the highest NDMI values were found for vegetated areas and water bodies in areas of the Kathmandu valley such as Nagarjun Forest Reserve, Shivpuri Nagarjun National Park, Godavari Forest, and Bhandarkhal Jungle, and Bhagmati river, respectively (Figure 9).

The correlations between LST and NDMI were also significantly negative for all time points (p < 0.001), declining from −0.78 in 1988, to −0.83 in 2018, but most decreased in 2008 to −0.84 (Table 4 and Figure A4). It was very apparent that reduction in moisture content resulted in higher LST and vice versa (Figure 4).

3.6. Analysis of Urban–Rural Gradient Pattern

Changes in LST, NDBI, NDVI, and NDMI across the city/rural gradient for at each time-point (i.e., 1988, 1998, 2008, and 2018) are shown in Figure 10 (see Figure 1d for demarcation of buffer zones). Gradient patterns are reasonably similar across the years for both NDVI and NDBI with the latter mirroring the former. However, there was much greater variation in LST between years. In 1988 and 1998 mean NDVI increased gradually with increasing distance from the city center (roughly 23–24 °C) up to around 12 km, remaining at roughly 26–27 °C until about 17 km when it gradually declined to 25 °C. In 2008 and 2018, mean NDVI was initially higher, rising more rapidly with increasing distance from 28–29 °C to around 32 °C at 17 km before declining gradually to 30–31 °C at the city periphery. Mean NDBI mirrored NDVI. In 1988 and 1998 it declined from around 23–24 °C to around 21–22 °C between 12 and 17 km then rose to around 23 °C at 22 km before declining again. In 2008 and 2018 mean NDBI at the city center was higher at 27–28 °C and declined to 25.5 °C at around 17 km before increasing to 26–27 °C at 22 km followed by a slight decline. Mean LST in the city center appeared to increase with each decadal transect. In 1988 it dropped from 25.5 °C to around 22.5 °C at 17 km and then showed a sharp increase to 26.3 °C at about 22 km before then sharply decreasing at the city boundary because of greater agricultural activities between 17–24 km (See top square box at North-Eastern side in Figure 2a–d, Figure 4a–d and Figure 9a–d). In 1998 mean LST at the city center was higher (around 28.8 °C) and declined 24.3 °C at 15 km before rising again to 28.8 °C at 22 km before then again sharply declining. In 2008 at the center of city, mean LST was around 31.4 °C, declining to 27.8 °C at 17 km then rising sharply to 30.2 °C at 22 km and again falling to around 28 °C at 25 km due to agriculture activities, and NDMI (<0 means the lowest moisture condition) also showed lowest moisture there (Figure 9a–d). In 2018 mean LST was highest at around 32 °C within 1 km of the city center rapidly declining to 25 °C at 17 km, and as in previous decades, a further spike of 30–31 °C at around 22 km before declining again.

Consequently, it is clear that mean LST mean was dramatically greater from the city center to 15 kms compared to the surrounding area to the city periphery. Mean NDBI value was also higher in the city center to 15 km compared to the surrounding area (Figure 10), while the mean NDVI trend was the opposite, as NDVI value is lower in the area of city center up to 15 km than in the surrounding area for each decadal time-point (Figure 10). It was also very evident that mean NDVI mean was low in the areas between city center and 10 km and high between 12 km 26 km (Figure 10). Therefore, it seems clear that the highly developed VL resulted in low LST and highly developed IL resulted in higher LST and vice versa for both. Interestingly, in the 20–25 kms area, LST was again higher due to increased agriculture (Figure 2 and Figure 10). These patterns, we believe, reflect the changes in development of built areas in the innermost parts of the city and the higher levels of vegetation towards the city periphery.

4. Discussion

We have shown that in the 30 year period between 1988 and 2018 there was a massive change in land use in the Kathmandu valley with extensive growth of IL at the cost of AL and to a lesser extent VL due to intense urbanization. At the same time, LST and NDBI have increased considerably while NDVI has reduced in core urban areas. Here we discuss the implications of our results in the context of previous studies.

4.1. Urbanization from the Perspective of LULC Transformation and Population Explosion

The capital of Nepal, Kathmandu City, is one of the fastest growing cities in South-East Asia [46,47], located in a valley with high mountains on all sides. It is only the central part of the valley which is urbanizing fast [46,47] and this rapid urbanization is due to religious and tourism attractions and the location of the capital administration. The rapid growth of built up areas in the form of transport networks, residential, commercial, industrial buildings, and associated parking lots have resulted in losses of some LULC classes, particularly AL, as we have shown here. This widespread transformation in LULC has led to elevated LST across a large part of the landscape [10,38,47,70]. The rapid urban expansion is the result of the three-fold population increase over the last three decades in Kathmandu City (Figure A5) since it was only 0.35 million in 1988 but has increased to 1.33 million by 2018. Similar rapid population growth has occurred across the whole of the Kathmandu Valley, including the districts of Bhaktapur, Kathmandu, and Lalitpur, where the total population was 0.77 million in 1981 but had increased to 2.88 million by 2016 [49] (Figure A6). Similar levels of rapid population growth and associated urban expansion have been observed in many other cities around the world including Chennai in India [71]; Lucknow in India [3,72]; Agra in India [73]; Xuzhou in China [74]; Baguio in the Philippines [75]; Kandy City in Sri Lanka [48]; Tokyo in Japan [11]; Tehran in Iran [10]; Istanbul in Turkey [76]; Sobotka in Poland [77]; Santiago de Chile [78]; Mekelle in Ethiopia [79]; Bucharest city in Romania, Budapest city in Hungary, Prague city in Czech Republic, Sofia city in Bulgaria, and Warsaw city in Poland [80]; and São José dos Campos in Brazil [81].

We have shown that IL in the last three decades, especially most recently in 2008 and 2018, has swiftly extended towards the eastern side (e.g., Madhyapur Thimi, Katunje, Bhaktapur, Liwali, Gapali, Jhaukhel, and Duwakot), the western side (e.g., Chandragiri, Kirtipur, Sitapalia, Gairi Gaun, Boharatar, and Ramkot), the northern side (e.g., Kapan, Dharampur, Bhangal, Thapa Gaun, Mahankal, Gongabu, Paiyatar, and Hiledol), and the southern side (e.g., Patan, Bhaisepati, Kusunti, Mitra Tole, Sunakothi, Dhapakhel, Thaiba, Thecho, and Bullu) (Figure 2). Such growth is predicted to continue as the United Nations has projected that Kathmandu City population will grow to 2.2 million by 2035 [4], roughly doubling the 2018 population (Figure A5), whereas NPCS has predicted the population for the whole of the Kathmandu Valley to increase to 3.85 million by 2031 (Figure A6).

4.2. Phenomena of SUHI and Sustainable Planning

Our results suggest that the rapid expansion of built-up areas and their influence in increasing SUHI is occurring at a large scale in Kathmandu Valley. This is largely due to the conversion of AL and, to a lesser extent, VL into IL in the form of transport networks, industrial, commercial, residential, parking lots, and other paved surfaces. We found that both VL and WB had lower mean LST than IL, AL, and OL at all studied sequential time-points. This contrasts with the results of other studies which indicated that IL had higher LST than others classes of LULC especially VL and WB, in tropical montane cities, such as, Tehran in Iran [10], Kandy in Sri Lanka [48], and Baguio in Philippines [75]. We noted that LST was higher in 2008 and we recognize the limitation in our study that in selecting single dates for our Landsat images we cannot account for daily variations in LST that are naturally likely to occur. However, overall mean LST increased over the whole 40-year time period. However, there could be other factors (like, elevation) influencing LST intensification which could be examined in the future.

We found that LST intensified throughout the central city area of the Kathmandu Valley as the IL rapidly increased. In many western cities, attempts to reduce the effects of SUHI have included strategies such as introducing more street trees or growing plants on roofs, so-called green roofs, as well as developing materials that cool roofs and cool pavements, and the use of light materials [82]. Other cooling mechanisms included improving wind flow, by carefully designing the size, shape, and orientation of buildings [83]. New cities such as Kathmandu which have grown with little planning need to recognize the value of such greening designs with the opportunity to enhance environmental sustainability in the city [10,18,41,68].

We also found that the mean LST was always higher in the city center than its periphery and that in general LST has risen over the 40-year period. In between 12–18 km from the city center, LST level increased due to increased NDBI because of greater population pressure resulting in urban expansion. Kathmandu Valley has also observed an increase in SUHI from the city center towards the city periphery. The reason behind this intensification was due to dense urban development (especially in the areas of Kaldhara, Chhetrapati, Sanepa, Kalimati, Siddhitol, Bhimsengola, Sinamangal, Kumaripati, Ekantakuna, and Narephate) and at the periphery of the city with the loss of green spaces (especially in the areas of Budhanilkantha, Tarakeswar, Sitapalia, Charghare, Bhatkepati, Suyel Gaun, Dadhikot, Chhaling, and Taudol).

4.3. Urban Sustainability Implication

Typically, in the relatively unplanned growth of built-up areas micro- to macro- level features develop to reflect economic and social needs [10,32,45,84]. Greater employment opportunities and higher economic development are important reasons why people are drawn to urban areas but at the same time the resulting reduction in AL and VL leads to intensification of SUHI and consequent harm to the local environment and the urban population [10,85,86].

With so few open areas available for development in the city center we found that Kathmandu valley has had higher IL development at the periphery than in the city center between 1988 and 2018. As a result, from 6–8 kms, IL expanded largely at the cost of AL. This partly explains some of the observed changes in LULC statistics (Figure 2 and Table A2).

Our results suggest that SUHI has intensified at all-time points resulting from the increase in IL at the cost of VL and AL which is likely to have severe consequences for the local environment. Sustainable land management practices can decrease the negative effects of stressors like climate change [7]. Therefore, to minimize some of the severe effects of SUHI novel greening strategies need to be developed, as discussed above, as well as reducing runoff and enhancing availability of freshwater through the creation of ponds/lakes, and use of rain water harvesting. Such structures can enhance the resilience of local ecosystems [9,10,18]. Increased public awareness regarding the severe effects of SUHI phenomena may result in pressure at the local and national level on public or private authorities to reduce the SUHI phenomena through adoption of remedies to minimize its effects as well as improving future sustainable planning for new urban areas [10].

5. Conclusions

We have shown a significant negative and positive relationship between LST and NDVI, and LST and NDBI, respectively. We interpret this to show that this indicates that vegetation had a very significant role in decreasing LST by nearly 0.82–5.94 °C compared to other LULC classes over the last 30 years in the study area. At the same time, we are concerned about the increase in mean LST of 0.12–5.94 °C for built-up areas. Our results show that there was a sharp decrease of mean LST with NDBI from city center to 15 kms with NDVI showing the opposite pattern sharply increasing from the city center to 15 kms, clearly demonstrating the creation of SUHI. Further, at around 20–25 kms from the city center, the mean LST again rose due to the intensive agriculture there as we have found mean NDVI and as well as mean NDMI also declined there, which reveals agriculture activities’ effects on higher LST distribution. It was very apparent that because of explosive increase in IL there was a resulting loss in AL and VL. One noteworthy observation is that OL had greater mean LST than all other classes of LULC at all consecutive time points. We conclude that it is essential to measure thermal state for cities over time to depict LULC and SUHI creation because it can explain the consequences of historical changes in the city’s landscape. The patterns of LULC and LST observed in this study we hope will be useful for future urban planning and policy making in the Kathmandu Valley.

The enormous changes in urbanization resulting from rapidly increasing population growth in the Kathmandu valley and their consequences for LULC and LST change are of great concern particularly as population growth is predicted to continue. To avoid severe consequences of SUHI, strong planning policies and actions need to be taken to protect current urban spaces, reduce vegetation depletion, and open space reduction. At the same time, efforts are needed to reduce SUHI by improving building design using green city technologies. At the same time, reduction in runoff and improved rainwater harvesting will be essential at both local and large scales through the participation of individuals, private organizations, and local to national government.