Evaluation of Change in Land Usage and Land Cover in Karaj, Iran

. In this study, classification results were derived from remote sensing data and the Support Vector Machine (SVM) algorithm used in this process, which classifies Landsat land-cover images. The accuracy of image classifications was evaluated by calculation of the Kappa coefficient. The area of study is Karaj, the capital of Alborz province, in north-central Iran. It is situated in the foothills of the Alborz Mountains and occupies a fertile agricultural plain. Landsat data used in the classification of land cover were collected from USGS websites, and multi-temporal images from the data were geometrically corrected. After this process, we calculated 11 metrics at the landscape and class-level scales: five metrics of class level and six metrics of landscape. The results showed that the landscape patterns in Karaj were changed due to the process of urbanization over an 11-year period. At the class level, for all classifications, the AI metric increased and the PD and NP metrics decreased. At the landscape level, the PD, ED, NP, and SHDI metrics decreased, and LPI and AI increased. These results provide insights about urban development policies and about whether the expansion of urban areas is beneficial for environmental sustainability in Iran and elsewhere in the world.


Introduction
Land is not only the locus of terrestrial natural ecosystem functionality but has also been used by humans in numerous ways (Song & Deng, 2017).Land use change has been recognized as the key human-induced effect on ecosystems (Kaczorowska et al., 2016).Land-use and land-cover mapping are currently important in the areas of science, research, planning, and management (Singh & Dubey, 2012).Remote sensing data, processed using geographic information system (GIS) software, can significantly help in the creation of land-use maps (Haas & Ban, 2017).Identifying and monitoring changes in land use and land cover is a complex process (Sun & Zhou, 2016).Therefore the understanding of recent land-use changes and how such changes will occur in the future is of fundamental importance (Rounsevell et al., 2006).This understanding is important for decision-support procedures to recognize appropriate land-use policies and for environmentalists and planners in the development of plans to tackle environmental issues (Romano et al., 2018).In this regard, the most important way to understand land-use and land-cover (LULC) change is to analyse changes in the landscape pattern (Fan & Ding, 2016).Urban expansion causes a change in population and landscape patterns (Zhang & Su, 2016).Investigating the relationship between urban expansion and landscape patterns can provide support for urban environmental management (Huilei et al., 2017).Landscape patterns can be quantified by landscape metrics, which is one of the key tools for landscape assessment and management (Li & Wu, 2004).Studies have shown that high-intensity human activity has led to LULC change (e.g., Tolessa et al., 2017;Sun et al., 2018;Hu et al., 2019).In this study we investigated the impact of LULC changes on landscape patterns in the city of Karaj, the capital of Alborz Province, in north-central Iran, with the combined use of satellite remote sensing, GIS, and FRAGSTATS software.

Location and description of the study area
Karaj is Iran's fourth largest city and is located in Alborz Province, to the west of Tehran.The Karaj city landscape covers an area of 117520 ha and is located between 35 • 46ʹ-36 • 09ʹ N latitude and 50 • 46ʹ-51 • 21ʹ E longitude (Figure 1).The altitude of Karaj is between 900 and 1700 m above sea level.Karaj is undergoing rapid urbanization.Urban expansion has directly affected Karaj's urban structure, so that the rural population has been attracted to urban areas, and the growing urbanization has led to environmental changes.

