Monitoring of Urban Landscape Ecology Dynamics of Islamabad Capital Territory (ICT), Pakistan, Over Four Decades (1976–2016)

Abstract

In the late 1960s, the Islamic Republic of Pakistan’s capital shifted from Karachi to Islamabad, officially named Islamabad Capital Territory (ICT). In this aspect, the ICT is a young city, but undergoing rapid expansion and urbanization, especially in the last two decades. This study reports the measurement and characterization of ICT land cover change dynamics using Landsat satellite imagery for the years 1976, 1990, 2000, 2010, and 2016. Annual rate of change, landscape metrics, and urban forest fragmentation spatiotemporal analyses have been carried out, along with the calculation of the United Nations Sustainable Development Goal (SDG) indicator 11.3.1 Land Consumption Rate to the Population Growth Rate (LCRPGR). The results show consistent increase in the settlement class, with highest annual rate of 8.79% during 2000–2010. Tree cover >40% and <40% canopy decreased at an annual rate of 0.81% and 0.77% between 1976 to 2016, respectively. Forest fragmentation analysis reveals that ‘core forests of >500 acres’ class decreased from 392 km² (65.41%) to 241 km² (55%), and ‘patch forest’ class increased from 15 km² (2.46%) to 20 km² (4.54%), from 1976 to 2016. The LCRPGR ratio was 0.62 from 1976 to 2000, increasing to 1.36 from 2000 to 2016.

1. Introduction

The conversion of one land cover type to another is one of the most visible and rapid changes that the earth is experiencing and these conversions have profound social and environmental impacts at multiple scales. In the past, the drivers of land cover changes have been classified into two categories: proximate and distant, or indirect [1]. Recently, the phenomenon of tele-coupling between places has been recognized as key to understanding how distant drivers relate to proximate drivers and influence local landscape changes [2,3,4]. Some of the key drivers of landscape changes include economic, technological, institutional and policy, cultural, and demographic factors [5].

The growth of urban areas leads to land cover change in many parts of the world, especially in developing countries [3,6]. Intense urbanization and increase in anthropogenic activities reflect the scope, intensity, and frequency of human interference, and the changes they cause in ecological processes and systems in the urbanized areas [7]. The urban areas consist of only 1–6% of the earth’s land surface, yet they have enormous impacts on the functioning and service of local and global ecosystems, by modifying local climate conditions, eliminating and fragmenting native habitats, generating anthropogenic pollutants, etc. [8]. The spatial pattern of an urban landscape is a result of the interaction between various driving forces including natural and socioeconomic factors [9]. Increasing trends in industrialization and urbanization, along with the migration from rural to urban areas, are the most dominant factors influencing the land cover transformation. In the rural areas, employment opportunities and income are insufficient, which contributes to large differences in income and facility levels between urban and rural areas [10]. In developing countries, new cities are being developed due to human migration, infrastructure development, and growing job opportunities [10,11].

The United Nations (UN) Sustainable Development Goals (SDGs) are a collection of 17 global goals, which include 232 indicators set by the UN General Assembly for the year 2030. These goals are an urgent call for action by all countries—developed and developing—in a global partnership. Pakistan is a signatory to the UN SDGs, and this study analyzes the land use change dynamics of Islamabad Capital Territory (ICT) in the context of the relevant SDGs and indicators. The Goal 11 of the SDGs, “Make cities and human settlements inclusive, safe, resilient and sustainable”, with indicator number 11.3.1 “Ratio of Land Consumption Rate to Population Growth Rate (LCRPGR)” [12] is an important parameter for analyzing sustainability of land use and land change with population growth. The LCRPGR parameter is actually based on the previously defined parameters of LCR [13,14] and PGR [15], and several studies of urban areas have been carried out before on the basis of these parameters. The LCRPGR parameter and terminology has been given prime importance now after being linked with the UN defined SDGs. LCRPGR is vital to understand the rate of land change as compared to the population boom, to understand historical land consumption traditions, and to guide decision and policy makers on the planned expansion of the city along with the protection of environmental, social, and economic assets. Another SDG, Goal 15, “Protect, restore and promote sustainable use of terrestrial ecosystems, sustainably manage forests, combat desertification, and halt and reverse land degradation and halt biodiversity loss” emphasizes the protection of tree cover and sustainable management of forests. The SDG indicator number 15.2.1, "Progress towards sustainable forest management" [16] calls for sustainable management strategies for conserving the forest cover, enhancing environmental education, and engaging a wide range of stakeholder institutions, policies, regulations, and considerations that promote sustainability and utilization of natural resources at multiple spatial scales [17].

Assessment and monitoring of land cover dynamics is essential for the sustainable management of natural resources, environmental protection, biodiversity conservation, and developing sustainable livelihoods. Therefore, the development of applicable and systematic methods for producing and updating land cover databases are considered an urgent need [18]. Around the globe, urban land expansion rates are higher than or equal to urban population growth rates [19,20,21]. Many research studies focus on big cities and metropolises, where increase in population through analysis of census statistics is directly linked with the urban land expansion; however, these statistics do not provide information regarding spatial distribution, pattern, and scale of urban land use change. Multi-temporal land cover change analysis and simulation based on coarse to very high resolution satellite remote sensing images is becoming a well established technique for quantifying changes occurring on the earth’s surface, and multi-temporal aerial and satellite datasets are now widely and continuously being used for urban growth mapping, monitoring, and modeling with a focus on the spatial dimension and structures [22,23,24,25]. In urban expansion studies, spatiotemporal analysis of land cover and land use changes has helped towards understanding the underlying natural and socio-economic factors and drivers. For instance, Seto et al. [26] have presented a meta-analysis of 326 studies which used temporal satellite images to map urban land conversion. A total of 58,000 km² increase in urban land area was reported in thirty years (1970 to 2000) and by 2030, global urban land cover is expected to increase between 430,000 km² and 12,568,000 km², with an estimate of 1,527,000 km² more likely. According to Seto et al. [26], across all regions and for all three decades, urban land expansion rates are higher than or equal to urban population growth rates. Yang et al. [27] studied and reported the evidence of urban agglomerations through satellite images in four major bay areas of US (San Francisco and New York), China (Hong Kong-Macau), and Japan (Tokyo), from 1987 to 2017.

