A remotely sensed tracking of forest cover and associated temperature change in Margalla hills

ABSTRACT

Worldwide, forest degradation is a serious environmental issue, and in Pakistan, forest wealth is depleting at the highest rate in South Asia. To ensure sustainable development goals of environmental stewardship, social development and economic growth, a sound monitoring and regulatory mechanism is essential for tracking forest cover changes. This study aims to quantify the decline of forest reserves and associated temperature variations in a relatively unexplored biodiversity hotspot of Islamabad, Margalla Hills National Park (MHNP). Imagery acquired by Landsat TM (Thematic Mapper) for the year 1992, 2000 and 2011 are used to assess the spatial and temporal changes occurred over the last two decades (from 1992 to 2011). A robust hybrid-classification routine is implemented to monitor the changes in forest cover and ANOVA along with Tukey’s HSD (Honestly Significant Difference) test is used to test the significance of temperature variation associated with a shift in land cover classes. The results showed a significant growth in settlements, agricultural area and barren soil whereas water body, lower vegetation, scrub and pine forest are diminishing. In both decades, the temperature alteration associated with a change in land cover classes are statistically significant (confirmed by ANOVA and Tukey’s HSD tests) for most of the land use/land cover classes. Based on these findings, this study concludes that forests are dwindling at MHNP and the degrading condition of the forest is below par and necessitates the promotion of conservation practices to minimize ecological disturbances.

KEYWORDS:

• Margalla hills
• Landsat thematic mapper (Landsat TM)
• forest cover change
• hybrid classification
• land cover classes
• ANOVA
• Tukey’s HSD

1. Introduction

Forests as renewable natural resources are the most potential sources of wealth that provide various economic, social and environmental benefits (Haque, Murty, and Shyamsundar 2011). Forests host most of the earth’s animal and plant species and play a crucial role in the balancing of an ecosystem (Peh et al. 2014). Forests ecosystem provide food, medicines to benefit humankind and is a major sink of carbon dioxide (CO₂) emissions, thus plays an essential role in regulating global weather (Chiesi et al. 2005; UNFCCC 2009; Baccini et al. 2012; Chen et al. 2016). According to UNDP, at least 25% of the land should be under forest for a healthy balance economy of a country (FAO 2012). Because of its prime importance, the policy-makers are enforcing policies and conventions to protect forests, sustain life and earth environment.

The world’s forest covering areas are shrinking at an alarming rate because of unsustainable logging and clearing ground for agriculture and settlement purposes (Hansen et al. 2013; IPCC 2013). Every year, the world is losing about 13–15 million hectares of forests due to deforestation (Chakravarty et al. 2012; FAO 2012). Scientists project that if deforestation continues at the current alarming rate, in the next 100 years, the entire world forests will be finished. The decline in forests is considered the major reason for climate change as it is a major natural sink, second only to oceans in their ability to minimize atmospheric carbon dioxide (CO₂) by sequestering it in photosynthesis (Klein and Perkins 1987; Allen et al. 2010; Ghasemi, Sahebi, and Mohammadzadeh 2013). Large-scale conversion of forested land to other land covers results in climate warming (Nobre, Sellers, and Shukla 1991). Scientists report that the forests can store up to 1.5 metric tons of carbon per hectare per year (Gough et al. 2008). The world forests lock more than 1 × 10¹⁵ metric tons of carbon in biomass, detritus and soil (Dixon et al. 1994). Deforestation annually releases 1400–2000 teragrams of carbon into the atmosphere which equals 17–25% of global fossil fuel emissions every year (Ryan et al. 2007). The hike in CO₂ concentration has led to the global temperature rise and has threatened the world environment (Fearnside 1996; Hansen et al. 2006). Intergovernmental Panel on Climate Change (IPCC) estimates a 5.2–9.5°F (2.9–5.3°C) increase in the Earth’s average surface temperature by the year 2100 (Nakicenovic and Swart 2000). This alarming rate has stirred concern among scientists world over to conduct researches on the decline of forest reserves, and have driven governments, regulative agencies and environmental protection authorities to frame regional and global policies (Chakravarty et al. 2012; FAO 2012; IPCC 2013).

