## 3 Data and Methods
### Remote Sensing Images Processing  
The two land uses/land covers, for 1989 and 2019, were extracted from remote sensing imagery. This method of data acquisition was employed in numerous studies, focusing on highlighting the effects of land-use changes on different phenomena [59; 60; 61; 62]. In the literature, there are many different methods used to extract the land cover from remotely sensed imagery. Thus, the unsupervised classification is one of the methods used to extract the land uses/land covers. This method can achieve very good results when the values belonging to a certain land-use/land-cover class are very close in terms of the spatial distribution [63]. Regarding the supervised classification, one of the most known methods is modified k-nearest neighbor [64], which allows the users to validate the value of each pixel according to its neighbors. Before the development of the PCs processing power, the minimum-distance-to-means classification method was considered faster than other methods because it is based on a simple mathematical algorithm [65]. The advantage of the support vector machine classification method lies in the fact that it allows for the analysis of high dimensional datasets [66]. Maximum likelihood is another very popular classification method that is based on a statistical rule that analyses the probability of each pixel belonging to a particular class [67]. In a study carried out on Istanbul city, [PERSON] et al. [68] showed that the maximum likelihood algorithm can obtain a higher classification accuracy than other methods (learner vector quantization) or a lower classification accuracy (than multilayer perceptron) depending on different factors. Nevertheless, given the strengths of the maximum likelihood algorithm, such as the advantages of both the mean vectors and the multivariate spreads of each class and the identification of those elongated classes [69], we decided to also use the maximum likelihood in the present analysis. The primary source for land-use/land-cover classification was the two aforementioned Landsat scenes (Figure 3a,b) (spatial resolution \(30\times 30\)m), acquired on 18 August 1989 and 17 May 2019. The scenes, provided by USGS Earth Explorer, were pre-processed and then classified using ENVI 5.0 software, made by L3 HARRIS Geospatial, Broomfield, Colorado--United States of America. Only the reflective bands (1-5 and 7) of the sensor were used. In a first stage, in the image pre-processing procedure, the radiometric and relative atmospheric corrections were performed in order to eliminate the effects of multiple factors, such as the changes in sensor characteristics, atmospheric condition, solar angle, and sensor view angle [70; 71; 72; 73]. In terms of radiometric corrections, it was mandatory to convert the digital number of pixels into spectral radiance values, and further, the radiance into reflectance [74; 75; 76]. Furthermore, the dark object subtraction method was used for the relative atmospheric correction [77; 78]. After pre-processing, the supervised maximum likelihood image classification was performed [79]. The first step included the identification of many ground truth regions of interest belonging to seven land-use/land-cover categories detected in Landsat imagery. In this respect, we tried to extract, based on the satellite image, the best training and testing areas to extract the seven land-use/land-cover datasets. We correlated the information found in the image with different land-use/land-cover information available online (we used the official Corine Land Cover dataset [80]), and we assigned the different reflectance values of the satellite imagery to the land-use/land-cover class that they overlaid/corresponded with. In this way, different reflectance values that corresponded with the same land use/land cover were assigned correctly in a more accurate final classification. In accordance with [PERSON] et al. [81], the regions of interest were divided into training areas (80%) and testing areas (20%).