Clarke et al. [28] proposed a framework to combine remote sensing and spatial metrics for improved understanding and representation of urban dynamics to come up with alternative conceptions of urban spatial structure and change. Particularly with regards to urban forestry, a few studies have focused towards the assessment, mapping, and monitoring of urban forest parameters in fast-growing cities of developing countries. For example, Gong et al. [29] carried out a 30-year forest fragmentation study over Shenzhen Special Economic Zone (SEZ), a city which was established in 1979 in Southern China. Huang et al. [30] utilized satellite images of 77 metropolitan areas in Asia, US, Europe, Latin America, and Australia to calculate and analyze seven spatial metrics (area weighted mean shape, area weighted mean patch fractal dimension index, centrality, compactness index, compactness index of the largest patch, ratio of open space, and density). According to their analysis, the compactness, density, and regularity of urban areas in developing regions generally exceeded the levels reported throughout developed countries [30]. Dewan et al. [31] studied the dynamics of land use/cover changes though landscape fragmentation analysis in Dhaka Metropolitan, Bangladesh, computing and analyzing the following metrics: Number of patches, Patch density, Landscape shape index, Largest patch index, Mean patch size, Area-weighted mean fractal dimension, Interspersion and juxtaposition, Contagion, and Shannon’s diversity index.

In Pakistan, like other developing countries, most urban development is haphazard, typically lacking appropriate planning strategies [10]. Pakistan’s urbanization rate is the highest in South Asia, and by 2030, Pakistan will have more people in cities than in rural areas. Growing population and rapid development is causing prime agricultural land to be encroached and also causing loss of tree cover [32,33,34,35,36]. In the late 1960s, the capital of the Islamic Republic of Pakistan was shifted from Karachi to Islamabad (officially named Islamabad Capital Territory (ICT)). The masterplan of ICT was developed by the famous Greek architect and town planner C. A. Doxiadis [37]. In terms of a planned new capital, ICT is similar to planned new post-colonial capitals/relocations as in the cases of Brasilia (Brazil), Nur-Sultan (named as Astana from 1998 to 2019, and Akmola previously) in Kazakhstan, and Canberra (Australia) [38,39,40]. In the recent decades, with various ongoing development activities, ICT has been struggling with rapid urbanization and gigantic levels of pollution from industrial, residential, and transportation sources. In terms of population, ICT is considered as the most diverse city of Pakistan with a large percentage of immigrants and foreigner population [41]. Unprecedented influx of migrants and population increase has resulted in urban sprawl and conversion of fertile agricultural land and green cover into concrete—a clear deviation from the original ICT master plan [39,40]. Uncontrolled population growth in ICT due to rapid urbanization has deteriorated the living environment, and increased the adverse ecological impacts on human health, flora, and fauna [42].

1.1. Literature Review—ICT Mapping and Monitoring

In the last 20 years, several studies have been conducted on the ICT, regarding land cover change, biomass estimation, water quality monitoring, and temperature increase using satellite datasets. In this section, we present a synthesis of published work regarding land cover change dynamics in the ICT.

Adeel [43] identified urban growth potential through land use for ICT zone IV (Figure 1), based on SPOT-5 2.5 m panchromatic dataset and population census data, and found that nearly 63% of zone IV carries a ’High’ to ’Very High’ future growth potential, which is mainly located close to Islamabad Expressway. This work uses satellite imagery and field data from one year (2007) and does not report spatiotemporal change dynamics. Butt et al. [35] studied the metropolitan development in ICT, based on growth direction and expansion trends from the city center, for the period 1972–2009 using Landsat satellite images. Using Principal Components Analysis (PCA), band ratios, and supervised classification methods, they found that the urban development had expanded by 87.31 km² in 38 years. Butt et al. [44] conducted a study on land cover change analysis over Simly dam watershed, ICT; the results derived from maximum likelihood supervised classification showed tree cover loss of up to 26% and 6% increase in settlements from 1992–2012, based on Landsat 5 TM and SPOT-5 imagery, respectively. Similarly, the watershed analysis of Rawal dam, ICT using Landsat 5 TM imagery, showed 3% degradation of tree cover and 2% gain of settlement from 1992–2012 [45]. Another ICT land cover change dynamics study conducted by Hassan et al. [46] utilized 30 m Landsat 5 TM data for 1992 and 2.5 m SPOT-5 data for 2012, using the maximum likelihood algorithm for image classification. The study revealed a decrease in forest cover of approximately 49% and over 213% gain of settlement area from 1992–2012. Sohail et al. [47] conducted a study to assess the water quality index and analyze the major change in land cover types, vegetation cover, rate of urbanization and its possible impact on groundwater resources, vegetation, and barren land. They used Landsat images for the years 1993, 1997, 2002, 2007, 2013, and 2017 for the assessment and mapping of land cover dynamics; according to their findings, from 1993 to 2017, vegetation areas decreased by 101.77 km², surface water was reduced by 1.10 km², barren land was reduced by 2.90 km², while built-up lands expanded by 105.77 km².

A comparison of Beijing, China and ICT for the role of vegetation in “controlling the eco environmental conditions for sustainable urban environment” was performed by Naeem et al. [42], where they used Gaofen-1 (GF-1) and Landsat-8 Operational Land Imager (OLI) satellite imagery with 8 m and 30 m spatial resolution, respectively. They evaluated various scenarios and models for future development to predict future spatial patterns in both cities. Another study was conducted by Naeem et al. [48] to study the association between green space characteristics, analyzed through landscape metrics, and land surface temperature for sustainable urban environments comparing Beijing, China and ICT.

Khalid et al. [49] conducted a study to quantify the decline of forest reserves and associated temperature variations in a relatively unexplored biodiversity hotspot of ICT, the Margalla Hills National Park (MHNP). In this work, Landsat satellite imagery from 1992, 2000, and 2011 was used to monitor the changes in forest cover and statistical significance tests were used to determine the significance of temperature variation associated with a shift in land cover classes. The study finds that deforestation and forest degradation by local communities is an ongoing practice in MHNP; this necessitates the promotion of conservation practices to minimize ecological disturbances here [49]. Batool and Javaid [50] carried out a study on the assessment of Margalla Hills forest by using Landsat imagery for years 2000 and 2018, and report that the forest cover has decreased from 87% in 2000 to 74% in 2018, whereas built-up area has increased from 5% in 2000 to 7% in 2018, and open land in the study area increased from 2% in 2000 to 7% in 2018. Mannan et al. [51] conducted a study using Landsat imagery, Markov Chain, and Cellular Automata on Margalla Hills, focusing on the quantitative assessment of spatiotemporal land use and land cover changes during 1998, 2008, 2018, and a simulation of 2028. In addition, a forest inventory survey was conducted for biomass and carbon sink estimations. This work shows that the forest area has reduced from 409.36 km² to 392.31 km² and settlement area has increased from 14.97 km² to 39.66 km² from 1998 to 2018. The average yearly biomass and carbon losses were 50.34 Gg/ha/yr and 31.33 Gg C/ha/yr, respectively.