Processing of the satellite images (classification)
The changes in Karaj's urban land-cover structure have directly produced environmental changes.For this study, we used two years of Landsat satellite images (2006 and 2017) and a classification method (Table 1).The spatial resolution of the Landsat satellite images was 30 m.The data were obtained free from the United States Geological Survey (USGS) platform (US Geological Survey, 2019).Because clouds are an obstacle for the interpretation and classification of satellite images (Silva et al., 2019), the chosen images from Landsat-5 (Thematic Mapper [TM] sensor) and Landsat-8 (Operational Land Imager [OLI] sensor) do not contain clouds over the study area.Supervised and unsupervised classification are the two main classification methods used for spectral images.Supervised classification uses training samples for classification, resulting in greater accuracy (Wang et al., 2018).Land-cover changes were identified by an image classification method using the Support Vector Machine (SVM) algorithm with EVNI 5.1 software.In a majority of geospatial studies, SVM has been applied to the classification of remotely sensed data, and it has been used for modelling LULC changes (e.g., Samardžić-Petrović et al., 2017;Karimi et al., 2019;Heydari & Mountrakis, 2019).SVM is a supervised classification method derived from statistical learning theory that often yields good classification results from complex and noisy (raster) data.We selected 10 land-cover classes from the classification parameters within the software: human made, agriculture, garden, water body, low dense grassland, dense grassland, barren, rocky outcrop, green space, and river.A kernel function mathematically represents the weights of nearby data points in estimating the target classes.The Radial Basis Function (RBF) kernel type works well in most cases.The listed mathematical definitions of each kernel function are as follows (Chang & Lin, 2001;Wu et al., 2004;Hsu et al., 2007;Chih-Wei et al., 2010): (2) where i x represents scattered data at point i; j x represents scattered data at point j; xj T denotes the value of function at point xj ; g is the gamma term in the kernel function for all kernel types except linear; d is the polynomial degree term in the kernel function for the polynomial kernel; r is the bias term in e kernel function for the polynomial and sigmoid kernels; and g, d, and r are user-controlled parameters, as their correct definition significantly increases the accuracy of the SVM solution.
In the study, we used the RBF kernel type.Therefore, gamma in the kernel function = 0.250.This value is greater than zero.We considered using the value of the Penalty Parameter (100) for the SVM algorithm.In addition, the default had value 0 and was set for Pyramid Levels under the classification probability threshold.
The accuracy (quality) of classification comes from the calculation of the kappa coefficient (Silva et al., 2019).This was obtained from the following equation (Cohen, 1960;Silva et al., 2019): where I is the kappa coefficient; ii x is the value in row i and column i, x  is the sum of row i, x  is the sum of column i of the matrix, n is the total number of observations, and c is the total number of classes.

