Remote Sensing-Based Quantification of the Relationships between Land Use Land Cover Changes and Surface Temperature over the Lower Himalayan Region

Population growth and population inflow from other regions has caused urbanization which altered land use land cover (LULC) in the lower Himalayan regions of Pakistan. This LULC change increased the land surface temperature (LST) in the region. LULC and LST changes were assessed for the period of 1990–2017 using Landsat data and the support vector machine (SVM) method. A combined cellular automata and artificial neural network (CA-ANN) prediction model was used for simulation of LULC changes for the period of 2032 and 2047 using transition potential matrix obtained from the data years of 2002 and 2017. The accuracy of the CA-ANN model was validated using simulated and classified images of 2017 with correctness value of 70% using validation modules in QGIS. The thermal bands of Landsat images from the years 1990, 2002 and 2017 were used for LST derivation. LST acquired for this period was then modeled for 2032 and 2047 using urban indices (UI) and linear regression analysis. The SVM land cover classification results showed a 5.75% and 4.22% increase in built-up area and bare soil respectively, while vegetation declined by 9.88% during 1990–2017. The results of LST for LULC classes showed that the built-up area had the highest mean LST as compared to other classes. The future projection of LULC and LST showed that the built-up area may increase by 12.48% and 14.65% in 2032 and 2047, respectively, of the total LULC area which was ~11% in 2017. Similarly, the area with temperature above 30 ◦C could be 44.01% and 58.02% in 2032 and 2047, respectively, of the total study area which was 18.64% in 2017. This study identified major challenges for urban planners to mitigate the urban heat island (UHI) phenomenon. In order to address the UHI in the study area, an urban planner might focus on urban plantation and decentralization of urban areas.