The ICT is a relatively new and spatially heterogeneous city surrounded by the Himalaya mountainous dense forest as compared to other fast growing and expanding cities like Dhaka, Bangladesh [31,52,53,54,55], New Delhi, India [56], Beijing, China [42,48,57,58], Shanghai, China [59,60,61], Tokyo, Japan [27,62], etc. Based on the literature review, we observed that most studies of the ICT land cover dynamics have used different remotely sensed data, methods, definitions, and classification schemes, and have provided diverse results. Most studies which have analyzed the land change dynamics in the ICT focus on the overall analysis of land-cover and land-use change, and a detailed analysis of the landscape ecology and urban forestry characteristics is missing. There has further been very little focus in these studies towards urban landscape metrics and indicators of sustainable urban growth such as LCRPGR.

1.2. Study Objectives

In this paper, well established, proven, and articulated research methodology, satellite datasets, and definitions of features were adopted with the goal to systematically achieve the following defined objectives:

• Detection, measurement, and characterization of land cover features using Landsat medium resolution freely available satellite data (1976, 1990, 2000, 2010, and 2016) and determine the annual rate of change in land cover classes at 10 years interval.
• Landscape metrics and forest fragmentation spatiotemporal analysis, to estimate and report changes in the ICT urban ecosystem over forty years.
• Calculation of SDG indicator number 11.3.1 “Land Consumption Rate to the Population Growth Rate (LCRPGR).”

2. Study Area

The ICT is the capital city of Islamic Republic of Pakistan (Figure 1) located in the Potohar plateau. It comprises an area of 906 km² including mountains and uneven plains exceeding 1,175 m in height above the mean sea level [63]. According to the 2017 national population census, the total population of ICT is approximately two million, which makes it the ninth largest city of Pakistan. It has a humid and sub-tropical climate with four distinct seasons: autumn, spring, summer, and winter. The temperatures vary from 13 °C in January to 38 °C in June. ICT consists of five planning zones: Zone I, II, and V are reserved for planned urban development, while the remaining two zones, III and IV, are managed as National Park and the rural fringes. Marble and chemical factories, steel mills, flour mills, oil units, pigments, paints, and pharmaceutical manufacturing plants are some of the main industries in the city [34].

3. Materials and Methods

3.1. Datasets

To study the land cover dynamics and spatial analysis of ICT, various datasets have been collected from primary and secondary sources. The data collected from the primary sources include Landsat medium resolution satellite imagery and field observations using Geographical Positioning System (GPS) receiver, while the secondary or ancillary data consist of population census and ground truth data from Very High Resolution Satellite (VHRS) imagery.

3.1.1. Landsat Satellite and Digital Elevation Model (DEM) Data

For land cover mapping, we used 30 m spatial resolution orthorectified and cloud free Landsat Multispectral Scanner System (MSS), Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+), and Operational Land Image (OLI) sensor images. The TM, ETM+, and OLI sensors have the same 30 m spatial resolution, while the MSS has a spatial resolution of 57 m. The overall ICT region is covered within one Landsat scene (185 × 185 km). A total of five images for years 1976, 1990, 2000, 2010, and 2016 were downloaded from the United States Geological Survey (USGS) Earth Resources Observation and Science (EROS) archive. All satellite data was collected for the months of October and November, to avoid cloud cover (

Table 1. Landsat images used for land cover mapping.

Acquisition Date	Satellite	Sensor	Path/Row	Spatial Resolution
19 October 1976	Landsat 3	MSS	161/37	57 m
01 October 1990	Landsat 5	TM	150/37	30 m
05 November 2000	Landsat 7	ETM+	150/37	30 m
09 November 2010	Landsat 5	TM	150/37	30 m
24 October 2016	Landsat 8	OLI	150/37	30 m

).

A 30 m spatial resolution Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM) from USGS was used to understand the topography of the area and to prepare the study area maps.

3.1.2. Population Census Data

For the calculation of SDG Goal 11 indicator 11.3.1 LCRPGR, we used the data collected from secondary sources including demographic data (primary census abstracts for the year 1972, 1998, and 2017) from the Pakistan Bureau of Statistics. We used actual census data of 1972 and 1998 to determine the shape of the population curve, which was exponential, as expected. Using the fitted exponential regression function, we derived the growth rate through the exponential function derivative—Equation (1), and then Equation (2) was utilized to determine the approximate population for the year 1976. The population data used for the calculation of LCRPGR is given in

Table 2. Population data of ICT for 1972, 1976, 1998/2000, and 2016/2017.

Year	1972	1976	1998/2000	2016/2017
Population	238,000	291,888	805,235	2,006,572

.

$$\mathrm{Growth}\mathrm{rate}=\left(\mathrm{ab}\right){\mathrm{e}}^{\mathrm{ac}}$$

(1)

$${\mathrm{Pop}}_{1976}={\mathrm{Pop}}_{1972}+4\times \mathrm{Growth}\mathrm{rate}$$

(2)

where a and b are regression coefficients, c represents the population variable, POP₁₉₇₆ = Estimated Population in 1976, and POP₁₉₇₂ = Past Population in 1972.

3.1.3. Ground Truth Data

In this study, we collected ground truth data from two sources: the primary ground data was collected using a Global Positioning System (GPS) receiver in October 2016 during field survey at various locations of ICT, while for the secondary ground data we relied on sub-meter VHRS imagery from Google Earth. A total of 205 samples were gathered, which were utilized as reference data for the satellite image classification and accuracy assessment of results. Using the Google Earth temporal VHRS data, we visually identified and geo-tagged the sites witnessing land degradation, urban development, and tree loss. The geo-tagged dataset helped to validate Landsat-based land cover changes.

3.2. Methodology

Our methodology (Figure 2) on land cover assessment and spatial analysis of ICT consists of (i) pre-processing of Landsat images, (ii) classification of Landsat images, (iii) accuracy assessment of land cover products, (iv) temporal land cover assessment, (v) urban landscape metrics analysis, (vi) urban forest fragmentation analysis, and (vii) LCRPGR calculation.