Landscape metrics features analysis
The maps with their classifications were used to evaluate the landscape metrics.The metrics display numerical information about landscape composition, configuration, and dimension (Martinez del Castillo et al., 2015).Composition is a non-spatially-explicit characteristic.It does not measure or reflect the patch geometry or geographic location (Leitao & Ahern, 2003).This includes the metrics SHDI, NP, CA, and the dimension index.dimension (e.g.FRAC metric), which is appealing because it reflects shape complexity across a range of spatial scales.Landscape configuration relates to spatially explicit characteristics of land-cover types in a given landscape, namely those associated with patch geometry or with the spatial distribution of patches (Leitao & Ahern, 2003).It includes the metrics ED and PD.Composition and configuration are two basic components of landscape structure (Lacoste et al., 2014).Landscape metrics can be applied at three different scales: landscape, class, and patch level (McGarigal & Marks, 1995;Martinez del Castillo et al., 2015).We used landscape-level and class-level metrics.The patch-level metrics are not useful for our purposes.
The class-level metrics are computed for every patch type or class in the landscape.The resulting output file contains a row (observation vector) for every class, where the columns (fields) represent the individual metrics.
The landscape level is computed for the entire patch mosaic.The resulting landscape output file contains a single row (observation vector) for the landscape, where the columns (fields) represent the individual metrics.The landscape metrics features were calculated and analysed by FRAGSTATS software 4.2 (McGarigal et al., 2002), and the shape file format of land use data (image land use in years 2006 and 2017) was converted to a raster format and brought into FRAGSTATS (Zhang & Gao, 2016).The landscape structure was represented in the landscape mosaic model (Zhang et al., 2015).
In this study, we calculated six landscape levels (PD, NP, ED, LPI, AI, and SHDI) and five class levels (PD, NP, ED, LPI, and AI).The metric SHDI was calculated only at the landscape level (Table 2).These metrics are defined as follows:  Largest Patch Index (LPI)areas of the largest patch in the landscape divided by total landscape area. Patch Density (PD)number of patches of corresponding patch type (class) per unit area. Edge Density (ED)sum of the lengths (m) of all edge segments in the landscape divided by the total landscape area (m 2 ). Shannon's Diversity Index (SHDI)popular measure of diversity in community ecology, applied here to landscapes.SHDI reflects the landscape heterogeneity and is sensitive to the non-equilibrium distribution of patch types in the landscape; SHDI = 0 when the landscape contains only one patch (i.e., no diversity).SHDI increases as the number of different patch types (i.e., patch richness, PR) increases and/or the proportional distribution of area among patch types becomes more equitable. Number of Patches (NP)simple measure of the extent of subdivision or fragmentation of the patch type.
NP is calculated for one selected class.NP = 1 when the landscape contains only one patch of the corresponding patch type, that is, when the class consists of a single patch.NP ≥ 1, without limit.
 Aggregation index (AI)ratio of the observed number of like adjacencies to the maximum possible number of like adjacencies given the proportion of the landscape consisting of each patch type.These metrics, which can reflect the composition, shape of the patches, and aggregation of the landscape in Karaj, have been used frequently by previous researchers (e.g., Su et al., 2011;Zhang & Gao, 2016;Huilei et al., 2017;Zarandian et al., 2017).The accuracy of the classification is denoted by the kappa coefficient, I = 0.93% ( 2006) and I = 0.95% (2017).According to the USGS standard, the minimum acceptable kappa coefficient is 85% (USGS, 2019).The accuracy obtained from the classification is acceptable.The maps were used for the next analyses.Human-made changes increased from 8 343.18 ha ( 2006) to 14 478.66 ha (2017).The area defined by this classification will increase in the future.The area classified as agriculture decreased from 11 181.42 ha ( 2006) to 6 511.86 ha (2017).These changes reflect rapid urbanization in Karaj, which can lead to numerous changes in ecosystem services (Hu et al., 2019).In fact, human activities can cause irreversible changes in the environment (Li et al., 2016) as natural ecosystems become semi-natural and artificial, which is a major threat to the structure and function of ecosystems (Song & Deng, 2017).Urban expansion causes problems such as destruction of agricultural land, reduction in green space, water pollution, soil erosion, and loss of environmental quality.
Low dense grassland grew from 20 254.The lowest number of patches was observed in the river class (NP = 1) in both years.The decrease in the number of patches for human-made, agriculture, garden, water body, low dense grassland, dense grassland, Barren, rocky outcrop, and green space is due to the increase their size and integration from 2006 to 2017.
PD characterizes the fragmentation level of the landscape (Zang et al., 2017).The PD for all types of classes declined LPI is used as an indicator of landscape fragmentation (Sun & Zhou, 2016;Hassan, 2017).LPI approaches 0 when the largest patch of the corresponding patch type is increasingly small and 100 when the entire landscape consists of a single patch.The LPI in the agriculture class sharply declined but increased in the human-made class.In 2006 the dominant cover of the study landscape was agriculture, and in 2017 it was human-made.Due to rapid urbanization, this metric declined in the agricultural areas.The AI metric increased from 2006 to 2017 for all classes except river.The highest AI metric in 2006 and 2017 was in the human-made class.The reason for the increase in the human-made class was the destruction of agriculture and dense grassland for construction.