In Pakistan, the forest covered area is not more than 5% of the total land (IUCN 2005; Shahbaz, Ali, and Suleri 2007). Majority of the country forests are located in Khyber Pakhtunkhwa (a province in the north and western part) which accounts for over 40% of national forest wealth and is continuously dropping over the years (IUCN 2005). Pakistan is reported to have the highest rate of deforestation in South Asia (World Bank 2013). This alarming situation demands detailed research efforts to quantify the decline of forest reserves. Various researchers have attempted to study forest attributes using old conventional techniques which are based on human observations and have rarely employed latest technologies like satellite remote sensing. As such, conventional methods are costly, more time consuming, demand large human resources and lack continuous spatial coverage, which hardly help the government in enforcing regulations and policies to mitigate deforestation. Remote sensing (especially satellite remote sensing) offers an effective alternative to monitor forests wealth at a low cost and more precisely to quantify forest cover changes at different scales. Moreover, using remote sensing technology, results can be reproduced. Remote sensing applications in forestry not only reduce the survey crew and cost, but can observe areas which are otherwise inaccessible by conventional approaches. Thus, remote sensing has emerged as a vital decision-support tool for forest managers to monitor forests over large areas at reparative intervals.

Several researches clearly validate the potential of integrating high resolution multi spectral remote sensing data and aerial photography for observing land pattern changes (Green, Kempka, and Lackey 1994; Shalaby and Tateishi 2007; Hansen et al. 2008). By means of remotely sensed data, forest structure has been effectively modeled in the state of Nepal over a span of 14 years (Gautam, Webb, and Eiumnoh 2002). In Changbai Mountain range of China, changes made in forest cover due to management-induced activities have been successfully investigated using Landsat TM images (Li et al. 2006). Quantifying and exploring vulnerability in forests of Kenya indicated their decline with a subsequent increase in agricultural parcels (Wachiye, Kuria, and Musiega 2013). Choropleth maps of the forested area in Pakistan have been formulated which highlight the dwindling coniferous forest over the last two decades (Ahmad et al. 2012).

For the efficient management of forest resources, and in particular, deforestation, forest cover change must be monitored at both the regional and global scale. National inventories and international treaties (United Nations Framework Convention on Climate Change (UNFCCC), REDD+, Kyoto Protocol, etc.) also demand an accurate assessment of forest attributes at large scale which is possible through the analysis of remote sensing imagery. Hence, this study aims to monitor forest pattern on Margalla hills from 1992 to 2011 using visible to short-wave-infrared (VSWIR) data of Landsat TM sensor. Besides tracking land cover changes, the impact of class shifts on surface temperature is also assessed using thermal infrared (TIR) band of Landsat TM.

2. Materials and methods

2.1. Study area

The study area is comprised of Margalla Hills and the adjacent Rawal Lake regions (Figure 1) of Islamabad, Pakistan. Margalla Hills National Park (MHNP: covers an area of 173.9 Sq. Km) is a hill range immediate north of Islamabad capital territory and is situated at the foothills of lesser Himalayan. MHNP, forming the northern boundary of the Potohar plateau is located between 33°40′01″–33°42′43″ N, 72°45′01″–72°52′32″E and is at an elevation varying from 1347 to 3907 feet above the mean sea level. MHNP was declared as a National Park in 1980 under Islamabad Wild Life Ordinance (Malik, Husain, and Benes 2003). The landscape is furrowed, rocky and mainly composed of limestone (Nasir and Akhter 1987). The weather is harsh in nature with hot summers and cold winters. In summer season (May and June), the temperature escalates up to 42°C, whereas in winter (months of December and January) the temperature falls below zero (Hussain 1986; Khalid et al. 2015). The Margalla Hills occurs at the far end of monsoon zone and the precipitation mostly falls in July and August (in these months the mean rainfall is 254 mm). The mean relative humidity varies between 59% and 67% during monsoon period (Masroor 2011). The area is rich in plant biodiversity and hosts 616 species of flora. The dominant ecological zones are subtropical coniferous pine (evergreen) and subtropical deciduous scrub (deciduous) forests. The common plant species found in Margalla hills are Chir (Pinus roxburghii), Kao or Kahu (Olea ferruginea), Sanatha (Dodonaea viscosa), Phulai (Acacia modesta), Toot (Morus alba) and Shisham (Dalbergia sissoo). MHNP also provides shelter to a large number of animal biodiversity including 30 mammals, 27 fish, 21 amphibians and reptiles, 238 birds, 39 butterflies and several other invertebrate species (WWF 2009). MHNP inhabits 100,000 people dwelling in 37 small villages and hamlets. About 60% of the population residing in MHNP premises is practicing rain-irrigated agriculture and domesticating animals for milk and meat supply. Besides Sorghum, corn and wheat are the main crops of the area (WWF 2009).