3.2.1. Pre-Preprocessing Landsat Images

Pre-processing steps for Landsat images included Top of Atmosphere (ToA) reflectance conversion, terrain illumination correction, and layer stacking of spectral bands. A reduction in between-scenes variability was accomplished through normalization for solar irradiance with a two-step process. First, we converted all pixels’ Digital Numbers (DNs) to radiance values using the bias and gain values, which were scene-specific and given in the metadata file of the respective scene. Second, we converted radiance data to ToA reflectance [64]. Terrain illumination correction was conducted using an empirical rotation model proposed by Tan et al. [65]. Layer stacking was performed to construct multi-spectral images at 30 m spatial resolution. The false color composite (FCC) of bands 542 from the multi-spectral images was used to visually interpret land cover features.

3.2.2. Classification of Landsat Satellite Images

The image classification procedure is to instinctively categorize all pixels in an image into land cover classes [66]. In this study, the maximum likelihood classification method was used for land cover-mapping from Landsat images. The maximum likelihood classification method is a supervised classification technique, which works on the basis of multivariate normal probability density function of categories, and utilizes both the variance and co-variance of the spectral response of each unknown pixel to assign it to a particular category [67,68]. In this study, five land cover classes were specified: Tree cover >40% canopy, Tree cover <40% canopy, Settlement, Soil, and Water. Almost 35 training samples were collected for each land cover class to classify the various Landsat images. The classified objects with an area smaller than the Minimum Mapping Unit (MMU) (i.e., 1 ha ~ 3 × 3 pixels) were fused with the neighboring land cover classes [69]. We adopted on-scene digitization technique for land cover change detection. First, we overlaid 2016 base land cover map on the multi-spectral Landsat image of 2010. We traced patches where land changes had occurred, leaving unchanged patches unmodified for consistency [70]. We followed the similar approach to detect land cover change between 2000 and 2010, using 2010 land cover map as a base layer. This process was repeated to generate land cover change maps between 1976–1990, 1990–2000, 2000–2010, and 2010–2016.

3.2.3. Accuracy Assessment of Land Cover Products

Assessment of classification accuracy of 1976, 1990, 2000, 2010, and 2016 land cover maps was carried out to determine the quality of information derived from the data. Random sampling method was adopted to assess the accuracy of satellite derived land cover products. Accuracy was assessed using 125 reference points, based on ground truth data and satellite visual interpretation. The comparison of classification results and reference data was carried out statistically using error matrices. In addition, the non-parametric kappa statistic was computed for each classified map to measure the accuracy of results, as it not only accounts for diagonal elements but also for all elements in the confusion matrix [71].

3.2.4. Temporal Land Cover Assessment

In this study, within the 906 km² total area of ICT, the land cover maps (1976, 1990, 2000, 2010, and 2016) were compared in terms of the area. According to Puyravaud [72], the annual rate of change is based on the compound interest law, considering non-linear change across the timeline to estimate the percentage change per year. In this study, the annual rate of change was calculated using Equation (3) proposed by Puyravaud [72].

$$r=\left(\frac{1}{t_2-t_1}\right)\times (\ln \frac{A_2}{A_1})\times 100$$

(3)

where r is the annual rate of change in percentage, and $A_1$ and $A_2$ is the area at earlier time $t_1$ and later time $t_2$, respectively, and ln denotes the natural logarithm function.

3.2.5. Urban Landscape Metrics Analysis

Landscape metrics are used in this study to quantify the spatial patterns of land cover categories. The landscape metrics can be defined as quantitative and comprehensive measurements showing spatial diversity at a specific scale and resolution [73,74]. In this study, landscape metrics analysis was performed at land cover class level to quantitatively analyze spatial structures and patterns of the topographically and biophysically heterogeneous and diverse landscape of the ICT [74].

Under the three habitat categories Patch Density and Size, Shape and Edge, and Proximity/Isolation, a total of nine landscape metrics were used to quantify structural changes: Number of patches (NP), Mean Patch Size (MPS), Largest Patch Index (LPI), Mean Radius of Gyration (MRG), Mean Shape Index (MSI), Edge Density (ED), Mean Perimeter to Area Ratio (MPAR), Mean Euclidean Nearest Neighbor Distance (MED), and Mean Proximity Index (MPI) (

Table 3. Land cover class level description of landscape metrics used in this study.

Category	Metrics	Unit	Description
Patch size and density	Number of Patches (NP)		Gives the total number of patches in the landscape of a particular category.
	Largest Patch Index (LPI)	Percentage (%)	Ratio of area of largest patch to total landscape area.
	Mean Patch Size (MPS)	Square Kilometer (km2)	Sum of area of all the patches divided by the number of patches of the class.
	Mean Radius of Gyration (MRG)	Meter (m)	Mean distance for each cell in one patch to the patch centroid. It measures connectivity inside habitat patches.
Shape and edge	Edge Density (ED)	Meter per Square Kilometer (m/km2)	Measures total length of edge per unit area. It explains complexity of patch shape.
	Mean Shape Index (MSI)		Measures complexity of patch shape compared to a standard shape of the same size. Its value increases with complexity of shape.
	Mean Perimeter to Area Ratio (MPAR)		Patch shape complexity, based on perimeter length to patch area. It explains shape complexity without standardization to a standard Euclidean shape (square).
Proximity/Isolation	Mean Proximity Index (MPI)		Measure of connectedness of a habitat class. Considers size and proximity of all patches with the same habitat type inside a specified search radius.
	Mean Euclidean Nearest Neighbor Distance (MED)	Meter (m)	Measures minimum edge-to-edge distance to the nearest neighboring patch of the same type. It explains connectedness or isolation in landscape or habitat class.

). The landscape metrics were calculated using the FRAGSTATS v4 software tool. The output of the landscape metric analysis depends upon the spatial resolution of the data [75]; in this particular study, 30 m spatial resolution data has been chosen from 1990 to 2015 and 60 m spatial resolution for 1976.

3.2.6. Urban Forest Fragmentation Analysis

Forest fragmentation is the splitting up of large contiguous forest fields into smaller or less contiguous areas. A number of events or activities can lead to forest fragmentation including road formations, woodcutting, forest conversion to agriculture, forest fires, and human conflict over forest patches [25]. To assess forest fragmentation in ICT, the land cover map is divided into two major categories: forest and non-forest. Forest class consists of tree cover greater and less than 40% canopy cover while non-forest class comprises of settlement, soil, and water land cover classes. The outcome of forest fragmentation analysis was represented into six categories: Patch, Edge, Perforated, Core (<250 acres), Core (250–500 acres), and Core (>500 acres). These categories are signs of forest ecosystem quality and can be used to estimate the amount of fragmentation present in a landscape and the potential habitat impacts [25,76].

