Modeling Spatio-temporal Land Transformation and Its Associated Impacts on land Surface Temperature (LST)

: Land use land cover (LULC) of city regions is strongly affected by urbanization and affects the thermal environment of urban centers by influencing the surface temperature of core city areas and their surroundings. These issues are addressed in the current study, which focuses on two provincial capitals in Pakistan, i.e., Lahore and Peshawar. Using Landsat data, LULC is determined with the aim to (a) examine the spatio-temporal changes in LULC over a period of 20 years from 1998 to 2018 using a CA-Markov model, (b) predict the future scenarios of LULC changes for the years 2023 and 2028, and (c) study the evolution of different LULC categories and investigate its impacts on land surface temperature (LST). The results for Peshawar city indicate the significant expansion in vegetation and built-up area replacing barren land. The vegetation cover and urban area of Peshawar have increased by 25.6%, and 16.3% respectively. In contrast, Lahore city urban land has expanded by 11.2% while vegetation cover decreased by (22.6%). These transitions between LULC classes also affect the LST in the study areas. Transformation of vegetation cover and water surface into built-up areas or barren land results in the increase in the LST. In contrast, the transformation of urban areas and barren land into vegetation cover or water results in the decrease in LST. The different LULC evolutions in Lahore and Peshawar clearly indicate their effects on the thermal environment, with an increasing LST trend in Lahore and a decrease in Peshawar. This study provides a baseline reference to urban planners and policymakers for informed decisions.

have dramatically influenced the environment [65] and often threatens the sustainable urban development [66,67].
In recent years, Lahore, in the province of Punjab (see Figure 1), has experienced massive growth because of the high demand for the urban area. Most of the agriculture and barren land was bought by local developers and the Lahore Development Authority (LDA) and turned into different housing societies and other projects [68,69]. As a result of this urban expansion and over-cutting, the vegetation cover of Lahore has decreased substantially. Moreover, natural land has been replaced mainly by built-up areas, i.e., urbanized areas including buildings, roads, and different other infrastructures [68,69].