Figure 1. Study area map showing Margalla Hills National Park and Rawal Lake.

2.2. Dataset and preprocessing

The datasets comprised of multispectral Landsat TM imagery for the year 1992, 2000 and 2011. The VSWIR bands of Landsat TM were used to categorize various land cover classes. Landsat TM TIR band for the same dates were acquired to compute land surface temperature (LST). The satellite imagery acquired was from a nearly identical time of the year (i.e. from 21 May 1992, 31 May 2000 and 5 June 2011). Besides satellite images, ground truth and ancillary data of maximum, minimum and average temperature were also acquired from the field and the Pakistan Meteorological Department (PMD), respectively.

Atmospheric correction was performed to eliminate atmospheric effects (such as Rayleigh scattering, Mie scattering, etc. which is the dominant atmospheric affect in Landsat TM) (Song et al. 2001) and the retrieval of normalized reflectance across different bands. This is necessary to carry out quantitative analysis and change detection with multi-temporal or multi-sensor data (Brunn et al. 2003).

Pixels of Landsat TM Level-1 data expressed as digital numbers (DNs) were radiometrically calibrated to retrieve absolute at-sensor Radiance (W m⁻² sr⁻¹ μm⁻¹) using gain and bias values (see Equation (1)).(1) $L_{\lambda }=\text{gain}\ast \text{DN}+\text{bias},$(1)

where, ‘L_λ’ is the radiance pixel value, ‘DN’ is the digital number cell value, ‘gain’ is the gain value for a specific band and ‘bias’ is the bias value for a specific band.

Radiance was then converted to atmospherically corrected spectral Reflectance using the FLAASH (Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes) algorithm. Standard radiance rescaling factors for each band provided in the metadata were used. With initial visibility set to 40 km, a 2 band Kaufman-Tanre Aerosol retrieval method was used on a Tropospheric Aerosol model and a Tropical atmospheric model. Multispectral sensors (like Landsat TM) do not have appropriate bands to perform water retrieval but the average water amount was set to 1.096 atm-cm. The cirrus cloud masking was performed using wavelength range 1.55–1.75 µm.

The VSWIR bands of satellite images were then stacked and re-projected in Universal Transverse Mercator (UTM) and World Geodetic System (WGS) 1984. The study area (MNHP and Rawal Lake) were extracted from the images by subsetting the images using shapefile of the region of interest. The TIR bands were also cropped according to the size of the study area for the year 1992, 2000 and 2011.

2.3. Categorization of land cover classes and change detection

A classification procedure called hybrid maximum-likelihood routine was adopted to categorize the study area into various land cover classes and to assess the shift in forest cover. The hybrid classification encompasses the advantages of both supervised and unsupervised classification and is considered a robust routine (Kuemmerle et al. 2006; Li et al. 2013; Zhang et al. 2014). First, the subset images were subjected to unsupervised classification for spectral clustering and self-organization. In the spectral clustering, if a class of significant size was absent after unsupervised classification, a supervised routine was followed to add their signature. All the spectral signatures were then combined to form one signature file. The final signature file (ensuing both unsupervised and supervised classification) was then used to classify the entire image in different land use and land cover classes such as agricultural land, barren soil, coniferous forest, scrub forest, built-up area, lower vegetation and water body. The evaluation of classification results is an essential process in the procedure of satellite image classification to find the degree of accuracy. In doing so error/confusion matrices were used. This is the frequently employed approach for evaluating per-pixel classification (Lu and Weng 2007). To assess the accuracy of results, 120 random points (reference data) were created on the unclassified subset image. Confusion matrix was constructed for the three years and the points were given reference on the basis of best guess of their representative class. Classification was followed by post-classification procedures like vectorization, class extraction, merging and area calculation. To quantify the shift that has occurred in various classes, land cover maps for the three years (1992, 2000 and 2011) were overlaid and changed maps were generated to quantify the shift between the years from 1992 to 2000 and then 2000 to 2011.

2.4. ANOVA and Tukey’s HSD