3.2.7. Ratio of Land Consumption Rate to the Population Growth Rate (LCRPGR) Calculation

One of the goals of this study is the calculation of the SDGs indicator 11.3.1 Land Consumption Rate to the Population Growth Rate (LCRPGR), which aims at monitoring and measuring urban development by comparing the urban expansion rate with the population growth rate on similar temporal and spatial scales [12,77]. If the LCRPGR ratio value lies between 0 ≤ LCRPGR ≤ 1, it shows the simultaneous increase of population growth rate (PGR) and land consumption rate (LCR), but the land consumption rate is much slower than the population growth rate. On the other hand, if LCRPGR > 1, it reflects the simultaneous increase of PGR and LCR, with a faster LCR than the PGR. To estimate the LCRPGR, satellite data derived land cover maps were used for the years 1976, 1998/2000 and 2016/2017, and census data was used for the specific years (Table 2). The indicator 11.3.1 was assessed at the local level for ICT using the mathematical expressions currently proposed by UN-Habitat, given below in Equations (4)–(6) [77]:

$$\mathrm{LCRPGR}=\frac{\mathrm{Land}\mathrm{Consumption}\mathrm{Rate}}{\mathrm{Population}\mathrm{Growth}\mathrm{Rate}}$$

(4)

$$\mathrm{Land}\mathrm{Consumption}\mathrm{Rate}=\frac{\ln (\frac{{\mathrm{Urb}}_{\mathrm{t}+\mathrm{n}}}{{\mathrm{Urb}}_{\mathrm{t}}})}{\mathrm{n}}$$

(5)

$$\mathrm{Population}\mathrm{Growth}\mathrm{Rate}=\frac{\ln (\frac{{\mathrm{Pop}}_{\mathrm{t}+\mathrm{n}}}{{\mathrm{Pop}}_{\mathrm{t}}})}{\mathrm{n}}$$

(6)

where, ln = Natural logarithm, $\mathrm{Ur}{\mathrm{b}}_{\mathrm{t}+\mathrm{n}}$ = Surface occupied by urban areas at the final year (t + n). $\mathrm{Ur}{\mathrm{b}}_{\mathrm{t}}$ = Surface occupied by urban area at the initial year (t), $\mathrm{Po}{\mathrm{p}}_{\mathrm{t}+\mathrm{n}}$ = Population living in urban areas at the final year (t + n), $\mathrm{Po}{\mathrm{p}}_{\mathrm{t}}$ = Population living in urban areas at the initial year (t), and n = Number of years between the two time intervals.

4. Results

4.1. Accuracy Assessment of Land Cover Products

According to the confusion matrix analysis, 90% overall accuracy and kappa coefficient value of 0.85 was attained for the 2016 classified map. Similarly, overall classification accuracy levels of 88% with kappa coefficient of 0.84 for 2010, 86% with kappa coefficient of 0.82 for 2000, 85% with kappa coefficient of 0.81 for 1990, and 83% with kappa coefficient of 0.79 for 1976 image classifications were achieved (

Table A1. Accuracy assessment of land covers classification for 1976, 1990, 2000, 2010, and 2016.

2016’s Land cover accuracy assessment report
	Tree cover > 40% canopy	Tree cover < 40% canopy	Settlement	Soil	Water	Total	User’s accuracy (%)
Tree cover >40% canopy	22	3	0	0	0	25	88
Tree cover <40% canopy	3	21	0	1	0	25	84
Settlement	0	0	23	3	0	26	88
Soil	0	1	2	21	0	24	88
Water	0	0	0	0	25	25	100
Total	25	25	25	25	25	125
Producer’s accuracy (%)	88	84	92	84	100
Overall accuracy: 0.90, Kappa: 0.85
2010’s Land cover accuracy assessment report
Tree cover >40% canopy	22	5	0	0	1	28	79
Tree cover <40% canopy	3	20	0	0	0	23	87
Settlement	0	0	23	4	0	27	85
Soil	0	0	2	21	0	23	91
Water	0	0	0	0	24	24	100
Total	25	25	25	25	25
Producer’s accuracy (%)	88	80	92	84	96
Overall accuracy: 0.88, Kappa: 0.84
2000’s Land cover accuracy assessment report
Tree cover >40% canopy	22	3	0	0	2	27	81
Tree cover <40% canopy	3	20	1	1	0	25	80
Settlement	0	0	22	3	0	25	88
Soil	0	2	2	21	0	25	84
Water	0	0	0	0	23	23	100
Total	25	25	25	25	25	125
Producer’s accuracy (%)	88	80	88	84	92
Overall accuracy: 0.86, Kappa: 0.82
1990’s Land cover accuracy assessment report
Tree cover >40% canopy	21	2	0	0	1	24	88
Tree cover <40% canopy	4	20	1	1	0	26	77
Settlement	0	0	21	4	0	25	84
Soil	0	3	3	20	0	26	77
Water	0	0	0	0	24	24	100
Total	25	25	25	25	25	125
Producer’s accuracy (%)	84	80	84	80	96
Overall accuracy: 0.85, Kappa: 0.81
1976’s Land cover accuracy assessment report
Tree cover >40% canopy	21	3	0	0	1	25	84
Tree cover <40% canopy	4	21	1	2	0	28	75
Settlement	0	0	20	5	0	25	80
Soil	0	1	4	18	0	23	78
Water	0	0	0	0	24	24	100
Total	25	25	25	25	25	125
Producer’s accuracy (%)	84	84	80	72	96
Overall accuracy: 0.83, Kappa: 0.79

in Appendix A).

4.2. Temporal Land Cover Assessment

The overall settlement in the ICT increased while tree cover classes (greater and less than 40% tree canopies) decreased in forty years (1976 to 2016). The estimated land cover area for each image is summarized in

Table 4. Land cover 1976, 1990, 2000, 2010, and 2016 assessment based on satellite images.

Land Cover Classes	1976	1990	2000	2010	2016
	km2 (%)
