Two Decades (2000–2020) Measuring Urban Sprawl Using GIS, RS and Landscape Metrics: a Case Study of Municipality of Prishtina (Kosovo)

Since the appearance on Earth, human has been constantly operating in nature, exploiting its riches, but also adapting it to its own needs. Both developing and developed countries are constantly concerned about the urbanization process. Urbanization, in order to be positive, must be developed correctly. If such a thing does not happen, then this development will negatively affect both the environment and human health. In order to develop adequate strategies and policies for the most sustainable and effective land use management, it is necessary to quantify, monitor, determine the factors that have influenced this change in land use and the spread mapping of the urban environment. In this study, Landsat satellite images were used to determine the spatial-temporal characteristics of the urban sprawl environment in the Municipality of Prishtina for a period of 20 years (2000–2020). To map the land cover for Prishtina from 2000 to 2020, the Supervised maximum likelihood classification was used using the Landsat ETM + and OLI data archives in ArcGIS 10.5 software. Using landscape metrics and detection techniques after the classification of satellite images, enabled and assisted in the evaluation and analysis of trends and patterns of urban sprawl. The determination of the changes and the analysis made revealed that during the period 2000–2020, in Prishtina, there was an increase of the built areas by 16.46 km2 at the expense of the unbuilt areas. That there has been an increase in urban areas was also confirmed by the results of landscape metrics.