The statistical tests such as one-way analysis of variance (ANOVA) followed by Tukey’s HSD (Honest Significance Difference) were applied to the LST data retrieved from TIR band of Landsat TM. ANOVA was used to determine whether the temperature variance among various group means differs statistically. The null hypothesis stated that the variance of temperature means for all changed classes is the same, i.e. there is no significant difference between them. ANOVA was tested at confidence level of 99% (p < .01). Once it was established that spatial temperature means were statistically different for various classes, post hoc Tukey’s HSD tests were implemented for pairwise comparison to ascertain that altered (shift) classes were statistically different from their ancestor class (parent class). Pairwise comparisons for temperature were first made for the year 1992 and 2000 and then 2000 and 2011. The schematic flow diagram of the methodology followed for the research is shown in Figure 2.

Figure 2. Schematic flow of the implemented methodology. Preprocessing included band stacking, reprojection and subsetting. And post-classification processing included reclassification (or recoding), vectorization of classes, class extraction, merging and geometry computation.

3. Results

3.1. Multi-temporal land covers classification

Based on hybrid classification of sensor data for the years 1992, 2000 and 2011, Table 1 specifies the area (percentages and hectares) of various land use and land covers in thorough detail. In the year 1992, the lower vegetation was the dominant class enveloping an area of 5700 ha (equal to 29% of the total MHNP). This class encompassed range lands, shrubs, herbs and grass lands. Coniferous pine (composed of Pinus roxburghii) forest was the second major class having total area coverage of 5007 ha (occupying 25% of the total study site). Water bodies covered the rest of the area (594 ha) which was 3% of the total land (Figure 3). In the year 2000, lower vegetation was again the dominating land cover class covering an area of 4398 ha (22% of the total area). The second major class was settlements with an area of 3524 ha (18% of the study area). The spatial extent of Water was 316 ha (2%-of the study area) and remained the lowest class in the year 2000. In 2011, lower vegetation covered an area of 4131 ha (21% of the total surface of MHNP) and was the major class. Next major classes were the settlements and the scrub forest covering 19% and 17% of land, respectively (Figure 4– 6).

Figure 3. Pie-charts showing percentage of area coverage of land cover classes on MHNP in 1992, 2000 and 2011.

Figure 4. Land use and Land cover map of Margalla Hills National Park for the year 1992.

Figure 5. Land use and Land cover map of Margalla Hills National Park for the year 2000.

Figure 6. Land use and Land cover map of Margalla Hills National Park for the year 2011.

Table 1. The area coverage (both in hectare and percentage) of different land cover classes in MHNP for the year 1992, 2000 and 2011.

Land cover classes	1992		2000		2011		Percentage change
	Area (Ha)	Area (%)	Area (Ha)	Area (%)	Area (Ha)	Area (%)	1992–2000 (%)	2000–2011 (%)
Agriculture	2124	11	2529	13	2613	14	9	4
Settlements	1921	10	3524	18	3656	18	83	4
Lower vegetation	5700	29	4398	22	4131	21	−22	−6
Pine forest	5007	25	2336	12	2980	15	−53	27
Scrub forest	3167	16	3489	17	3269	17	10	−6
Barren soil	1199	6	3120	16	2752	13	160	−11
Water body	594	3	316	2	311	2	−46	−1.5
Total	19,712	100	19,712	100	19,712	100

Notes: The last column specifies the percentage change in area of classes from 1992 to 2000 and 2000 to 2011. Percentage change = [(b−a)/a]×100, where a is the area of the class belonging to the erstwhile year and b is the area of the class from the more recent year. The positive sign shows an increase while the negative sign shows decrease in different years.

The changing trends in land use and land cover classes over the years are illustrated in Figure 7. The built-up region and agricultural land consistently increased. While the graph of the water body and lower vegetation reveals a continuous decreasing trend from 1992 to 2011, for other LULC classes such as barren soil, scrub forest and pine forest, the pattern did not transform with uniformity from 1992 to 2000 and then from 2000 to 2011.

Figure 7. LULC changes in Margalla Hills from the year 1992–2011.

3.2. Temporal shifts in the LULC classes

The post-classification comparison between the years revealed spatial and temporal shifts among classes and also pointed out classes that are relatively stable over time (see Figures 8 and 9). The major LULC classes which the forests transformed into were the agricultural land, barren soil and built-up areas. However, the drop in lower vegetation class contributed to an increase in forest class (i.e. lower vegetation shifted mainly to pine and scrub forest).