Changes in landscape level
The landscape level was calculated for PD, NP, ED, LPI, AI, and SHDI (Figure 5).The SHDI metric was calculated only on the landscape level.From 2006 to 2017 the NP metric at the landscape level significantly decreased, from 3 574 to 1 532.The PD metric decreased, from 3.04 to 1.3.During the study period, construction land expanded with the corresponding loss of cultivated land and grassland, which reduced the PD metric and made the whole landscape more agglomerated and regular.The ED metric at the landscape level also decreased.ED and PD metrics both measure the complexity of the shape of patch and represent the spatial heterogeneity of a landscape mosaic (Wei & Zongyi, 2012).
For the SHDI metric there was a very slight decrease (from 1.72 to 1.69) due to increased human interference.Shannon index metric decrease, Shows heterogeneity decrease in Karaj city land use Generally, the urbanization process caused the PD and SHDI decrease and resulted in the landscape patches becoming more regular and aggregated.The LPI metric increased significantly, from 5% in 2006 to 20% in 2017.The AI metric increased during this period due to construction and continuous filling in of cultivated lands.

Conclusions
Significant LULC changes were observed in Karaj from 2006 to 2017, especially in human-made, agriculture, and garden.Urbanization, population increase, and human activities were the main causes of land use changes.These changes have had a direct impact on the value of ecosystem services, and they also threaten ecosystem structures and ecological processes.Both class-level and landscape-level metrics computed in this study show diverse landscape patterns and land fragmentation processes.The direct cause of changes in landscape patterns in Karaj is the urbanization of land use.Spatial changes in land cover over time also affect the provision of ecosystem services in the landscape and the fragmentation of green areas.The evaluation of changes in land use and land cover is vital to support decision-making processes at different levels.These study results also provide information about the expansion of urban areas and thus are beneficial for environmental sustainability and urban development policies in Iran and elsewhere in the world.

F
. Mohammadyari et al.Evaluation of change in land usage and land cover in Karaj, Iran 2 1. Materials and methods

Figure 1 .
Figure 1.Map of the study area (source: authors) proportional abundance of each patch type in the landscape % PD Number of patches of the corresponding patch type divided by total landscape area (lengths (m) of all edge segments involving the corresponding type, divided by the total landscape area (m 2 ) m/ha SHDI Sum, across all patch types, of the proportional abundance of each patch type multiplied by that proportion No units NP Number of patches of the corresponding patch type (class) No units AI Number of like adjacencies involving the corresponding class, divided by the maximum possible number of like adjacencies involving the corresponding class, which is achieved when the class is maximally clumped into a single, compact patch % 2. Result and discussion 2.1.Land-cover changes in Karaj from 2006 to 2017 The classification of images from Landsat-5 (2006) and Landsat-8 (2017) by SVM classification shows 10 land-use classes in the study area by ENVI software, representing visually significant changes in the Karaj area (Figure 2).

Figure 2 .
Figure 2. Land-cover classification, 2006 and 2017 (source: authors) 23 ha (2006) to 22 586.96ha (2017), and dense grassland decreased from 43 221.87 to 40 167.09ha.These changes indicate that Karaj grasslands are under threat.A reduction in grasslands can increase climate fluctuations.Rocky outcrop areas declined from 21 189 ha (2006) to 24 393.78 ha (2017).This parameter is influenced by human-made changes.Green space increased from 780.57ha (2006) to 802.42 ha (2017).Areas classified as water body underwent a helpful increase from 305.55 ha to 326.97 ha.Those classified as garden fell from 10 674.99 ha to 7 036.67 ha and barren decreased from 578.97 ha to 227.43 ha.All classification results are shown in Figure 3.

Figure 4 .
Figure 4. Spatio-temporal land cover metrics at the class level (PD, NP, LPI, ED, AI) (source: authors) . The highest PD in 2006 was observed in the garden class (PD = 0.78), but by 2017 it had declined (PD = 0.29).Thus all classes, especially the garden class, became fragmented over time.The decrease in the PD and NP metrics indicates the continuity of the urban landscape in Karaj at the class level.The ED metric declined from 2006 to 2017 for all classes except human-made.The ED for the human-made class was 4.19 in 2006 and increased to 4.48 in 2017, indicating an increase in the dispersion and irregularity of the human-made patches.