Tree cover >40% canopy	182.87 (20.18)	192.19 (21.21)	178.56 (19.70)	140.76 (15.53)	132.27 (10.86)
Tree cover <40% canopy	417.03 (46.02)	342.94 (37.85)	399.31 (44.07)	360.79 (39.82)	306.53 (35.72)
Settlement	29.99 (3.31)	44.51 (4.91)	59.73 (6.60)	143.82 (15.87)	170.40 (18.80)
Soil	270.27 (29.83)	317.91 (35.08)	257.00 (28.37)	254.34 (28.02)	289.45 (33.82)
Water	5.86 (0.65)	8.44 (0.94)	10.51 (1.16)	6.29 (0.69)	6.51 (0.76)
Total	906 (100)

and spatial distributions are presented in Figure 3.

The settlement class in ICT increased from 29.99 km² (3.31%) in 1976 to 170.40 km² (18.80%) in 2016. The tree cover >40% canopy class showed a slight increase, 182.87 km² (20.18%) in 1976 to 192.19 km² (21.21%) in 1990 and then declined to 132.27 km² (14.59%) in 2016. Tree cover <40% canopy land cover class faced overall decline from 1976 to 2016: 417.03 km² (46.02%) in 1976 to 342.94 km² (37.85%) in 1990, increase from 1990 to 2000 to 399.31 km² (44.07%), while again decreased from 2000 to 2016 to 306.53 km² (35.72%). The other two land cover classes, soil and water, observed fluctuating trends in the forty year period. The soil area of 270.27 km² (29.83%) in 1976 increased to 317.91 km² (35.08%) in 1990 and then reduced to 257 km² (28.37%) in 2000, then further decreased to 254.34 (28.02%) in 2010, and again increased to 289.45 km² (33.82%) in 2016. The change in water area was not significant though it has shown fluctuations during the study period; it increased from 5.86 km² (0.64%) in 1976 to 8.44 km² (0.94%) in 1990, and to 10.51km² (1.16%) in 2000, then decreased to 6.29 km² (0.69%) in 2010, and then increased slightly to 6.51 km² (0.76%) in 2016 (Table 4).

Based on the temporal land cover data in the ICT, urban landscape settlement increased at an annual rate of 4.34% since 1976, with the highest annual rate of 8.79% during 2000–2010. The tree cover >40% canopy decreased at an annual rate of 0.81% between 1976 to 2016 and tree cover >40% canopy declined at an annual rate of 0.77% in forty years (1976–2016). The tree cover >40% canopy witnessed the highest loss rate at 2.38% per annum between 2000–2010, and tree cover <40% canopy experienced tree cover loss at 2.72% during 2010–2016 (

Table 5. Land cover annual rate of change (% change per year) from 1976–1990, 1990–2000, 2000–2010, 2010–2016, and 1976–2016.

Land Cover Classes	1976–1990	1990–2000	2000–2010	2010–2016	1976–2016
Tree cover >40% canopy	0.36	−0.74	−2.38	−1.04	−0.81
Tree cover <40% canopy	−1.40	1.52	−1.01	−2.72	−0.77
Settlements	2.82	2.94	8.79	2.83	4.34
Soil	1.16	−2.09	−0.14	2.16	0.17
Water	2.61	2.19	−5.13	1.59	0.42

).

4.3. Urban Landscape Matrix Analysis

Landscape matrix analysis was performed at class level, and the results for each class are reported in Figure 4 and

Table A2. Land cover class level landscape matrix parameters results for 1976, 1990, 2000, 2010, and 2016.

Land Cover Class	Year	Patch Density and Size				Shape and Edge			Proximity/Isolation
		NP	LPI (%)	MPS (km2)	MRG (m)	ED (m/km2)	MSI	MPAR ×103	MPI	MED (m)
Tree cover >40% canopy	1976	1072	12.50	16.98	71.58	21.08	1.27	0.55	869.65	228.47
	1990	5590	9.95	3.44	34.13	44.78	1.24	1.05	1233.12	109.63
	2000	8195	10.31	2.18	26.89	54.25	1.19	1.14	1247.77	98.59
	2010	6121	8.59	2.30	27.34	41.94	1.20	1.15	1401.38	104.15
	2016	2223	7.10	5.95	40.30	23.54	1.28	1.03	1407.03	99.17
Tree cover <40% canopy	1976	1872	30.62	22.29	75.52	66.35	1.32	0.53	4809.76	145.76
	1990	8856	5.89	3.87	40.00	125.30	1.35	1.02	1837.19	77.20
	2000	10145	16.92	3.94	31.18	146.14	1.26	1.11	7773.81	73.00
	2010	9253	29.83	3.89	31.39	122.42	1.23	1.06	18114.32	79.52
	2016	8728	14.97	3.51	36.00	107.08	1.28	1.02	5386.70	83.08
Settlement	1976	4296	0.02	0.70	35.83	17.73	1.09	0.64	1.06	193.66
	1990	7752	0.22	0.57	22.92	26.69	1.13	1.16	12.70	116.02
	2000	6310	0.30	0.95	28.28	30.58	1.19	1.09	33.80	115.14
	2010	11694	5.16	1.23	27.39	68.28	1.20	1.11	784.86	86.70
	2016	12754	8.07	1.34	25.56	71.35	1.18	1.14	1382.81	87.29
Soil	1976	2180	13.76	12.40	65.86	50.43	1.29	0.54	1763.86	159.37
	1990	9271	15.94	3.42	32.34	102.90	1.25	1.06	6172.97	81.38
	2000	12060	3.15	2.13	33.12	103.57	1.27	1.05	564.47	80.48
	2010	8885	4.02	2.86	33.49	84.81	1.26	1.05	530.05	85.99
	2016	12436	5.72	2.33	30.23	107.68	1.26	1.1	1197.42	76.53
Water	1976	8	0.60	68.27	154.90	0.25	1.15	0.54	162.66	2192.65
	1990	242	0.62	3.49	55.79	1.90	1.39	1.02	29.98	362.05
	2000	449	0.61	2.35	43.54	2.93	1.25	1.07	26.36	322.78
	2010	34	0.59	18.51	69.64	0.47	1.21	0.97	65.20	1230.49
	2016	120	0.58	5.43	36.36	0.71	1.18	1.09	46.78	620.00