Figure 8. Post-classification comparison between the classified images of 1992 and 2000 presenting the shift among classes.

Figure 9. Post-classification comparison between the classified images of 2000 and 2011 presenting the shift among classes.

Considering the entire area covered by all the classes of MHNP in 1992, a total of only 22% remained unchanged in 2000 while 77% of the area was shifted to other classes (see Table 4). This shift that occurred was a combined effect from all classes, i.e. 11.2% scrub forest area transformed while 21% of the change was contributed by Pine forest, etc. However, in the year 2011 in MHNP, the change occurred was 49%. An area of 10.2% out of the total land continued to exist as pine forest in the year 2011, while 1.78% transformed to other classes. Estimations revealed that 2.2% of Scrub forest persisted as so, while 20.4% shifted to other feature classes (see Table 4).

3.3. Accuracy assessment and kappa statistics

The assessment of classification accuracy of the images for the years 1992, 2000 and 2011 was carried out to determine the quality of information derived from the data. The results attained through this method comprised of the overall classification accuracy report, Kappa statistics (Table 2) and Cell report. The accuracy results indicated a high rate of classification accuracy that was 99.17%, 98.33% and 98.33% again for the years 1992, 2000 and 2011, respectively. Kappa accounts for all elements of the confusion matrix and excludes the agreement that occurs by chance. Consequently, it provides a more rigorous assessment of classification accuracy. The overall accuracy and Kappa statistics achieved for the classification of the images for 1992, 2000 and 2011 are reasonably acceptable in order to conduct the analysis and change detection.

Table 2. Overall classification accuracy and kappa statistics.

	1992	2000	2011
Overall classification accuracy (%)	99.17	98.33	98.33
Overall Kappa statistics	0.9843	0.9627	0.9679

3.4. Change in temperature associated with class transition

The results of ANOVA (one-way analysis of variance) shows that the variances of the mean temperature of various classes were statistically significant (p = .000) indicating that the mean temperature of at least one group was statistically different (see Table 3). Once the null hypothesis (i.e. classes mean temperatures do not differs from each other) was rejected, post hoc Tukey’s HSD was used for pairwise comparisons. The Tukey’s HSD reveals that the mean temperatures of changed classes are statistically different from their ancestor class in most cases between the years 1992–2000 and 2000–2011 (see Figures 10 and 11).

Figure 10. Boxplots for temperature difference and graph of P-values for the change classes from 1992 to 2000.

Figure 11. Boxplots for temperature difference and graph of P-values for the change classes from 2000 to 2011.

Table 3. Summary of ANOVA of the changed classes land surface temperature in the years 1992–2000 and 2000–2011.

Years	Degree of freedom	F-value	p-Value
1992–2000	36	7957	0.000
2000–2011	17	373.8	0.000

The boxplots (Figure 10) reveals that from the year 1992 to 2000, the transition in classes such as Water to Agriculture (W to A), Scrub forest to Soil (S to So) show a wider range of temperature differences. Some classes such as Settlements to Water (Se to W), Soil to Pine (So to P) and Soil to Water (So to W) show very little range of temperature difference. For the majority of the classes, the P-values resulting from the Tukey’s HSD were statisticaly significant (p < 0.05).This means that for these classes the temperature change is statistically significant with the transition in cover classes. The pairwise comparison of few classes such as transitions from Soil to Pine (So to P) and Water to Pine (W to P) rendered p-values higher than 0.05, which means that the mean temperature difference of these transitions classes is not statistically significantly different from their parent/ancestor classes. Low vegetation to Water (L to W), Settlements to Water (Se to W), Soil to Water (So to W) and Water to Scrub (W to S) forest also yielded relatively higher p-values.

For the year 2000–2011, the boxplots show that some transitions such as Low vegetation to Agriculture (L to A), Scrub to Pine (S to P) and Low vegetation to Settlements (L to Se) cover a wide range of temperature differences (see Figure 11). Other classes such as Settlements to Pine (Se to P) and Soil to Low Vegetation (So to L) had minimum ranges of temperature difference. The p-values for majority of the classes were statiscally significant (Figure 11) except a few, such as Soil to Scrub (So to S), Water to Settlements (W to Se) and Low vegetation to Agriculture (L to A).

