*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.

**Figure 2.** Methodological flow chart.

#### 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<sup>2</sup> 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 \left(\ln \frac{A\_2}{A\_1}\right) \times 100\tag{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). 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.


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

#### 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]:

$$\text{LCRPGR} = \frac{\text{Land Consumption Rate}}{\text{Population Growth Rate}} \tag{4}$$

$$\text{Land Consumption Rate} = \frac{\ln(\frac{\text{Urb}\_{l+n}}{\text{Urb}\_l})}{\text{n}} \tag{5}$$

$$\text{Population Growth Rate} = \frac{\ln(\frac{\text{Pop}\_{t \to n}}{\text{Pop}\_{t}})}{\text{n}} \tag{6}$$

where, ln = Natural logarithm, Urbt<sup>+</sup>n = Surface occupied by urban areas at the final year (t + n). Urbt = Surface occupied by urban area at the initial year (t), Popt<sup>+</sup>n = Population living in urban areas at the final year (t + n), Popt = Population living in urban areas at the initial year (t), and n = Number of years between the two time intervals.