in Appendix A. The parameters Number of patches (NP), Largest Patch Index (LPI), Edge Density (ED), and Mean Proximity Index (MPI) for the settlement class overall increased from 1976 to 2016. Largest Patch Index (LPI) increased to 16% in 1990 from 13% in 1976, and then reduced to less than 10% in 2000, 2010, and 2016. For tree cover greater than and less than 40% canopy, Number of patches (NP) and Edge Density (ED) initially increased from 1976 to 2000, and then declined from 2000–2016. For tree cover greater than and less than 40% canopy, soil, and settlement classes, the Mean Patch Size (MPS) and Mean Euclidean Nearest Neighbor Distance (MED) sharply declined from 1976 to 1990 and remained lower. In the Mean Euclidean Nearest Neighbor Distance (MED), tree cover greater than and less than 40% canopy, soil and settlement classes remained less than 250 m throughout from 1976 to 2016. From 1976 to 2010, the Mean Radius of Gyration (MRG) decreased for tree cover greater than and less than 40% canopy, and then slightly increased during the years 2010 to 2016. For tree cover >40%, the canopy Mean Proximity Index (MPI) slightly increased between 1976 to 2016 while for tree cover >40%, the canopy MPI decreased from 1976 to 1990, grew between 1990 to 2010, and then again sharply reduced from 2010 to 2016. In all the land cover classes, the Mean Shape Index (MSI) value remained greater than 1 with slight fluctuations. In the Mean Perimeter to Area Ratio (MPAR), the settlement class increased from 640 in 1976 to 1000 in 1990 with slight increase till 2016.

4.4. Urban Forest Fragmentation Analysis

The forest fragmentation analysis results are presented in Figure 5 and Figure 6, and

Table A3. Forest fragmentation results for 1976, 1990, 2000, 2010, and 2016 and percentage assessment based on satellite images).

Forest Fragmentation Classes	1976	1990	2000	2010	2016
	km2
Patch	14.74 (2.46)	15.67 (2.93)	15.46 (2.68)	17.84 (3.56)	19.92 (4.54)
Edge	105.58 (17.62)	119.48 (22.33)	123.42 (21.36)	94.06 (18.75)	96.83 (22.07)
Perforated	35.21 (5.88)	31.52 (5.89)	41.79 (7.23)	24.39 (4.86)	22.74 (5.18)
Core (<250 acres)	41.4 (6.91)	79.94 (14.94)	68.04 (11.77)	52.05 (10.38)	52.24 (11.91)
Core (250–500 acres)	10.99 (1.83)	17.18 (3.21)	15.74 (2.72)	5.54 (1.10)	5.63 (1.28)
Core (>500 acres)	391.98 (65.41)	271.34 (50.71)	313.42 (54.24)	307.67 (61.34)	241.44 (55.02)
Total	585.16 (100)	535.13 (100)	577.87 (100)	501.55 (100)	438.80 (100)

in Appendix A. The analysis shows that the overall core forests of >500 acres decreased from 391.98 km² (65.41%) in 1976 to 241.44 km² (40.29%) in 2016. In the forty years’ time span, the patch class increased to 20 km² (4.54%) in 2016 from 14.74 km² (2.46%) in 1976. The perforated forest fragmentation class was 35.21 km² (5.88%) in 1976, which increased 41.79 km² (7.23%) in 2000, and then again declined to 22.74 km² (5.18%) in 2016. For the edge class, a minor increase was observed from 1976 to 2000 and then a drop in 2016.

4.5. Ratio of Land Consumption Rate to the Population Growth Rate (LCRPGR)

For the years 1976–2000, the LCRPGR ratio (SDG indicator 11.3.1) was obtained to be 0.62, which highlights the simultaneous increase of Population Growth Rate (PGR) and Land Consumption Rate (LCR), but the LCR is much slower than the PGR. For the years 2000–2016, the LCRPGR ratio was obtained to be 1.36, which reflects the simultaneous increase of PGR and LCR, with a faster LCR than PGR.

5. Discussion

The present study was conducted in order to quantitatively analyze the landscape change dynamics in the ICT from 1976 to 2016 using medium spatial resolution Landsat images. In this study, we observed an increase in settlements over the last forty years. The dramatic land cover change in settlements is exerting severe pressure on other land cover classes, particularly tree cover and soil. Already existing urban areas can be seen expanding through rapid construction of residential blocks in the form of housing societies, industrial blocks and road expansions, leading to horizontal and vertical developments in the city and development of lavish farm houses in the vicinity areas (Figure 7). This rapid increase in urbanization is linked also with migrations. Urban growth may have positive or negative impacts on the environment but unplanned growth of urban areas has negative effects. For the economic development of the country, necessary planning is required to make urbanization helpful, as social, health, and environmental issues often accompany the process of urbanization.

In the time span of forty years, an overall decrease has been observed in the area of tree cover classes (i.e., greater than and less than 40% tree canopy) in the ICT, and most of the tree loss has occurred after the year 2000 with a corresponding increase in built-up areas. Apart from tree cutting, another important reason behind rapid tree loss or degradation is forest fires. In the Margalla Hills, forest fires usually take place during dry hot climate conditions when there is no rain for months and temperature goes up to 45 °C. According to Khalid and Ahmad [78], a total of 320 forest fires were recorded from 2002 to 2012 and approximately 8 km² area got burnt as a result. In the ICT, due to uncontrolled urbanization and lack of awareness, huge tree loss has been observed in the last sixteen years, i.e., 2000–2016. Another factor for forest degradation is uncontrolled grazing of livestock [78]. As such, there are no adequate plans or manageable methods to stop grazing activities. For ecotourism and public awareness drives, a number of jogging and hiking trails have been formed in Margalla Hills National Park, and a large number of visitors has severely affected the Margalla Hills National Park by dumping waste. These illegal and unmonitored activities in the forested area cause threats to the forest ecosystem [50,79].

Massive migrations have occurred over the past few years from rural to urban areas, mostly due to low cultivated land output, landlessness, sub-division of land, poor economy, and better educational and health opportunities in urban areas. The rapid increase in population has contributed towards natural resource depletion and rapid deforestation close to settlements [80].

The LCRPGR parameter is an indicator of urban sustainable development, whether urban expansion is in balance with population growth or not. According to literature review, limited scientific peer reviewed studies have been reported on the monitoring and mapping of LCRPGR. Under each SDG, a number of targets and indicators have been defined, which countries have to quantify, but most of the developing countries do not have comprehensive databases through which they can compute, quantify, and report the SDGs indicators. To the best of the authors’ knowledge, only Nicolau et al. [77] and Wang et al. [81] have computed and reported scientific results of SDG proposed LCRPGR over urban areas. Nicolau et al. [77] based their study over the mainland of Portugal, while Wang et al. [81] carried out their study over mainland China, using earth observation and population census data, and reported the increase of LCRPGR value from 1.69 in 1990–2000 to 1.78 in 2000–2010. The LCRPGR related research findings from these studies show that in most cities, both horizontal and vertical urban expansions are carried out in an unplanned manner which has already effected the equilibrium of land consumption versus population increase to attain effective development goals by 2030 [77,81]. In this study over the ICT, the LCRPGR ratio was 0.62 from 1976 to 2000, which increased to 1.36 from 2000 to 2016. Based on studies conducted on the global scale, in the most of the cases LCR is higher than or equal to PGR due to high demand of luxurious occupancies in the urban areas [26,30].