4. Discussion

4.1. Transitions in cover classes

In Pakistan, the forests cover is over an area of nearly 4.6 million hectares which is almost 2% of the total land (Rich 2013). The forest cover is continuously dropping over the years and the country is regarded to have a very high percentage of annual deforestation as a function of the forest area (Allen and Barnes 1985). Since the creation of Pakistan (1947), about 61,000 ha of the valuable forested land has been deforested. There are many factors responsible for the decline of forest resources, one of the principal reasons influencing it is a competitive and flourishing world economy. Although the forest reserves of Pakistan are meager but they contribute considerably to the economy (Triana et al. 2013). Margalla Hills provide fodder, timber, non-timber forest products (NTFPs) and wood for fuel (Haque, Murty, and Shyamsundar 2011). The findings specified in the cross-tabulation (see Tables 4 and 5) summarized the changes and shifts in land cover classes in a more simplistic way. In Margalla Hills, the major transitions in land cover classes happened between the years 1992–2000. During this period, 77% of the total land of transformed to other classes while only 22% of the area remained unchanged. From 2000 to 2011, the proportion of changed and unchanged area was 49% and 51%, respectively.

Table 4. Summary of the class shift from the year 1992 to 2000.

Class (1992)	% Unchanged	% Changed	Changed Class in 2000 in Hectares (ha)							Total (1992)
			Agriculture	Water	Soil	Pine forest	Scrub forest	Settlements	Lower vegetation
Water	1.5	1	47	305	89	2.3	2.3	41	32	518.6
Settlements	3.5	3	43	0.15	309	7	84	700	195	1338.15
Low vegetation	3	25	1332	1.3	119	1148	1811	650	576	5637
Agriculture	0.2	10	41	–	1015	4	71	669	307	2107
Soil	4.4	5.8	8	0.15	872	1	14	228	899	2022.15
Scrub forest	4.6	11.2	294	–	292	477	910	454	702	3129
Pine forest	3.4	21	736	–	338	671	1560	660	996	4961
Total (2000)	22	77	2501	306	3034	2310	4452.3	3402	3707	19712

Table 5. Summary of the class shift from the year 2000 to 2011

Class (2000)	% Unchanged	% Changed	Changed class in 2011 in hectares (ha)							Total (2000)
			Agriculture	Water	Soil	Pine forest	Scrub forest	Settlements	Low vegetation
Water	1.5	0.7	20	305	12.54	–	–	1	8	346.54
Settlements	10.4	13.9	–	–	63	0.3	2649	2058	30	4800
Low vegetation	3.9	3.96	117	–	0.3	–	743	22	776	1658.3
Agriculture	9.89	2.8	1951	–	–	424	–	1.68	133	2509.7
Soil	13.4	4.68	–	–	2652	–	0.3	926	0.22	3578.5
Scrub forest	2.2	20.4	120	–	–	725	438	13	3167	4463
Pine forest	10.2	1.78	309	–	43	2012	––	–	–	2364
Total (2011)	51	49	2517	305	2770.8	3161	3830	3021	4114.2	19712

Scrubs forest is one of the dominating land cover class in MHNP. Some scrub forests are perennial (evergreen) while others are deciduous or semi-evergreen. Margalla Hills mainly populate the deciduous type of scrubs which shed their foliage and thrive well in tropical wet and dry climate. These forests shed leaves and are adaptive to the scarcity of water in long dry seasons (Khalid et al. 2015). In this study, the scrub forest initially showed an increase of 10% from the year 1992 to 2000 followed by a decrease of 6% in the next decade (from 2000 to 2011). Since the year 1992, 13.9% of the land occupied by scrub forest transformed to other classes while 10.4% sustained as previous. The clearance of scrub forest could be attributed to the spread of settlements, timber mafia, population boost and urbanization. A prevailing practice in MHNP is the removal of strips of forest to grow crops by farmers (Khalid et al. 2015). The periodic wild forest fires which turn down the vegetation to ashes could be another cause possible for the reduction in forests of MHNP (PARC 2009).