Introduction
The land use land cover (LULC) changes have become a source of concern because of its role in reducing biodiversity and altering regional climate and creating an urban heat island (UHI) in cities [1]. Rapid urbanization caused LULC changes which increased land surface temperature (LST) [2]. The land cover indicates the biophysical properties of the land surface and land use is the anthropogenic activity placed on a specific land; due to their interconnectivity they are studied together [3]. Human migration to cities causes urbanization which brings rapid changes in the ecosystem, biodiversity and the environment [4]. Rapid urbanization has both positive and negative impacts. The benefits include additional job opportunities and better life quality. The negative impacts are environmental pollution, burden on city infrastructure and health problems. The increase in built up areas in the urban center could cause an increase in LST which could be a concern for geographers, climatologists and urban planners [5,6].
The land surface consists of vegetation, water surfaces, impervious surface materials, soil and rocks. The higher trend of land surface temperature of urban areas was due to impervious surfaces which contribute towards the formation of an urban heat island (UHI) [7]. Tropical and sub-tropical climate urban areas LST depend on LULC, where built-up areas have the highest LST in comparison with arid environments in which bare soil had the highest LST [8,9]. Urban warming increased air temperature in the Himalayan region by affecting humidity and temperature of the lower boundary layer of air [10,11]. Also, the increase in LST regimes would cause the formation of an UHI [12]. The urban areas were projected to increase as they occupied only 3% of total land cover in 1950 but expected to cover 66% by 2050 [13] globally. Such rapid urban expansions may replace vegetation with built-up areas which increases LST. This results in UHI creation in urban areas of the world [14][15][16].
Many methods and algorithms were developed for LULC modeling. The most common LULC prediction models are the Markov chain [17], cellular automata (CA) [18,19] and artificial aeural network (ANN) [20]. The Markov chain is a stochastic-based model and commonly used for short term projection of a large area [21][22][23]. The CA model is a model with set transition rules and work on the principle of a neighborhood where the current state of each cell depends on the previous one and its functionality pattern is similar to LULC [24][25][26]. The validation of the CA model is done through Kappa coefficient by comparing observed and simulated LULC map [27]. The simulation of LULC changes using a cellular regression model in the lower Himalayan region for the year 2032 and 2047 showed 84% accuracy [28]. The combination of models such as CA-Markov, CA-Stochastic and CA-ANN could predict multi-directional change and give better results as compared to previously used models because the latter only deal with one-way transitional change [29,30].
Many studies indicated the relationship between LULC and LST [25,[31][32][33][34]. Some researchers studied long-term and seasonal changes in temperature in relation to urban growth [35,36]. Other studies showed the quantitative relationship between LULC and temperatures without future temperature prediction [37][38][39]. The relationship between LST and a variety of LULC indices were known to be strong. Therefore, LULC indices such as normalized difference built-up index (NDBI) and fractional vegetation cover can be used for future LST projection. However only a few researchers used LULC indices for future LST projection.
Although strong correlation between LST and LULC indices were reported in several studies, only a few researchers used LULC indices for LST projection. For example, normalized difference vegetation index (NDVI) was used for simulation of urban LULC indices and LST in Dhaka Bangladesh [40]. It was reported that NDVI is not a good predictor of LST in comparison to NDBI and the percentage of impervious surface area (ISA) [32,41]. The combination of indices such as normalized difference bareness index (NDBaI), soil adjusted vegetation index, NDBI, NDVI, urban index (UI), water index, bareness index and enhanced built-up was also used for LST projection [42]. The results showed that a combination of indices in regression analysis affected model predication accuracy due to multicollinearity [40]. Several indices which included both vegetation and nonvegetation were used for LST prediction and it was found that the UI is a better predictor of LST as compared to other indices [43].
There was huge influx of population in the study area from the neighboring districts in the past three decades which caused LULC changes. Therefore, it was important to detect the past patterns of LULC and LST in this region and predict the future changes for a proper urban planning and natural resource management. This study was designed to estimate the current and future LULC and LST using a combination of remote sensing data, urban indices, modeling techniques and linear regression methods. The specific objectives of this study were to: (i) assess LULC and LST changes using Landsat data and support vector machine (SVM) method, and (ii) simulate/predict the future changes in LULC and LST regimes using the CA-ANN model and indices-based regression analysis, respectively, for the years 2032 and 2047.

Study Area
This study considers parts of Mansehra and Battagram districts ( Figure 1) situated in Northern Pakistan, at 34°15′ to 34°45′ N latitude and 72°45′ to 73°20′ E longitude. It is located in the lower Himalayan region on the route of the China-Pakistan Economic Corridor. It has a total area of 1802 km² with a population of 1.8 million [44]. According to Köppen and Geiger, the climate of the study area is classified as Cfa which is humid subtropical or mild temperate climate [45]. The mean annual temperature is 18.5 °C and the total amount of annual rainfall is 1445 mm [46]. The study area is a tourist spot and works as a transit point to Northern Pakistan. The Karakoram Highway (KKH; Figure 1) passes through this region, which was a part of the famous ancient silk road. The people migrated from earthquake-affected districts and Afghanistan due to the Soviet-Afghan War. This migration brought significant LULC changes in the study area.

Data
In this study, three Landsat images available from the United States Geological Survey were used during 1990, 2002 and 2017 (see Table 1 for details) for evaluating changes in LULC and LST regimes. The entire Landsat scene cloud cover for the years of 1990, 2002 and 2017 was about 1%-13% but it was less than 1% in the study area (Table 1). Before LULC classification, a pretreatment process was applied in order to remove the atmospheric effects from Landsat images [47]. For LULC classification accuracy, a ground survey was conducted to collect 40 points for each land cover class.
The classification accuracy was above 90% for all three classified images. The spatial variables data for LULC simulation such as road shape file and DEM was collected from local administration and the United States Geological Survey online data portal.

Classification Schema for LULC and Assessment of Its Accuracy
The Anderson classification scheme (Level 1) was used for LULC classification which is a commonly used classification system for Landsat data [48,49]. The acquired Landsat images were classified for the years of 1990, 2002, and 2017 using the SVM method [50,51] to generate the following land cover target classes: (i) build-up (consisting of industrial, residential, commercial infrastructure, and transportation networks); (ii) bare soil (fallow land, sand, and vacant land); (iii) vegetation (parks, trees, grasslands, and playgrounds); and (iv) water body (coastal water, lakes, canals, and streams). In order to evaluate the classification accuracy, ground truth data was obtained from the 40 land cover reference points for validation and these were overlaid on high resolution satellite images of Google Earth. For better accuracy both spectral and spatial profiles and ancillary information such as Google Earth images were used for training site development. The accuracy of the classification was assessed with kappa coefficient and overall accuracy was obtained using a confusion matrix which is a widely used method for classification accuracy [52].

Estimation of LST
The LST was derived from geometrically corrected Landsat satellite images for the years of 1990, 2002 and 2017. The following set of equations was used for LST calculation [53]. Firstly, the pixelspecific digital number (DN) was transformed into the radiance (Lλ) as follows:

255
(1) where Lmaxλ and Lminλ are the maximum and minimum radiance values, respectively; and their values were available from the metadata of the Landsat images. Secondly, the Lλ values were converted into brightness temperature (TB) as follows: where K1 and K2 are constant and available from the United States Geological Survey. Thirdly, the TB values were translated into LST as follows: where λ is the center of the wavelength, σ is the Boltzmann constant (= 1.380 649 × 10 −23 J/K), h is the Planck's constant (= 6.626 × 10 −34 J/s), c is the velocity of light (= 2.998 × 10 8 m/s), and ε is the emissivity that is being derived as follows [54]: where PV is the proportion of vegetation and calculated as a function of NDVI values as follows: (5) where NDVImax and NDVImin are the image-specific maximum and minimum of the NDVI values, respectively. Also, the calculation of NDVI is found in Table 2.
Note that the Landsat 8 has two thermal bands centered at 10.9 μm and 12 μm. Therefore, LST for Landsat 8 was calculated by averaging LST values from both of the bands.

Normalization of LST
As the LST images were generated for three different years, normalizing them prior to further analysis was required. The LST image of 2017 was used as the reference in this normalization process. The z-score method was employed for this purpose with the following expression as modified after Salama et al. [55]: Several prediction models were employed previously for modeling the future LULC [56]. In the current study, modeling of LULC changes for 2032 and 2047 was completed using a combined CA-ANN model in QGIS. It first modeled the transition potential matrix using ANN, then predicted future LULC changes using the CA model. A combination of models such as CA and ANN (CA-ANN) can predict multi-directional change and give better results than the models which deal only with one-way transitional change [29,30]. The LULC for the years 2032 and 2047 was predicted using the 2002 and 2017 year-specific classified images. A similar method was used by [57] for prediction of LULC for the year 2028 using classified images from 2002 and 2015 in the Coastal Assasuni Upazila area in south-western Bangladesh. Variables of elevation and distance from main roads were used as land use dynamics or drivers. The land use drivers and classified images from 2002 and 2017 were used to generate a transition potential matrix using the QGIS MOLUSCE tool. Once the transition potential matrix was generated then this was modeled using ANN. Finally, LULC projection for 2032 and 2047 was obtained by the CA model. The validation of the CA-ANN model is important; therefore, the CA-ANN model was validated by comparing simulated LULC for 2017 with estimated LULC of 2017 using MOLUSCE QGIS validation.

Modeling of LST for 2032 and 2047
After simulating the LULC for 2032 and 2047, the LULC indices-based methodology proposed by [43] was modified in this study to project LST for the years 2032 and 2047. The different indices which were used for forecasting LST change were NDBaI, NDVI, NDBI and UI. The indices-specific equations are provided in Table 2.
After derivation of LULC indices, correlation (bivariate) analysis between indices and LST was performed in order to select the indices that had strong correlation with LST. The correlation analysis between indices and LST showed that the UI significantly correlated with LST (p < 0.05). The values of correlation coefficient (R) and Root Mean Square Error (RMSE) were 0.410 and 0.094, respectively.
A regression model with a selected index (i.e., UI), was then used to project the LST for 2032 and 2047. The regression analysis was conducted between UI and LST in order to formulate the regression equation. The following regression equation was formulated using UI and LST from the year 2017. LST = 24.25 + (12.41 UI) (6) (6) and predicted/simulated UI indices were employed to calculate LST for the years 2032 and 2047.

Past Pattern of LULC Classes
LULC maps were classified using the SVM method, and its accuracy was evaluated by kappa coefficient using the confusion matrix (see Table 3 for details). The kappa coefficient value was above 0.87 for all three classified images. The overall classification accuracy was 94.96%, 92.26% and 91.35% for the years of 1990, 2002 and 2017, respectively ( Table 3). The results are given in Figure 2 and Table  4. The outcomes indicated that most of the built-up area increase happened in the urban parts of the districts of Mansehra and Battagram ( Figure 2). The built-up area increased from 5.75% to 10.87% of total area during the period of 1990-2017 (Table 4). Table 4 indicated that vegetation declined significantly by 175.77 km 2 (−9.88%) during the study period. The built-up area was expanded by 103.87 km 2 (5.75%). The amount of bare soil increased by 77.8 km 2 (4.22%) from the year of 1990 to 2017.

Past Pattern of LST
The mean LST of individual LULC classes are given in Figure 3. Overall its trend from higher to lower range was for built-up, bare soil, vegetation and water body.  The areal distribution of LST among the different temperature ranges for 1990, 2002 and 2017 is given in Figure 4. The results indicated that in 1990 only 0.8% area had a temperature above 30 °C which increased to 18.64% in 2017. The LST increased from 10.53% to 16.54% for the temperature of 27 °C to 30 °C. The decreasing trend was observed for temperature below 24 °C. For example, the temperature ranged from 21 °C to 24 °C decreased from 31.52% to 20.04%. A similar trend was observed with temperature range of 15 °C to 21 °C. The increase in LST from 2002 to 2017 was greater than from 1990 to 2002. The spatial distribution of the LST in the study area for 1990, 2002 and 2017 is given in Figure 5. It showed that high temperature regimes were mostly found in the urban area (Tehsil Mansehra) as compared to sub-urban area such as Oghi and Battagram (Figure 1).   Figure 6 demonstrates the simulated/predicted LULC patterns for the years 2032 and 2047 and areal statistics are given in Table 5. The accuracy of the CA-ANN model showed %-correctness value of 70% which is an acceptable accuracy level. The results showed that the total built-up area could be ~12.5% and ~15% in 2032 and 2047, respectively, which was ~11% in 2017. The projected results indicated that vegetation could be 45.03% and 42.61% in 2032 and 2047, respectively, which was ~47% in 2017. The bare soil could change to almost 42.35% and 42.63% of the total area in 2032 and 2047, respectively, which was 42.14% in 2017. The bare soil showed an increasing trend for the projected time period. There was no significant change in water body from the base year of 2017. The water content will be 0.12% and 0.08% in 2032 and 2047, respectively, which was 0.23% of the total land cover in 2017.   Figure 7 shows the simulated LST for 2032 and 2047 and areal statistics are provided in Table 6. The LST range above 30 °C was 44.01% and 58.02% in 2032 and 2047, respectively, which was 18.64% in 2017. The mean LST for 2032 and 2047 was above 30 °C for the study area which was 18.64 °C in 2017. The results indicated that the area in the LST range less than 30 °C decreased while above 30 °C showed an increasing trend in the modeling forecast.

Change in LULC
The current research study examined LULC changes during the 1990-2017 period. The results indicated a 2.31% increase in built-up area from 2002 and 2017. The built-up area was increased by 5.75% from 1990 to 2017. The built-up and bare soil increased at the expense of vegetation which declined by −9.98% during the 1990-2002 period, and the built-up area expanded from 93.11 km 2 to 134.85 km 2 which was less as compared to the built-up during the period of 2002-2017. Similar decadal growth was noted in other South Asian cities such as Kathmandu [61] and Dhaka [62]. Several factors contributed to urbanization such as population growth, the political and economic situation, and the migration of people from surrounding areas especially Balakot after the 2005 earthquake [63]. The natural increase in Pakistan's general population could also be a major reason for urbanization. The population of Pakistan has been growing by 2% on a yearly basis (four children per mother); if this remains constant, Pakistan's population could be nearly 300 million by 2030, and 450 million by 2050. The Soviet-Afghan War from 1979-1989 in Afghanistan led to a massive influx of migrants from across the border into urban areas of Pakistan. Nearly four million Afghans migrated into north-western Pakistan during the Soviet-Afghan War. Initially they settled in border refugee camps, but later on they migrated into cities. As suggested by different authors, a combination of political and economic factors [64] and high population growth rate [65] were the major reasons for the urban growth in general, with our study area as not an exception, but rather a hot spot. These factors caused a speedy growth of urban areas and a combination of urban and rural LULC categories became the key features in our area of interest. The newly built-up areas were more frequently witnessed, specifically in the sub-urban zones. Built-up areas and bare soil mostly replaced the areas previously covered by vegetation but located at a lower altitude near the built-up areas.

Change in LST
Our study found that built-up areas had the highest land surface temperature followed by bare soil, vegetation and water. Similar findings were reported by [64] who assessed the LULC impact on the land surface temperature in Zhujiang (China), in order to assess LULC and thermal signature relationship. These outcomes were similar to [66], who found that bare soil and built-up area had a higher LST as compared to water bodies and vegetation categories. The current study showed that the LST regimes expanded due to conversion of vegetative surfaces into non-vegetated areas. The findings of the current study also indicated that LST increased for all LULC classes from 1990 to 2017 (Figures 3 and 5) which was probably due to global warming and surface modification. According to Terando et al. [67], LST could increase when urban expansion cascaded with climate change [68]. Similar findings in the lower Himalayan region were also reported [28]. The increase in LST in the study area was higher due to speedy urbanization from 1990 to 2017. Another fact, such as different temporal resolution of Landsat images could contribute for the LST increase beyond the built-up area; however, the outcomes showed that the area was affected by both local warming and/or climate change.

Modeling of LULC Changes for the Year 2032 and 2047
The LULC of the study area is expected to change over the next three decades, as 12.48% and 14.65% of the current land use land cover could be transformed into a built-up class in 2032 and 2047, respectively, which was 10.87% in 2017. A similar trend of built-up areas was also noticed in future forecasting for other growing cities, like Beijing in China [69]. The increasing built-up areas might be associated with the population growth during different decades since 1990, as found in the case of Bangladesh [70]. Rapid urbanization creates an urban sprawl which has both positive and negative impacts. The positive impacts include industrialization which increases employment opportunities to local residents and others living in the surrounding areas. The urban expansion could increase business opportunities and improve health facilities. The construction and establishment of new housing societies could increase in the study area as well. The construction of new luxurious residential societies could provide comfortable and improved living standards to people.
Some of the negative impacts of urbanization include unemployment, extra fuel consumption and several health-and environment-related issues. With the growing population, there is an increasing need of employment opportunities as Pakistan's population is projected to be 306 million in 2050, which could be 3.6% of the total world population [71]. Similarly, the population increase and city expansion increases traveling distance for residences which may cause additional fuel usage and traffic jams. It could increase air pollution and impact numerous health-related issues, especially for the children and the senior citizens in our area of interest. This could also increase pressure on public services. Urbanization also causes social disparities among residents [72].

Modeling of LST Changes for the Year 2032 and 2047
The LST modeled in this study for 2032 and 2047 demonstrated that the LST regimes could be over 27 °C for most of the study area. The predicted temperature increase due to rapid urban growth concurs with the variety of regional and global models which predicted similar scenarios for urban growth [73][74][75]. This is due to high energy demands caused by population increase. The urban expansion effects greenhouse gases which are mostly released from different human activities. In fact, enormous amounts of anthropogenic carbon dioxide emissions are usually the case in urban areas as a result of burning of fossil fuels for: (i) both cooling and heating; (ii) industrial production; (iii) vehicles. For example, Kahn [76] noted that the greenhouse gas concentration and anthropogenic heat emissions could intensify due to urbanization as the population may grow in the middle class. The income increment in the middle class could raise the consumption of gadgets such as air conditioners, fridges, stoves and heaters. Forest work as a major carbon sink which was affected by LULC changes primarily by deforestation at the regional scale, caused changes in surface characteristics which could be the reason for an increase in LST along with other factors such as land utilization for roads and cities and goods and services demands by urban residents [77]. The mean annual temperature in Pakistan has increased by 0.57 °C during the past 100 years and it will further increase by 3 °C to 5 °C as per central global emissions scenario, while as per higher global emissions scenario, it will rise by 4 °C to 6 °C [78]. Since air temperature and LST are correlated, therefore an increase in air temperature would lead to an increase in LST and land surface temperature would be greater than air temperature [79]. The increase in temperature could have negative impacts on plants and animal species as well as human health.

Conclusions
The main purpose of this study was to assess the effect of LULC on LST during the period 1990-2017 in the study area. Landsat data and support vector machine method were used to generate LULC classification maps. The past trends of LULC and LST were employed in CA-ANN and regression models to predict the future trends of these two variables. The results indicated significant changes in the LULC, in particular the built-up and bare soil area were amplified by 5.75% and 4.22% respectively, while vegetation declined by 9.88% during the past data period of 1990-2017. The zonal analysis between LST and LULC indicated that the mean LST for a built-up area was high as compared to other classes. The land surface temperature increased for all land cover classes which was due to both urban warming and climate change. The future simulation results showed that the built-up area is expected to encompass 12.48% and 14.65% of the total study area in 2032 and 2047, respectively, as compared to ~11% in 2017. Similarly, the area with an LST above 30 °C could increase to ~44% and ~58% of the total area in 2032 and 2047, respectively, as compared to ~18% in 2017. It was noted that lower temperature zones are shifting towards higher temperatures which could lead to large-scale UHI formation. This study identified major challenges for urban planners to mitigate UHIs in the study area. Urban plantation and decentralization of urban areas are two possible solutions in order to mitigate the UHI and its consequences. Future urban research could focus on the issue of public health and infrastructure burden associated with rapid urbanization.