The Shenzhen SEZ and ICT cities were both developed in late 1970s. In terms of urban forest cover change, Shenzhen SEZ had been restored to ~85% (1973–2005) [29] while in this study, we detected urban forest loss of ~27% (1976–2016) in the ICT. In both cities (SEZ and ICT), urban forest fragmentation results revealed losses in forest patches. In South Asia, the landscape of the ICT resembles strongly with the cities of Kathmandu (the capital of Nepal) and Thimphu (the capital of Bhutan), as they are all cities topographically surrounded by the pine trees. As compared to ICT, Kathmandu and Thimphu cities are older but highly migrant resident populated, vastly encroached, and over grazed. Although many studies have investigated the land cover dynamics of Kathmandu and Thimphu cities using temporal satellite data [82,83,84,85], detailed analysis using parameters such as landscape metrics, forest fragmentation, and LCRPGR has not been performed. In this aspect, this study serves as a methodological framework for application and analysis over other similar cities in developing countries.

Due to rapid tree loss and urbanization, the ICT has observed rapid spells of dust storms and soil erosion [86]. Soil erosion, dry temperatures, and consequent dust storms have a negative impact in the form of air pollution [87] and land degradation. The soil in ICT and the surrounding areas is shallow and has a clay composition [63]. The alluvial lands and terraces in the area tend to have low agricultural productivity and in the southern and western parts of the Potohar plateau, the soil is thin and infertile [88]. Streams and ravines cut the loose plain and cause erosion and steep slopes. This land is generally unsuitable for cultivation. However, large patches of deep, fertile soil are found in the depressions and sheltered parts of the plateau and these support small forests and agriculture [89]. Butt et al. [45] described the most significant land degradation issues in the Rawal watershed as soil erosion and loss of soil nutrients. They further showed that, between 1992 to 2012, the majority of land that was previously vegetation, bare soil, or water bodies was converted to agriculture and settlements, suggesting increased pressure on natural resources in the Rawal watershed [45].

The results of landscape analysis in this study reveal that the ICT urban landscape has become more heterogeneous, disproportional and diverse, and tree patches have declined. Alarmingly, core forests of >500 acres have declined almost 15% in forty years. Although at the individual level, the residents of ICT, civil society, and the local government are trying to recover tree cover loss by planting trees, but these initiatives should be continuous and ongoing on a regular basis to monitor the growth of trees without damaging the existing mature standing trees. The temporal forest fragmentation analysis shows that due to tree cover loss, the three categories of core forest fragmentations (i.e., <250 acres, 250–500 acres, and >500 acres) have decreased in forty years (1976–2016). The loss in forest fragmentation negatively influences the habitat and terrestrial biodiversity of ICT [79,90,91]. Based on landscape metrics analysis over the Dhaka metropolitan, a similar Asian developing city, Dewan et al. [31] revealed that cultivated areas and vegetated lands became highly fragmented with increasing anthropogenic disturbances and urban built up category became aggregated and convoluted.

While our analysis of ICT land cover change dynamics in this study is important and unique, there are a few caveats which are worth mentioning. First of all, for the 1976 land cover map, we have relied on approximately 57 m spatial resolution Landsat 3 Multispectral Scanner System (MSS) sensor data which is relatively coarser spatially as compared to Landsat 5 and 8 (i.e., 30 m), which may affect the spatial heterogeneity and accuracy of developed land cover maps. Second, as Landsat is a pioneer Earth Observation (EO) satellite program initiated in 1972, so acquiring remote sensing imagery of the years before that is not possible, and thus we cannot derive land cover maps before the 1960s, when the Islamic Republic of Pakistan’s capital shifted from Karachi to Islamabad. Third, from Landsat 30 m medium spatial resolution satellite data, we cannot detect and delineate the boundaries of built-up areas as well as can be done from sub-meter VHRS images available from the 2000s onward. However, the VHRS datasets come with a high cost, especially when acquisitions at multiples times have to be acquired, and may not be therefore feasible for study sites in developing countries. Of course, research in this domain is ongoing with regards to utilization of VHRS data for study of urban areas, using methods like object-based image analysis, machine learning, etc. Fourth, there may be some level of uncertainty for the field measurement data, as there is always a potential for human error especially when ground truth is collected over larger areas. Fifth, in this study, we utilized only temporal optical satellite data, which may cause optical signal saturation in closed canopy forests, and atmospheric effects (i.e., cloud coverage, haze, and smog).

6. Conclusions

This study sheds new light on systematic land cover dynamics of ICT over the period of four decades, underlining a major decrease in tree cover along with significant increase in settlements and soil. The findings from the analysis of land cover changes, landscape metrics, forest fragmentation, and LCRPGR, are aligned well and agree overall with previous published work on the regional and global scales. The computation of LCRPGR is vital to understand the rate of land change as compared to the population boom, to understand historical land consumption traditions, and to guide decision and policy makers on planned urban expansion while ensuring the protection of environmental, social, and economic assets. It is indeed a scientifically accepted phenomenon that in most of the cities, the land consumption rate overall is much faster than the population growth rate, due to high demand of luxurious occupancies in the urban areas.

Both natural and anthropogenic activities are responsible for land cover changes in ICT. In terms of landscape and fragmentations, the ICT urban landscape ecology is continuously under threat due to the encroachment of housing societies, and influential business community. The findings from urban landscape matrix and forest fragmentation analysis in this study could help ICT Capital Development Authority (CDA) and Islamabad Wildlife Management Board (IWMB) in making strategic decisions to prevent tree loss, forest degradation, and encroachment in the urban landscape of ICT, and also plan future urban growth keeping with the use of remote sensing imagery and geospatial analysis. The methodology outlined in this study is cost-effective and could be easily replicable to other countries/regions, especially in developing countries, through integrating freely available satellite images with a 5–10 years interval.