Subtropical pine forest is another major class in the hills. The chir pine (Pinus roxburghii) is the most abundant species growing at an elevation of 1000 m above the sea level. The pine forest class declined by 53% between the years 1992 and 2000. In the year 1992, the pine forest occupied a total land of 5007 ha that colossally declined to 2336 hectares in the year 2000. Within the premises of MHNP, 21% of the total pine forest changed to other land covers while 3.4% continued as pine forest class. The forest may have been influenced by prolonged drought, inclement weather and climate conditions, when there is no downpour for months and temperature soars up to 45°C (Sykes and Prentice 1996; Zonneveld et al. 2009; Lindner et al. 2010). Extensive cum recurrent cattle grazing and fuel wood retrieval by the local communities have warped the vascular plants (higher plants) to non-vascular plants or bushes (Shafiq et al. 1997). The clearing of valuable pine forest for agriculture practice could also be the reason of decline in pine forests (Khalid et al. 2015). Cement industries exploit Margalla Hills for the acquisition of limestone and shale and hence, contribute to deforestation (Khan 2003; Iqbal, Khan, and Malik 2013). The outcome of this study is in line with the findings of Sykes and Prentice (1996) and Lindner et al. (2010). Between 2000 and 2011, an increasing trend (27% increase) was observed in the pine forest. Estimations let out that 10.2% maintained as so, while 1.78% shifted to other classes. The increase in pine forest acreage between 2000 and 2011 could be accredited to the better management activities such as ban on cutting pine trees for fuel wood and timber purposes. The plantation campaign of pine forest in MHNP could also be a reason for the increase in pine forest area.

Over a period of two decades, a total of 13% increase was observed in agriculture area in Margalla Hills. Initially, there was 9% increase between the years 1992 and 2000, whereas after a period of 11 years (from 2000 to 2011), the agriculture area increased by 4%. About 60% of the population living in MHNP is dependent on agriculture for their livelihood (Ayaz 2005). The inhabitants are involved in agricultural practices and herding their cattle in the park area. Expansion of road network is one of the many reasons for the burgeoning agriculture in the study site. Roads provide farmers an easy access to encroach and fragment the forests for growing crops (Khalid et al. 2015).

With increasing human population and migration towards the cities, the demand for new settlements has risen in Islamabad and its surrounding areas. This unplanned urbanization has adversely affected the ecosystem of MHNP. Just in one decade (1992–2000), the settlements class increased by 83%, from 1921 ha to 3524 ha. In the next decade (2000–2011), the settlements class increased further by 4% (reached to an area of 3656 hectares). The surging built-up area could be associated with the ever-increasing population and conurbation, which led to the development of housing projects and residential colonies to accommodate them. Other setups include road networks, recreational pursuits (i.e. Daman-e-Koh, Lake View Park, PirSuhawa and Said-pur village), industrial and commercial units. These developments happened at the cost of forest decline, shrinking water bodies and increase in barren soil.

Lower vegetation in MHNP shows a decline of 22% from the year 1992 to 2000. In 2011, it further dropped by 6%. The locals are relying on MHNP for herding their livestock; hence, overgrazing could be a major reason behind the decline in lower vegetation. Wild fires also cause rapid desertification of range lands to the barren soil. Apart from overgrazing and forest fires, urban sprawl and construction of settlements are other factors causing the decline in lower vegetation acreage.

The 160% upsurge in Barren soil could be linked to anthropogenic activities like haphazard bush burning, stone crushing, over grazing of pastures, fuel wood collection and wild fires. Much of the abandoned land was then transformed into settlements after year 2000. The class thus exhibited a rise of 11% from 2000 to 2011.

In the year 1992, the water body (Rawal Lake) occupied an area 594 ha which shrank to 316 ha in the year 2000. The class further dropped to 311 ha in the year 2011. Rawal Lake is the chief supplier of water to the twin cities, supplying 23 Million Gallons per Day (MGD) for drinking and domestic usage to Rawalpindi city, and a limited amount to Islamabad capital for irrigational purposes (Hashmi, Qaiser, and Farooq 2012). The regime of water in Rawal Lake has altered because of the deforestation and rapid urban expansion in the catchment areas. The recharging capacity of the aquifer in Islamabad and Rawalpindi is adversely affected by the formation of road network, pavements and other structures. The concrete surfaces restrict infiltration that eventually affects the surface and subsurface hydrology of the region. The storage capacity of Rawal Lake that used to be 47,230 ha in the year 1960 has been dropped to approximately 35%, owing to the process of sedimentation (IUCN 2005).

4.2. Temperature change associated with the class shift