INTRODUCTION
One of the elements that is very clearly seen in developed and developing countries is urbanization -a global socio-economic phenomenonwhich has led to radical changes in land use. The arrival and increase of the number of urban inhabitants as a result of involuntary migration (poor living conditions) and those voluntary migrations is known as urbanization [Ji et al., 2006; Chikowore & Willemse, 2017]. As we speak, more than half the world's population lives in urban areas. There are reasons why migrations to urban areas in developing and developed countries take place. In the former, it is due to the increase of the population (high increase), while in the latter it is due to the migration of the population [Unhabitat, 2016].
It is well known that the urbanization process does not pass without any effect. The main challenge affecting urban areas from such a phenomenon is the development of informal settlements [Unhabitat, 2016].
Massive population growth in recent decades has caused urbanization-in most of the world-to undergo drastic expansion. Urbanization, in order to be positive, must be developed correctly. If such a thing does not happen, then this development will negatively affect both the environment and human health [Morefield et al. 2018].
Variability in consumption and production models due to the increase in the number of urban inhabitants leads to the emergence of climate change (increase in carbon emissions) [Unhabitat, 2016].
Two Decades (2000-2020) Measuring Urban Sprawl Using GIS, RS and Landscape Metrics: a Case Study of Municipality of Prishtina (Kosovo) Developing countries are more vulnerable (attacked) by the internal breakdown of the city, increased pollution, traffic jams, pressure/housing needs [Mubea et al., 2011;Lee & Chang, 2011;Kityuttachai et al., 2013].
In addition, another major (irreversible) loss due to unprecedented urbanization is the loss of biodiversity and the strengthening of urban morphology [Sudhira et  Even the cities of the Republic of Kosovo are not exempt from urbanization and the challenges it brings. Numerous challenges are highlighted in Kosovo's capital -Prishtina -which is heavily affected by migration, natural population growth, and economic developments which have led to urban growth.
There are different ways and techniques for tracking and detecting changes in the land cover over time. Previously, researchers went straight to the ground and took data and aerial photographs for mapping LULC changes (LULCC)of course, for small areas. However, over time, the need to cover larger areas arose, and the use of these methods was time-consuming, costly, and at the same time tedious [Roberts et al., 2003;Cingolani et al., 2004]. With the development of Remote Sensing, monitoring and detection of LULCC became very easy through satellite imagery. These images, in addition to covering large geographical areas, also include long time periods (historical land use/cover) -ideals for analysis. Remote Sensing also provides data on spaces that are inaccessible (or very difficult to reach) [Fonji & Taff, 2014]. To gain a more accurate and clear understanding of landscape dynamics, as well as to have better land management, it is necessary to detect and monitor LULCC.
In order to have good land management policies and strategies, it is more than necessary to make an assessment and quantification of urban sprawl [Sudhira et al., 2003a,b,c;Ji et al., 2006].
Measuring trends and patterns of urban sprawl requires the use of appropriate algorithms and detailed analysis of land cover in order to quantify the change of urban land cover and urban sprawl.
The advent and development of GIS and Remote Sensing technology have enabled experts to have historical land perspectives and through various techniques to discover changes in urban land use [Sudhira et al., 2003a,b,c;Sudhira & Ramachandra, 2007;Araya & Cabral, 2010]. Now, the assessment of spatial-temporal changes of land use is done through GIS and Remote Sensing techniques, which have proven to be reliable and cost-effective tools [Maktav et al., 2005;Hegazy & Kaloop, 2015].
In order to fully describe urban processes and models, the use of Remote Sensing alone is a bit incomplete. Therefore, the inclusion of landscape metrics complements it [McGarigal et al., 2002;Mallupattu & Sreenivasula Reddy, 2013].
To make the quantitative and complete measurement and description of the basic spacetime structures and models in each landscape, landscape metrics are used [Herold et al., 2005;Aguilera et al., 2011]. These metrics are used to have a clear and unambiguous understanding of the characteristics of the urban landscape in order to sustainably manage urban environments [Mc-Garigal et al., 2002]. In our study area, there is not even a single study on the use of landscape metrics to measure and evaluate urban dynamics.
In relation to landscape metrics, there are a total of 2 types of these metrics: configuration and composition metrics [Gustafson, 1998].
Indices related to the proportion, variety, richness, and presence of the patch without considering the spatial character of the patch are known as composition metrics [Gustafson,1998;Magidi & Ahmed, 2019]. Commonly used of these metrics are Patch richness density (PRD) and Shannon's Diversity Index (SHEI) [Aithal & Ramachandra, 2013;Gustafson, 1998].
Whereas, spatial arrangement, position, character, and orientation of patches in the landscape are elements that refer to configuration metrics [McGarigal & Marks, 1995;McGarigal et al., 2002]. In the landscape, the total area of demolished and unbuilt areas is measured through the Class Area (CA) index in acres/hectare/m 2 / km 2 [Araya & Cabral, 2010].
The measurement of the different dimensions of the urban land use structure is done by means of area-weighted mean patch fractal dimension Determination and quantitative evaluation of spatial models of urban sprawl in the Municipality of Prishtina using GIS, RS and landscape metrics between 2000-2020 is the main purpose of this study.

STUDY AREA
Prishtina is the capital and largest city in the Republic of Kosovo. Geographically, it is located in the north-eastern part of Kosovo, in the geographic latitude 42º40´00´´ and geographic longitude 21º20´15´´, which shows a very convenient central position in the Balkan Peninsula   Figure 2 shows the work steps in a general way, while the details for each step taken will be presented in the following sections.

Pre-processing
In this study, 2 Landsat satellite images of 2000 and 2020 were downloaded from United States Geological Survey (USGS) official website (https://earthexplorer.usgs.gov). In our study, we have taken into account the data of the sensor ETM+ for the year 2000, and the data of the sensor OLI for the year 2020 ( Table 1).
As can be seen from Table 1, the images used in our study were cloudless. Given that the above images are downloaded for free, it makes them more cost-effective and much more efficient.
This data, after downloading, was georeferenced and cut only for our study area (analysis is easier -less load in our software -ArcGIS 10.5).

Land cover/use classification
Based on expert knowledge, numerous training sites have been digitized into the false-color composite (FCC) of the images we downloaded and, eventually, the attributes were added.
To make the LULC classification for our study area, through the use of ArcGIS 10.5 software, the supervised classification method was applied with the maximum likelihood algorithm. The determination of the locations of the sampled sites for each class was performed by visual interpretation of Landsat images (2000 and 2020) based on the land construction information, Google earth, study area knowledge, observations in the field, and historical information. A total of 5 classes were defined: vegetation cover, built-up area, bare land, agricultural land, and water bodies.
The basis of the Maximum Likelihood classification-supervised method-lies in Bayes' theorem. The assignment of pixels in the class with the highest probability is done through the use of a discriminatory function. In this function, the main inputs are the covariance matrix and the class means vector, while their evaluation can be done by the training pixels of a particular class [Halimi et al., 2017].
After supervising the landscaping images, the next step was to reclassify them. This was done by turning them into 2 main classes, such as nonurban areas and urban areas. The characteristics of these 2 classes are described in Table 2.

Accuracy assessment
To understand the results obtained and to use these results in decision-making, accuracy must be assessed. In this study, to make the accuracy assessment, the following were taken into account: producer's accuracy, user's accuracy, overall accuracy, and Kappa coefficient. To determine if the land cover classification was correct, aerial photographs, topographic maps, Google earth, study area knowledge, field observations, and historical information were used. A non-parametric statistical measure that describes the relationships between estimators between categorical variables, not only for diagonal elements but also for variables in the confusion (error) matrix is the Kappa Coefficient. It is calculated using the formula be- where: is the number of rows in the matrix, X is the number of observations in row i and column i (the diagonal elements), !" and !! are the marginal totals of row r and column i, respectively, N is the total number of observations.
We have used the following formula to determine the overall accuracy [Congalton, 1991;Foody, 2002;Congalton, 2005; Verma et al., 2020]: The probability that each pixel in a certain category is correctly classified is defined as the producer's accuracy (PA). We have calculated it using the formula [Verma et al., 2020]: Artificial lakes, rivers, scrubs, forest areas, conifers, grass, and other plantation of different varieties, fallow lands, crop fields, pasture, bare ground, areas with no vegetation cover, uncultivated agricultural land, land with exposed soil, landfill sites, very sparse vegetation.

Built-up area
Areas that include residential, industrial, and commercial areas, mixed-use buildings, roads, and other transport facilities, communication and utilities and other manmade structures, active excavation lands.

Figure 2. Flowchart depicting methodology
The probability that a pixel classified in an image represents that particular category on earth is defined as the user's accuracy (UA). It is calculated as given below: where: !! is the number of correctly classified pixels, N is the total number of pixels, r is the number of rows, and !"#$ and !"#$ are the column and row total, respectively. Using ArcGIS 10.5 software, randomly generated ground checkpoints were generated and the data obtained were compared with reference data to assess the accuracy of each classified map.

Using landscape metrics to quantify urban sprawl
The landscape metrics of each map (2000 and 2020) were calculated. A total of 12 such metrics were used to make such a calculation using the Patch Analyst extension [ESRI, 2015]. This add-on Mean shape index (MSI), Area Weighted Mean Patch Fractal Dimension (AWMPFD), and Mean Perimeter-Area Ratio (MPAR) were used to measure the irregularity or complexity of the shape of urban areas. The degree of complexity of urban patch forms was determined through the AWMPFD [Pili et al., 2017]. The complexity of the shape was measured through MSI and has to do with the complexity of an urban plot in geometric terms [Vaz et al., 2017]. The measurement of the complexity of the urban area was done using MPAR [Yu & Ng, 2007;Pili et al., 2017]. The MPI is calculated to assess the isolation of urban patches. The degree of fragmentation of the pieces and their isolation was determined through MPI.

RESULTS AND DISCUSSIONS
Based on the supervised classification and visual interpretation, the differences in urban areas between the years under consideration are clearly seen. This is shown in the land cover maps of 2000 and 2020 (Figures 3 and 4 respectively).
The results of the LULC classification in which the accuracy assessment was made showed an overall accuracy of 0.90 for 2000 and 0.97 for 2020; producer's accuracy of built-up area for the year 2000 was 94.73 and for the year 2020 was 100; user's accuracy of built-up area for the year 2000 was 90 and for the year 2020 was 95; the Kappa coefficient for the year 2000 was 0.88 and for the year 2020 0.96.
In land cover maps, an increase (visually) in impermeable surfaces is clearly seen (Figures 3  and 4). The land cover maps visually show the spread of urban spaces in the Municipality of Prishtina and that the built-up areas are constantly expanding in environments that have not been built before.
From the results of Table 3, there is a significant change in the growth of the urban environment. The urban environment increased from 13 km 2 (2000) or 2.49% of the total area of the municipality to a total of 29.46 km 2 (2020) or 5.64% of the entire territory of the municipality.
The increase of urbanization, population, and migration towards the city of Prishtina, in our study, has influenced significant changes in LULC classes. Thus, the large population has moved to this area to live and earn more money.In terms of population concentration, on the one hand, and socio-urban chaos, on the other hand, after putting Kosovo under the international military and civilian protectorate, among the urban centers of Kosovo, without a doubt, will lead Prishtina as the   capital of the country. Acceleration of urban development of Prishtina was done through the concentration of other functions, such as economic function, especially in its immediate vicinity, health function, trade, communication, the concentration of entities, and other social, cultural, and scientific institutions. The movement of labor within the day from the place of residence to the place of work and vice versa, and the creation of households with two sources of income (agricultural and non-agricultural) in the vicinity of Prishtina enabled daily contacts with the city and the transfer of contemporary elements to the vicinity. Only a few hundred meters from the center of the capital we encounter many elements which belong to rural settlements, such as manure piles, sheet metal fences, animal stables, etc., while trash can be found in abundance in the center of Prishtina.
From the aspect of urban culture, it is very important how long is the tradition of living in the urban community. Uncontrolled immigration after the war in Prishtina had created an anarchic state. The great concentration, which has created violent urbanization, has caused great problems in the social field, where many households had come to the brink of extinction, hurting even more urban culture and the development of the city in general.

Using landscape metrics to detect changes in urban areas
Landscape metrics were generated using Patch Analyst Extension in ArcGIS software. Based on the results obtained, most of the landscape measurements indicate a positive trend, while some of the indices considered indicate a negative trend.
The results showed a continuous increase in the number of patches (NUMP). NUMP, in 2000 was 4183, while in 2020 they increased to 4848 ( Table 4). The index used to indicate aggregation or division in the landscape is NUMP. According to the results obtained, the cause of the urbanization process is that there were areas built in a non-continuous manner. There was an increase of 13.71% in the number of isolated areas between the period 2000-2020.
There are other determinants of the spread of the urban environment. One such is the size of the patches. The more the urban environment grows, the larger the size of the patches. In this study, we used a total of 5 indices to determine the size of urban patches: Total Edge (TE), Mean Patch Edge (MPE), Patch Size Standard Deviation (PSSD), and Patch Size Coefficient of Variance (PSCoV).
The change in patch size distribution detected by the Patch Size Variation Coefficient (PSCoV) increased from 1718.13 in 2000 to a total of 3049.64 in 2020 ( Table 4). The Average Patch Size (MPS) increased between 2000 and 2020 from 0.31 to 0.60 ha (Table 4). Based on this result we understand that the size of patches will increase with the accumulation of some of the urban patches, and all this is due to urbanization. The same thing happened with the PSSD, which increased from 5.33 in 2000 to 18.52 in 2020.
An increase from 770167.82 in 2000 to 1338411.10 in 2020 also had Total Edge (TE). MPE also increased from 184.11 in 2000 to 276.07 in 2020.
There has been a significant decline in Edge Density (ED) from 592.54 in 2000 to 454.56 in 2020. An increase was observed in Total Edge (TE) by 42.45% between 2000 and 2020. Such a thing happened due to the urban growth which has influenced the increase of the total length of the edge of urban patches.
The Area Weighted Mean Shape Index (AWMSI) also increased with a value from 4.07 in 2000, to 10.09 in 2020 ( Table 4). The increase of AWMSI to 6.02 indicates the complexity and irregularity of urban patches and that such an increase is in line with the growth experienced by the urban environment. Precisely the merging or transformation of non-built areas into built areas during the process of increasing the urban environment has caused the value of AWMSI to increase.
The Mean Shape Index (MSI) also increased for the 20-year period from 1.25 in 2000 to 1.27 in 2020. The complexity of urban patch forms is determined by MSI and, in our case, the value of this index is above number 1. Such a value indicates that urban patches have high geometric complexity.
Increasing the perimeter of the patch and determining the complexity of the shape of urban patches is determined by the Area Weighted Mean Patch Fractal Dimension (AWMPFD). In our case, the value of the AWMPFD increased from 1.40 in 2000 to 1.44 in 2020.
The index which showed a decline over a period of 20 years is Mean Perimeter-Area Ratio (MPAR) which was 1457.06 in 2000, while it fell to 1278.09 in 2020.
In order for researchers to have a clear, deep, and correct understanding of the growth patterns of urban environments, the interpretation of urban growth in the Municipality of Prishtina for the period of 20 years (2000-2020) allows such a thing to be done. For our study area, the main reasons that led to the changes of LULC include migration from rural to urban areas, rapid population growth, non-implementation of existing municipal and urban spatial plans, use of environmentally harmful (obsolete) technologies, etc.
Assessments of the development of the urban environment and the facts established in this development highlight all the changes (desirable and undesirable) that have occurred as a result of urbanization. The importance of determining these changes and the factors that have led to them serve as an extremely important and meaningful method for sustainable decision-making and the formulation of adequate policies.
The importance of the use of landscape metrics lies in defining and evaluating the spatialtemporal changes of urban growth in an easy, straightforward and clear way. In our study area, there was a significant increase in urban areas at the expense of non-urban areas for the period of 20 years.

CONCLUSIONS
In this study, to provide the most accurate and correct understanding of the spread of urban areas, geospatial technology (GIS and Remote Sensing) was used, as well as landscape metrics. In this study, the amount and degree of urban sprawl were determined through the use of landscape metrics.
GIS and Remote Sensing technology was used to quantify and monitor the impacts of urbanization, while landscape metrics were used to obtain spatial data and attributes of urban areas in the Municipality of Prishtina. Using Landsat satellite imagery, urban and non-urban areas are defined through supervised classification to determine spatial-temporal patterns and trends of land cover change in the Municipality of Prishtina for a period of 20 years, from 2000 to 2020.
To determine and understand the spatial trends and distribution patterns of urban areas, the shape index, dimensions, and landscape metrics of the constructed density were used. There was an increase of a total of 16.46 km 2 of urban areas from 2000 to 2020 at the expense of non-urban areas.
Clear views of the growth of urban areas were also observed in the land cover maps. From the obtained results, there is an increase in the complexity of the forms of urban patches. This was discovered using landscape metrics such as AWMPFD, AWMSI, and MSI.
That there was an increase in an urban area over the 20-year period was revealed through isolation, relative size, complexity, and absolute size.
Through the results placed in the tables and the use of landscape metrics, it was clearly shown that there was a significant increase in urban areas which were clearly illustrated in the land cover maps for the period of 20 years. For the period of 20 years, there are clear and significant changes of land cover in Prishtina, which through geospatial technology can be determined and monitored.
In order to be able to quantify the uncontrolled violation of the natural environment by the uncontrolled growth of urban areas, there is a need to apply sound urban planning strategies against rapid urban sprawl. The application of GIS and Remote Sensing technologies offers excellent opportunities to provide the right information for decision-makers, with the sole purpose of ensuring sustainable management of urban areas.
This study demonstrated that the analysis of the assessment of the amount of landscape change from urban sprawl in the Municipality of Prishtina can be done through the use of Landsat satellite images accompanied by landscape metrics. Thus, all of these can be used for the right information of decision-makers on uncontrolled man-made landscape changes, especially on illegal constructions which are among the main elements of this irregular growth of urban areasevidenced by the use of AWMPFDI and AWMSI.
Looking at the changes in the number of patches, Total Edge and MPS conclude that the change of spatial configuration is possible, while the further increase of urban areas will cause an increase in landscape diversity in general.