Data Collection and Preprocessing
For the analysis of the LULC dynamics and the effect of the LULC transition on the thermal environment in Lahore and Peshawar, satellite imagery including Landsat 5 Thematic Mapper (TM) and Landsat 8 (Operational Land Imager (OLI) was collected ( Table 1). The temporal resolution of the images is 16 days and the spatial resolution is 30 m. All images were freely downloaded from the USGS earth explorer website (http://earthexplorer.usgs.gov/, last access: 4 January 2020). After the preparation of the satellite images vector layer of the administrative boundary of the study area, Lahore and Peshawar were utilized as masks to subset the images for clipping the area of interest (AOI) from the Tagged Image File Format (TIFF). No atmospheric corrections were executed since the Landsat images were cloud-free [75,76]. After this, the thermal bands were used to derive the LST [77,78]. First, the digital number (raw data) of Landsat images thermal bands was converted into radiance and then surface temperature. For LULC maps, visible, near-infrared and shortwave infrared bands of Landsat products were stacked and mosaicked using ARDAS Imagine [79,80]. Supervised classification was used to classify the pixels into four different LULC classes [81]: vegetation, water bodies, built-up, and barren land. After the preparation of LULC maps for each year, the error matrix was constructed to evaluate the accuracy of LULC classes. To this end, ground truth samples for each land-use class were taken from random locations using Google Earth [82] and compared with the collocated LULC pixel. The percentage accuracy of each class was evaluated in ArcGIS using frequency analysis and error matrix [83]. The results are presented in Section 3. The LST retrieval is visualized in the flowchart of Figure 2, with details described below. Peshawar, in the province of Khyber Pakhtunkhwa (KPK, see Figure 1), has also experienced sudden changes in population growth due to the migration of people from Afghanistan [70]. Following the 1961 census, Peshawar represents 29% of the total population of the province of Khyber Pakhtunkhwa, whereas in 1998 the population in Peshawar had increased to 33% of the entire population of this province [71]. Additionally, in 2013, the provisional government of KPK started the "billion trees project" in the whole region [72], which aims to restore 150 million hectares of the world's degraded and deforested lands by 2020 (https://en.unesco.org/courier/2019-3/pakistangreen-again/, last access: 20 august 2020), and Pakistan hit its "billion trees" goal in August 2017 (https://www.weforum.org/agenda/2018/07/pakistan-s-billion-tree-tsunami-is-astonishing/, last access: 20 August 2020). Such changes in LULC influence the LST, air temperature, and topography of the neighborhood [73,74]. Therefore, there is a strong need to quantify the individual contribution of different LULC transformations to LST in and around urbanized regions of Pakistan.
In this study, the CA-Markov chain model is applied to simulate the current and future land use dynamics over two metropolitan cities in Pakistan, Lahore and Peshawar. As described above, land use management in these cities was very different, which resulted in contrasting effects on the thermal environment. Overall, this study was conducted to (a) examine the spatio-temporal LULC dynamics for the years 1998, 2003, 2008, 2013, and 2018, (b) predict the LULC in 2023 and 2028 using the Markov and CA-Markov models, and (c) examine the LULC transition into different LULC categories and to quantify its impacts on LST. The results clearly show the effect of different land-use policies on the urban environment and provide direction for reducing environmental problems in future climate scenarios using different choices of LULC. The study area, data collection and models used are described in Section 2, including the motivation for using the CA-Markov model to simulate LULC changes (Section 2.4). The results are presented and discussed in Section 3 and the conclusions from the study are provided in Section 4.

Study Area
The study area comprises two large cities (Lahore and Peshawar) in Pakistan ( Figure 1). Peshawar is the capital of KPK (Khyber Pakhtunkhwa) Province, along the Khyber Pass, near the border of Afghanistan. Peshawar is a magnificent and essential economic, political and military center of the province and Pakistan. Geographically Peshawar is located at 44 • 15' N, 71 • 42' E and covers an area of 1264 km 2 . According to the census of 2017, a total of 1.97 million people live in Peshawar (http://www.pbs.gov.pk/ last access: 13 July 2020). Lahore is the 2nd largest city in Pakistan and is the capital city of Punjab province. Lahore covers a total area of 1842 km 2 located between the latitudes of 31 • 15 -31 • 43 N and longitudes 74 • 10 -74 • 39 E. A total of 11.13 million people live in Lahore according to the statistics of the 2017 census (http://www.pbs.gov.pk/ last access: 13 July 2020).

Data Collection and Preprocessing
For the analysis of the LULC dynamics and the effect of the LULC transition on the thermal environment in Lahore and Peshawar, satellite imagery including Landsat 5 Thematic Mapper (TM) and Landsat 8 (Operational Land Imager (OLI) was collected ( Table 1). The temporal resolution of the images is 16 days and the spatial resolution is 30 m. All images were freely downloaded from the USGS earth explorer website (http://earthexplorer.usgs.gov/, last access: 4 January 2020). After the preparation of the satellite images vector layer of the administrative boundary of the study area, Lahore and Peshawar were utilized as masks to subset the images for clipping the area of interest (AOI) from the Tagged Image File Format (TIFF). No atmospheric corrections were executed since the Landsat images were cloud-free [75,76]. After this, the thermal bands were used to derive the LST [77,78]. First, the digital number (raw data) of Landsat images thermal bands was converted into radiance and then surface temperature. For LULC maps, visible, near-infrared and shortwave infrared bands of Landsat products were stacked and mosaicked using ARDAS Imagine [79,80]. Supervised classification was used to classify the pixels into four different LULC classes [81]: vegetation, water bodies, built-up, and barren land. After the preparation of LULC maps for each year, the error matrix was constructed to evaluate the accuracy of LULC classes. To this end, ground truth samples for each land-use class were taken from random locations using Google Earth [82] and compared with the collocated LULC pixel. The percentage accuracy of each class was evaluated in ArcGIS using frequency analysis and error matrix [83]. The results are presented in Section 3. The LST retrieval is visualized in the flowchart of Figure 2, with details described below.   The Landsat data are downloaded as raw data (digital numbers) which needs to be converted to spectral radiances (Lλ) by using the information in the LANDSAT metadata header file [16]: where L is the spectral radiance (Wm −2 sr −1 µm −1) , LMAXλ is the spectral radiance scaled to QCalmax (W m −2 sr −1 µm −1 ). LMINλ is the spectral radiance scaled to QCalmin (W m −2 sr −1 µm −1 ). QCalmax is the maximum quantized calibrated pixel value (DN = 255) that corresponds to LMAXλ. The minimum

Conversion of Raw Landsat Data into Radiance
The Landsat data are downloaded as raw data (digital numbers) which needs to be converted to spectral radiances (Lλ) by using the information in the LANDSAT metadata header file [16]: where L λ is the spectral radiance (Wm −2 sr −1 µm −1) , LMAX λ is the spectral radiance scaled to Q Calmax (W m −2 sr −1 µm −1 ). LMIN λ is the spectral radiance scaled to Q Calmin (W m −2 sr −1 µm −1 ). Q Calmax is the maximum quantized calibrated pixel value (DN = 255) that corresponds to LMAX λ . The minimum quantized calibrated pixel value (DN = 0) corresponds to LMIN λ . Q Cal is the quantized calibrated pixel value (DN) [84,85].
where r is the planetary reflectance (dimensionless), Lλ is the spectral radiance at the sensor aperture (W m −2 sr −1 µm −1 ), dr = 1 + 0.033cos (D × 2 × 3.14)/365), where D is the day of the year, Esun is the mean solar atmospheric irradiance (W m −2 µm −1 ), θ is the solar zenith angle (degree), θ = (90 -B), and B is the Sun elevation angle. dr is the inverse square of the earth-sun distance (astronomical unit).

Land Surface Temperature (LST) Retrieval
For the retrieval of LST, Landsat images were used with a temporal resolution of 16 days and a spatial resolution of 30m. The LST was derived from the Landsat thermal bands using the methodology recommended by [86]. In the first step, the radiance of the thermal band calculated using Equation (1) was converted into brightness temperature (TB) using Equation (3) [88]. Land surface emissivity ( Ɛ ) was estimated using the NDVI thresholds method as proposed by [77], using the following equations: where Ɛ is the vegetation emissivity, Ɛs is the soil emissivity, F = 0.55 is the shape factor (Lim et al., 2012), and Pv is the vegetation proportion, which was obtained using Equation (6) [87]: where NDVI is derived by using Equation (7) [89]: In Equation (7), NIR and RED are reflectances. For Landsat5 (TM), NIR is the reflectance measured in band 4 at wavelength λ = 0.76-0.90 µm and RED is the reflectance measured in band 3 with λ = 0.63-0.69 µm. In Landsat8 (OLI), NIR is the reflectance measured in band 5 with λ = 0.85-0.87 µm and RED refers to the reflectance measured in band 4 with λ = 0.64-0.67 µm. NIR and RED values were retrieved using Equation (2).
In this study, the focus is the effect of LULC transition on LST, which is stronger in the hottest months in Pakistan, i.e., May, June, and July. Therefore, LST trends over the study period were evaluated using averages over these 3 months.

Land Surface Temperature (LST) Retrieval
For the retrieval of LST, Landsat images were used with a temporal resolution of 16 days and a spatial resolution of 30m. The LST was derived from the Landsat thermal bands using the methodology recommended by [86]. In the first step, the radiance of the thermal band calculated using Equation (1) was converted into brightness temperature (TB) using Equation (3): where; K1 and K2 are conversion constants. For Landsat8 OLI, K1 = 774.89 mW cm −2 sr −1 µm −1 and K2 = 1321.08 K; Landsat5 TM, K1 = 607.76 mW cm −2 sr −1 µm −1 and K2 = 1260.56 K. In the second step, the land surface temperature was derived using the following relation [87]: where λ (≈ 11.5 µm) is the effective wavelength of the thermal bands.  In Equation (7), NIR and measured in band 4 at wavelen with λ = 0.63-0.69 µm. In Land 0.87 µm and RED refers to the r values were retrieved using Equ In this study, the focus is t months in Pakistan, i.e., May, evaluated using averages over t ) was estimated using the NDVI thresholds method as proposed by [77], using the following equations: d Remote Sens. 2020, 12, x FOR PEER REVIEW quantized calibrated pixel value (DN = 0) corresponds to LMINλ. QCal is the quantized calib value (DN) [84,85]

Conversion of Radiance to Reflectance
The radiances calculated using Equation (1) were converted into reflectance using E [85]: where r is the planetary reflectance (dimensionless), Lλ is the spectral radiance at the sens (W m −2 sr −1 µm −1 ), dr = 1 + 0.033cos (D × 2 × 3.14)/365), where D is the day of the year, Esun i solar atmospheric irradiance (W m −2 µm −1 ), θ is the solar zenith angle (degree), θ = (90 -B the Sun elevation angle. dr is the inverse square of the earth-sun distance (astronomical u

Land Surface Temperature (LST) Retrieval
For the retrieval of LST, Landsat images were used with a temporal resolution of 16 spatial resolution of 30m. The LST was derived from the Landsat thermal bands methodology recommended by [86]. In the first step, the radiance of the thermal band using Equation (1) [88]. Land surface emissivit estimated using the NDVI thresholds method as proposed by [77], using the following eq where Ɛ is the vegetation emissivity, Ɛs is the soil emissivity, F = 0.55 is the shape factor 2012), and Pv is the vegetation proportion, which was obtained using Equation (6) [87]: where NDVI is derived by using Equation (7) [89]: In Equation (7), NIR and RED are reflectances. For Landsat5 (TM), NIR is the measured in band 4 at wavelength λ = 0.76-0.90 µm and RED is the reflectance measured with λ = 0.63-0.69 µm. In Landsat8 (OLI), NIR is the reflectance measured in band 5 wit 0.87 µm and RED refers to the reflectance measured in band 4 with λ = 0.64-0.67 µm. NI values were retrieved using Equation (2).
In this study, the focus is the effect of LULC transition on LST, which is stronger in months in Pakistan, i.e., May, June, and July. Therefore, LST trends over the study p evaluated using averages over these 3 months.

Conversion of Radiance to Reflectance
The radiances calculated using Equation (1) were converted into reflectance [85]: where r is the planetary reflectance (dimensionless), Lλ is the spectral radiance at t (W m −2 sr −1 µm −1 ), dr = 1 + 0.033cos (D × 2 × 3.14)/365), where D is the day of the yea solar atmospheric irradiance (W m −2 µm −1 ), θ is the solar zenith angle (degree), θ the Sun elevation angle. dr is the inverse square of the earth-sun distance (astrono

Land Surface Temperature (LST) Retrieval
For the retrieval of LST, Landsat images were used with a temporal resolutio spatial resolution of 30m. The LST was derived from the Landsat thermal methodology recommended by [86]. In the first step, the radiance of the therm using Equation (1)  , ℎ is the Planck constant (6.6 is the speed of light (3.0 × 10 8 ms −1 ), ϵ is the land surface emissivity with values of 0 for vegetation, build-up and water surfaces, respectively [88]. Land surface e estimated using the NDVI thresholds method as proposed by [77], using the follo where Ɛ is the vegetation emissivity, Ɛs is the soil emissivity, F = 0.55 is the shap 2012), and Pv is the vegetation proportion, which was obtained using Equation (6 where NDVI is derived by using Equation (7) [89]: In Equation (7), NIR and RED are reflectances. For Landsat5 (TM), NIR measured in band 4 at wavelength λ = 0.76-0.90 µm and RED is the reflectance m with λ = 0.63-0.69 µm. In Landsat8 (OLI), NIR is the reflectance measured in ban 0.87 µm and RED refers to the reflectance measured in band 4 with λ = 0.64-0.67 values were retrieved using Equation (2).
In this study, the focus is the effect of LULC transition on LST, which is stro months in Pakistan, i.e., May, June, and July. Therefore, LST trends over the s evaluated using averages over these 3 months. The radiances calculated using Equation (1) were conver [85]: , θ is the solar zenit the Sun elevation angle. dr is the inverse square of the earth-su

Land Surface Temperature (LST) Retrieval
For the retrieval of LST, Landsat images were used with spatial resolution of 30m. The LST was derived from th methodology recommended by [86]. In the first step, the rad using Equation (1)  , ϵ is the land surface emissiv for vegetation, build-up and water surfaces, respectively [8 estimated using the NDVI thresholds method as proposed by where Ɛ is the vegetation emissivity, Ɛs is the soil emissivity 2012), and Pv is the vegetation proportion, which was obtaine where NDVI is derived by using Equation (7) [89]: In Equation (7), NIR and RED are reflectances. For L measured in band 4 at wavelength λ = 0.76-0.90 µm and RED with λ = 0.63-0.69 µm. In Landsat8 (OLI), NIR is the reflectan 0.87 µm and RED refers to the reflectance measured in band 4 values were retrieved using Equation (2).
In this study, the focus is the effect of LULC transition on months in Pakistan, i.e., May, June, and July. Therefore, LS evaluated using averages over these 3 months. The radiances calculated using Equation (1) were converted into reflectance using Equation (2)  [85]: where r is the planetary reflectance (dimensionless), Lλ is the spectral radiance at the sensor aperture (W m −2 sr −1 µm −1 ), dr = 1 + 0.033cos (D × 2 × 3.14)/365), where D is the day of the year, Esun is the mean solar atmospheric irradiance (W m −2 µm −1 ), θ is the solar zenith angle (degree), θ = (90 -B), and B is the Sun elevation angle. dr is the inverse square of the earth-sun distance (astronomical unit).

Land Surface Temperature (LST) Retrieval
For the retrieval of LST, Landsat images were used with a temporal resolution of 16 days and a spatial resolution of 30m. The LST was derived from the Landsat thermal bands using the methodology recommended by [86]. In the first step, the radiance of the thermal band calculated using Equation (1) was converted into brightness temperature (TB) using Equation (3): where; K1 and K2 are conversion constants. For Landsat8 OLI, K1 = 774.89 mW cm −2 sr −1 µm −1 and K2 = 1321.08 K; Landsat5 TM, K1 = 607.76 mW cm −2 sr −1 µm −1 and K2 = 1260.56 K. In the second step, the land surface temperature was derived using the following relation [87]: where λ (≈ 11.5 µm) is the effective wavelength of the thermal bands. ρ= (hc)/σ = 1.438 × 10 −2 mK, where σ is the Boltzmann constant (1.38 × 10 −23 JK −1 ), ℎ is the Planck constant (6.626 × 10 −34 Js) and c is the speed of light (3.0 × 10 8 ms −1 ), ϵ is the land surface emissivity with values of 0.95, 0.92 and 0.9925 for vegetation, build-up and water surfaces, respectively [88]. Land surface emissivity ( Ɛ ) was estimated using the NDVI thresholds method as proposed by [77], using the following equations: where Ɛ is the vegetation emissivity, Ɛs is the soil emissivity, F = 0.55 is the shape factor (Lim et al., 2012), and Pv is the vegetation proportion, which was obtained using Equation (6) [87]: where NDVI is derived by using Equation (7) [89]: In Equation (7), NIR and RED are reflectances. For Landsat5 (TM), NIR is the reflectance measured in band 4 at wavelength λ = 0.76-0.90 µm and RED is the reflectance measured in band 3 with λ = 0.63-0.69 µm. In Landsat8 (OLI), NIR is the reflectance measured in band 5 with λ = 0.85-0.87 µm and RED refers to the reflectance measured in band 4 with λ = 0.64-0.67 µm. NIR and RED values were retrieved using Equation (2).
In this study, the focus is the effect of LULC transition on LST, which is stronger in the hottest months in Pakistan, i.e., May, June, and July. Therefore, LST trends over the study period were evaluated using averages over these 3 months.

Conversion of Radiance to Reflectance
The radiances calculated using Equation (1) were converted into reflectance using E [85]: where r is the planetary reflectance (dimensionless), Lλ is the spectral radiance at the sens (W m −2 sr −1 µm −1 ), dr = 1 + 0.033cos (D × 2 × 3.14)/365), where D is the day of the year, Esun solar atmospheric irradiance (W m −2 µm −1 ), θ is the solar zenith angle (degree), θ = (90 -B the Sun elevation angle. dr is the inverse square of the earth-sun distance (astronomical u

Land Surface Temperature (LST) Retrieval
For the retrieval of LST, Landsat images were used with a temporal resolution of 16 spatial resolution of 30m. The LST was derived from the Landsat thermal bands methodology recommended by [86]. In the first step, the radiance of the thermal band using Equation (1) [88]. Land surface emissivit estimated using the NDVI thresholds method as proposed by [77], using the following eq is the vegetation emissivity, Ɛs is the soil emissivity, F = 0.55 is the shape factor 2012), and Pv is the vegetation proportion, which was obtained using Equation (6)  In this study, the focus is the effect of LULC transition on LST, which is stronger in months in Pakistan, i.e., May, June, and July. Therefore, LST trends over the study p evaluated using averages over these 3 months. s is the soil emissivity, F = 0.55 is the shape factor (Lim et al., 2012), and Pv is the vegetation proportion, which was obtained using Equation (6) [87]: where NDVI is derived by using Equation (7) [89]: In Equation (7 In this study, the focus is the effect of LULC transition on LST, which is stronger in the hottest months in Pakistan, i.e., May, June, and July. Therefore, LST trends over the study period were evaluated using averages over these 3 months.

Simulation of LULC Changes Using the CA-Markov Model
The CA-Markov model was used for LULC simulations [90] because it provides a powerful technique to study urban growth trends and future predictions [54,[91][92][93] and detailed information on a synoptic scale [50][51][52][53][54]. The CA-Markov model can be applied to simulate future LULC distributions. It consists of two parts, the CA model and a Markov chain model. The Markov chain model estimates the state of the LULC change between two time intervals to predict future changes. This model not only describes the conversion states between the land use types but can also account for the transfer rate among different land use categories. It can be used in spatial modeling for forecasting future LULC [94]. The CA model consists of a regular grid of cells, with properties that can change in time according to fixed rules and depend on the current state of a cell and its neighboring cells. Because of these characteristics, it has also been applied for simulating the LULC process [95,96]. In the CA model used here, the specific transition rules apply to 3 × 3 neighboring cells, i.e., a central cell with 8 neighbors which affect the properties of the central cell. For the application of the CA model to simulate the future LULC, the properties of LULC types should be considered, for instance, the built-up area cannot transfer into water in the near future [50]. The Markov model can be mathematically described as Equation (8): where S represents the land use status at time t, and S (t + 1) is the land-use status at time t + 1, while Pij is the transition probability matrix in a state which is calculated as Equations (9) and (10) [97] ||Pij|| = P 1,1 P 1,2 P 1,N P 2,1 P 2,2 P 2,N P N,1 P N,2 P N,N (0 ≤ Pij ≤ 1) where P is the transition probability; Pij stands for the probability of converting from current state i to another state j in next time; and PN is the state probability of any time. The Low transition will have a probability near (0) and high transition have probabilities near (1) [97].

Validation of the LULC Prediction Model
The Kappa Index of Agreement (KIA) approach was used to check the validity of the CA-Markov model. The Kappa statistic measures the accuracy of a classification relative to a completely random classification on a scale from 0 (utterly random assignment of class labels) to 1 (100% accuracy of class label assignment) [98,99]. KIA can be calculated using Equation (11): where Pr (a) represents the observed agreement, and Pr(e) represents chance agreement. The process is described in more detail in Section 3.2.

Relationship between LST and LULC:
The relation between LST and LULC was calculated by linear regression, using the linear regression expression in the Curve Fit tool by using ARC GIS; Curve Fit is an extension of the GIS application ArcMap that allows the user to run regression analysis on a series of raster datasets (geo-referenced images) (https://www.umesc.usgs.gov/management/dss/curve_fit.html/, last access: 19 January 2020). Later on, the zonal statistic was applied to get the values. The equation of linear regression is described below: where Y is the dependent variable, X is the independent variable, b is the slope of the line and a is the y-intercept.

Land Use Land Cover (LULC) Dynamics
Maps of the LULC in Lahore and Peshawar, for one selected day in each of the years 1998, 2003, 2008, 2013 and 2018 (Table 1), were obtained after preprocessing and supervised classification. The results are presented in Figure 3. The accuracy of the LULC classification was assessed using ground truth data obtained from Google Earth and the results are presented in Table 2. The data in Table 2 show that all LULC classes were classified with an accuracy of better than 70%. The classified images in Figure 3 for the years 1998, 2003, 2008, 2013, and 2018 capture the spatial and temporal characteristics of the LULC changes. Table 3 summarizes the LULC for each type for each year. These data show that in Lahore city, the area covered by water bodies, vegetation, and barren area decreased gradually, while the built-up area increased. However, the LULC in Peshawar city experienced an increase in vegetation and a decrease in the barren land from 1998 to 2018.

Land Use Land Cover (LULC) Dynamics
Maps of the LULC in Lahore and Peshawar, for one selected day in each of the years 1998, 2003, 2008, 2013 and 2018 (Table 1), were obtained after preprocessing and supervised classification. The results are presented in Figure 3. The accuracy of the LULC classification was assessed using ground truth data obtained from Google Earth and the results are presented in Table 2. The data in Table 2 show that all LULC classes were classified with an accuracy of better than 70%. The classified images in Figure 3 for the years 1998, 2003, 2008, 2013, and 2018 capture the spatial and temporal characteristics of the LULC changes. Table 3 summarizes the LULC for each type for each year. These data show that in Lahore city, the area covered by water bodies, vegetation, and barren area decreased gradually, while the built-up area increased. However, the LULC in Peshawar city experienced an increase in vegetation and a decrease in the barren land from 1998 to 2018.  Table 2. Accuracy assessment of land use classification in Lahore and Peshawar (%).    Spatial land transition analysis and change detection matrices were prepared by using the land change models for both cities to understand the land encroachment during the last decades. Figure 4 and Table 4 show the LULC transition of each LULC class from 1998 to 2018 in the cities of Lahore and Peshawar. The results show that in Lahore city, a water area of 23.33 km 2 (1.27%) has been turned into built-up land, and a water area of 11.84 km 2 (0.64%) area has been turned into barren land. Likewise, 150.13 km 2 (8.15%) of the vegetation area has been converted into built-up land. From 1998 to 2018, the built-up area in Lahore city has substantially increased, mostly from vegetation land (8.15%) and barren land (16.26 %). The transition between LULC classes in Peshawar city from 1998 to 2018 shows that most of the barren land has been turned into vegetation, while at the same time the city experienced an increase in the built-up area. In addition, barren land has been turned into built-up area. From 1998 to 2018, 381.80 km 2 (30.21 %) of barren land has been converted into vegetation, and 208.57 km 2 (16.50 %) area of barren land has been turned into built-up area.

Future Land Use Dynamics
LULC maps for 2003 and 2008 were used as a starting point for future predictions using the CA-Markov model (Section 2.4) to produce the transition probability matrix. The transition probability matrix for the period 2003-2008 is presented in Table 5. The numbers in Table 5 indicate the probability that one LULC class was converted into another one between 2003 and 2008. Table 5 was used to predict the LULC maps for Lahore and Peshawar city in 2013. The results are presented in Figure 5. The comparison with satellite-derived LULC maps shows that the spatial patterns of all land use categories in the predicted LULC maps are similar to those in the satellite-derived LULC maps. Kappa statistics were used to assess the accuracy of the predictions for each LULC category. The results in Tables 6 and 7 show that the kappa values are larger than 0.78 for each land use category and the average kappa value for Lahore is 0.95 and for Peshawar 0.89. These high kappa values indicate that the CA-Markov model is suitable for future predictions of LULC maps for the study region.

Future Land Use Dynamics
LULC maps for 2003 and 2008 were used as a starting point for future predictions using the CA-Markov model (Section 2.4) to produce the transition probability matrix. The transition probability matrix for the period 2003-2008 is presented in Table 5. The numbers in Table 5 indicate the probability that one LULC class was converted into another one between 2003 and 2008. Table 5 was used to predict the LULC maps for Lahore and Peshawar city in 2013. The results are presented in Figure 5. The comparison with satellite-derived LULC maps shows that the spatial patterns of all land use categories in the predicted LULC maps are similar to those in the satellite-derived LULC maps. Kappa statistics were used to assess the accuracy of the predictions for each LULC category. The results in Tables 6  and 7 show that the kappa values are larger than 0.78 for each land use category and the average kappa value for Lahore is 0.95 and for Peshawar 0.89. These high kappa values indicate that the CA-Markov model is suitable for future predictions of LULC maps for the study region.    For the prediction of the LULC for the year 2023, the transition probability matrix was calculated from the satellite-derived maps for 2013 and 2018; for 2028, the probability matrix calculated from the satellite-derived map for 2018 and the predicted map for 2023 were used. Figure 6 illustrates the LULC maps for Lahore and Peshawar city simulated for the years 2023 and 2028 and the LULC for each class are summarized in Table 8. The results for Lahore city show that the barren land area will gradually decrease in 2023 and 2028 due to its conversion to built-up area, while the area of vegetation land will remain almost the same. The area covered by barren land was 658.59 km 2 (35.7%) in 2018 and will decrease to 580.10 km 2 (31.49%) in 2028, while the area of built-up land was 755.91 km 2 (41.0%) in 2018 and will increase to 840.67 km 2 (45.63%) in 2028. For Peshawar city, the area covered by water bodies will decrease between 2018 and 2028, i.e., in 2018, the water bodies covered an area of 29.90 km 2 (2.4%), which decreases to 9.

Land Surface Temperature (LST) Variations from 1998-2018
The LST patterns for Lahore and Peshawar for three months (May-July) in the years 1998, 2003, 2008, 2013 and 2018 were calculated as described in Section 2.3. Figure 7a shows the increase in the three-months averaged (May-July) LST in Lahore city during these 20 years. Both the minimum and    In contrast, Figure 7b indicates that in Peshawar city the LST has decreased during the study period, for both the minimum and maximum values. For the years 1998, 2003, 2008, 2013

Correlation between LST and T (Air)
For evaluation of the Landsat-derived LST trend, it is compared to the corresponding 3 months (May-July) air temperatures (Ta) measured at meteorological stations in Lahore and Peshawar ( Table 9). The data in Table 9 show the increase in both LST and air temperature in Lahore, whereas in Peshawar, both the air temperature and the LST decrease during the study period. The correlation of 3 months (May-July) of mean LST and air temperature T (a). The results in Figures 8 and 9 show a good comparison with high correlation coefficients between LST and Ta.
Remote Sens. 2020, 12, x FOR PEER REVIEW 15 of 23 The correlation of 3 months (May-July) of mean LST and air temperature T (a). The results in Figures 8 and 9 show a good comparison with high correlation coefficients between LST and Ta.

Changes in LST in Response to Different Land-use Classes
The change in LST in response to the transformation between different land-use classes was studied by determining relationships between LULC transition from 1998 to 2018 and LST using the curve fit tool and zonal statistic in ArcGIS (Section 2.6).
The results presented in Table 10 and Figure 10 show the correlation between LST and LULC, i.e., the regression coefficient (r). The average change in LST (°C) from 1998 to 2018 due to the transitions between the four different LULC classes is presented in Table 10 for each type of change: the transition from built-up and barren land to vegetation and water cover results in a decrease in LST, whereas a change from water or vegetation cover to built-up or barren land results in an increase in LST. The results presented in Table 10 show that in Peshawar city, the LST was reduced by 0.25 and 0.27 °C in response to the transition from built-up and barren land to vegetation, respectively. In

Changes in LST in Response to Different Land-Use Classes
The change in LST in response to the transformation between different land-use classes was studied by determining relationships between LULC transition from 1998 to 2018 and LST using the curve fit tool and zonal statistic in ArcGIS (Section 2.6).
The results presented in Table 10 and Figure 10 show the correlation between LST and LULC, i.e., the regression coefficient (r). The average change in LST ( • C) from 1998 to 2018 due to the transitions between the four different LULC classes is presented in Table 10 for each type of change: the transition from built-up and barren land to vegetation and water cover results in a decrease in LST, whereas a change from water or vegetation cover to built-up or barren land results in an increase in LST. The results presented in Table 10 show that in Peshawar city, the LST was reduced by 0.25 and 0.27 • C in response to the transition from built-up and barren land to vegetation, respectively. In contrast, in Lahore city, the LST increased by 0.17 • C in response to the transition from vegetation to buildup land and 0.13 • C in response to the transition from water cover to built-up area.

Discussion
This research was carried out to quantify the individual contribution of different LULC transformations to LST in Lahore and Peshawar city (metropolitan areas of Pakistan) from 1998-2018, and future predictions of LULC for 2023 and 2028 by using CA-Markov model. Figure 3a highlights the land cover pattern for Lahore city and shows that the area covered by water bodies, vegetation, and barren land decreased gradually from 1998 to 2018, while built-up area increased. In Lahore city, the area covered by water bodies, vegetation, and barren land was 2.7%, 24.9% and 42.5% in 1998, which decreased to 0.6%, 22.6% and 35.7% in 2018, while built-up land gradually increased from 29.8% in 1998 to 41% in 2018, due to the transformation of different land-use categories into built-up

Discussion
This research was carried out to quantify the individual contribution of different LULC transformations to LST in Lahore and Peshawar city (metropolitan areas of Pakistan) from 1998-2018, and future predictions of LULC for 2023 and 2028 by using CA-Markov model. Figure 3a highlights the land cover pattern for Lahore city and shows that the area covered by water bodies, vegetation, and barren land decreased gradually from 1998 to 2018, while built-up area increased. In Lahore city, the area covered by water bodies, vegetation, and barren land was 2.7%, 24.9% and 42.5% in 1998, which decreased to 0.6%, 22.6% and 35.7% in 2018, while built-up land gradually increased from 29.8% in 1998 to 41% in 2018, due to the transformation of different land-use categories into built-up areas. These conversions were mainly due to different housing and other projects by the Lahore Development Authority (LDA) and local developers [68,69], mega-development projects (including Metro Bus, Orange Train) [68,69,100].
However, Figure 3b shows, in Peshawar city, that the area covered by water, vegetation, and built-up land was about 1.8%, 25.0%, and 68.1%, respectively in 1998, and gradually increased to 2.4%, 50.6% and 21.5% in 2018. Barren land was found to decrease gradually with a higher rate (from 68.1% to 25.5%) due to its transformation into other land use classes. The increase in vegetation area was promoted by the "Billion trees project" (2013) of the KPK local government [72]. The increase in built-up area was needed because of bulk migration of people from Afghanistan toward the KPK Province in 2000-2005 during war-time [70].
Changes in the LULC also cause changes in the temperature [73,74], because different classes of LULC have different properties of reflectance and evapotranspiration [101]. Results of LST (Figure 7) indicate that in Lahore city during these 20 years, both the minimum and maximum LSTs are higher in 2018 than in 1998. These data confirm the results from earlier studies indicating that the LST rises in Lahore city [101,102], which may be due to the increasing use of artificial materials, such as asphalt and concrete, for urban expansion [103,104]. Our results from Figure 3a highlight that in Lahore city, the built-up area is increasing. The LST increases with the increase in urban built-up and barren land [105]. Meanwhile, in Peshawar city, the LST has decreased during the study period, for both the minimum and maximum values. Our results from Figure 3b highlight that in Peshawar city, the vegetation cover of the city has gradually increased in 2018 compared to 1998. The increase in vegetation area was promoted by the "Billion trees project" (2013) of the KPK local government [72]. Conversion between different LULC classes, especially the transformation of vegetation into built-up land, can effectively influence the land surface temperature [8][9][10].
Climatological data can be developed for two kinds of surface temperatures: near-surface air temperature Ta and the skin temperature, or land surface temperature (LST) [106]. Ta is measured at official weather stations at 1.5 m above the surface with sensors protected from radiation and adequately ventilated [107]. Hence, LST and air temperature are different parameters that cannot be directly compared but do influence each other. For the evaluation of the LST, they were compared with Ta from station data. The results in Figures 8 and 9 show a good comparison with high correlation coefficients between LST and Ta and confirm our LST trends.
Furthermore, the change in LST in response to the transformation between different land-use classes was studied by determining the relationships between LULC transition from 1998 to 2018 and LST using the Curve Fit tool and zonal statistic in ArcGIS (Section 2.6). The results indicate that the average change in LST ( • C) from 1998 to 2018 due to the transitions between the four different LULC classes is presented in Table 10 for each type of change: the transition from built-up and barren land to vegetation and water cover results in a decrease in LST, while the transition from vegetation and water to built-up and barren land results in an increase in LST. LST variation can be associated with various factors including LULC transition, population growth, and urban expansion [108], because the surface reflectance and roughness of different LULC classes are different, and different LULC types have unique properties in terms of the emission, reflection and absorption of radiation [11,12]. The future predictions results indicate that urban growth in Lahore city continues, at the expense of vegetation and barren land (Table 8); these changes will also impact the LST because urban expansion has a profound impact on local and regional climate [108,109].

Conclusions
Land use land cover (LULC) in Lahore and Peshawar, Pakistan, was determined using Landsat data for selected dates with 5 year intervals between 1998 and 2018. The changes in LULC classes during these 20 years and their impacts on urban thermal environments were determined. Markov and CA-Markov models were used to predict future changes for 2023 and 2028. The results indicate that much of the barren and vegetation land in Lahore has been converted into the urbanized area. Additionally, Peshawar city has experienced an increase in the built-up area, but vegetation cover and water have replaced barren land, which significantly changed the adverse effects of urbanization observed in Lahore to cooling in Peshawar. The different LULC types have unique properties in terms of the emission, absorption and reflection of radiation. Therefore, the gradual urban expansion and other anthropogenic activities affect the surface temperature, which is reflected by the analysis of the LST in response to the transition between LULC classes. Conversion from water or vegetation into built-up or barren lands results in the increase in LST. Conversely, the transition of areas from built-up or barren land into vegetation or water has helped to lower LST. The LST in Lahore city has substantially increased, while the LST in Peshawar city has decreased over the 20 year study period.
The future scenarios indicate that severe environmental problems may be expected in Lahore city, whereas the future scenarios of Peshawar show a way to improve living conditions in an expanding city. These results indicate that Lahore policies are needed to control current environmental issues and avoid further increase in these problems. This research provides the necessary information to environmentalists and policymakers in developing greening strategies for megacities to minimize future eco-environmental threats.