Different land cover classes have a different temperature, due to evapotranspiration, forests and vegetated areas are always cooler than the surrounding barren soil and rocky grounds. The temperature drift associated with temporal transitions of land cover classes are tested by analyzing the spatial variation of LST derived from a thermal band of Landsat TM. The temporal images of TIR belong to the nearly identical time of the year from the years 1992,2000 and 2011 (Table 6). LST was derived from 21 May 1992, 31 May 2000 and 5 June 2011. It is assumed that the temperature remained similar during this time. The air temperature data recorded by PMD confirmed that the average air temperature at the time of image acquisition was nearly equal (see Table 6). So, it is hypothesized that if there is any change in class temperature that would be associated with the class shift over a period of two decades.

Table 6. Summary of air temperature recorded for the month at which the images were acquired.

Year	Date of image acquisition	Minimum monthly temperature (°C)	Maximun monthly temperature (°C)	Average monthly temperture (°C)
1992	21 May (day 141)	19.5	31.0	26.3
2000	31 May (day 151)	17.35	30.1	26.0
2011	5 June (day 156)	19.7	31.7	26.0

The majority of changed classes show a considerable variation in temperature (see boxplot of Figures 10 and 11) and the difference in mean temperature are statistically significant (see p-value plot of Figures 10 and 11). The changed classes such as Agriculture to Low vegetation (A to L), Agriculture to Pine (A to P), Low vegetation to Scrub (L to S), Pine to Agriculture (P to A), Scrub to Agriculture (S to A), Scrub to Low vegetation (S to L), Scrub to Pine (S to P), Scrub to Settlements (S to Se), Settlements to Low vegetation (Se to L), Settlements to Scrub(Se to S, Settlements to Soil (Se to So), Soil to Settlements (So to Se) all yielded results that were statistically significant (p < .05) in both decades. This reveals that the temperature change is statistically significantly related to the changed classes. However, some classes in 1992–2000 such as Settlements to Water (Se to W), Soil to Pine (So to P) and Soil to Water (So to W) and a few classes in 2000–2011 such as Settlements to Pine (Se to P) and Soil to Low Vegetation (So to L) showed very little range of temperature differences because of less spatial extent covered by the classes, i.e. the areas covered by the shifted classes were too small to give statistically significant temperature change with the classes. The switch of vegetation/forest cover class to other land cover types yielded a drift in temperature which is significantly different from their parent class. Previous research also found that deforestation and clearing of vegetation lead to an increase in temperature and dust storms (Ellis, Taylor, and Masood 1994) which is in line with the finding of this study. It has been established that changes in temperature and soil moisture is a critical determinant of vegetation spatial distribution (Lookingbill and Urban 2003). Sea and Whitlock (1995) have developed an inventory of the plant species and climatic history of a region, which suggested that vegetation flux has been influenced greatly by temperature variability over the last 14,000 years. Research focusing on the relationship between LST and vegetation abundance shows that the surface radiant temperature response is a function of changes in soil moisture content and vegetation cover. LST increases in sparsely vegetated areas such as urban areas (Weng, Lu, and Schubring 2004).

5. Conclusion

Earth observation data provides a sound monitoring and regulatory mechanism for tracking forest cover changes and associated temperature variations. The present study is conducted in a relatively unexplored biodiversity hotspot and reveals that land cover and land use practices have altered in MHNP since the last two decades. The transitions in land use and land cover classes have influenced the LST and are significantly different from their ancestor cover class. Based on the spatio-temporal analysis of earth observation data, this study concluded that the forest reserves are dwindling, whereas built-up and agriculture land is growing. The expansion of agriculture and settlements has subsequently affected the forest cover and water bodies. The findings of this study highlight the needs to conserve forest resources by controlling human interference in MHNP and its surrounding areas. Land use/cover must be monitored using continuous spatially explicit earth observation data and forest restoration activities should be initiated wherever needed. To ensure sustainable forest management, proper planning for urbanization is essential and areas with less ecological importance must be utilized for domestic, industrial and commercial settlements.

Acknowledgement

We thankfully acknowledge Pakistan Meteorological Department (PMD) for providing us data, World Wildlife Fund (WWF) for sending us their report on MHNP and Environment department of Capital Development Authority (CDA) for providing information vital to conduct this research.

Disclosure statement

No potential conflict of interest was reported by the authors.