Quantifying Land Cover Changes in a Mediterranean Environment Using Landsat TM and Support Vector Machines

The rapid advent in geoinformation technologies, such as Earth Observation (EO) and Geographical Information Systems (GIS), has made it possible to observe and monitor the Earth’s environment on variable geographical scales and analyze those changes in both time and space. This study explores the synergistic use of Landsat EO imagery and Support Vector Machines (SVMs) in obtaining Land Use/Land Cover (LULC) mapping and quantifying its spatio-temporal changes for the municipality of Mandra–Idyllia, Attica Region, Greece. The study area is representative of typical Mediterranean landscape in terms of physical structure and coverage of species composition. Landsat TM (Thematic Mapper) images from 1993, 2001 and 2010 were acquired, pre-processed and classified using the SVMs classifier. A total of nine basic classes were established. Eight spectral band ratios were created in order to incorporate them in the initial variables of the image. For validating the classification, in-situ data were collected for each LULC type during several field surveys that were conducted in the area. The overall classification accuracy for 1993, 2001 and 2010 Landsat images was reported as 89.85%, 91.01% and 90.24%, respectively, and with a statistical factor (K) of 0.96, 0.89 and 0.99, respectively. The classification results showed that the total extent of forests within the studied period represents the predominant LULC, despite the intense human presence and its impacts. A marginal change happened in the forest cover from 1993 to 2010, although mixed forest decreased significantly during the studied period. This information is very important for future management of the natural resources in the studied area and for understanding the pressures of the anthropogenic activities on the natural environment. All in all, the present study demonstrated the considerable promise towards the support of geoinformation technologies in sustainable environmental development and prudent resource management. Forests 2020, 11, 750; doi:10.3390/f11070750 www.mdpi.com/journal/forests Forests 2020, 11, 750 2 of 19


Introduction
Land Use and Land Cover (LULC) constitutes a key variable of the Earth's system that has in general shown a close correlation with human activities and the physical environment [1,2]. The fast development of human societies, particularly since the beginning of industrial revolution, has intensified different types of human activities that have resulted in a continuous and noticeable influence on LULC [3]. As natural and semi-natural habitats are continuously exposed to growing pressure due to various anthropogenic activities, conservation and sustainable land use have become a priority [4]. Quantifying the temporal and spatial patterns of LULC change and its corresponding consequences has also been recognized today as a highly significant topic in land change-related sciences [5]. Understanding the driving forces behind the development of land use in the past, managing the current situation and modeling LULC for the future aids the scientific community in developing plans for multiple uses of natural resources and nature conservation. The type of change that is of interest in each case is related to either short-term phenomena, such as floods and snow cover or the occurrence of a fire, or to longer-term phenomena, such as urban sprawl and desertification [6][7][8][9].
Geoinformation technologies, mainly Earth Observation (EO) and Geographic Information Systems (GIS), provide the preferred approach nowadays for analyzing, designing and modeling LULC. Among the key advantages of EO as compared to conventional field surveys is the provision of timely and often inexpensive imagery at an adequate spatial resolution at low or no cost from the local to the regional and global scales. Looking at the relevant literature, one can observe that many studies over time have used satellite data, and in particular Landsat satellite imagery, to identify and record LULC [10,11]. In comparison to other commercial medium spatial resolution radiometers (e.g., ASTER, ALOS, SPOT), or to freely-distributed coarser spatial resolution imagery (e.g., MODIS, MERIS, AVHRR), Landsat has a number of advantages that make it unique in mapping mining activity and reclamation. A series of platforms launched since 1972 (i.e., MSS, TM, ETM+) have allowed the collection at no cost today of spectral observations globally acquired from the visible to the thermal infrared parts of the electromagnetic spectrum at a spatial resolution ranging from 30 to 120 m and at a temporal resolution of 16-18 days.
The promise of the Landsat sensor in providing straightforward and accurate mining and reclamation mapping when combined with different types of image processing techniques was also underlined even in the early years of the sensor launch. Indeed, looking back at the techniques used to detect changes in land use, the pioneers were Angelici et al. [12], who developed the first technique for detecting land use changes using Landsat images. Initial investigations focused on finding changes in vegetation. A typical example is the study by Allum and Dreisinger [13], which used Landsat images from 1973 and 1983 to create maps of vegetation change using a methodology that is useful for highlighting abrupt rather than gradual changes. Li and Yeh [14] applied a combined methodology to highlight the changes of urban areas in the Pearl River Delta region of China. They first formed principal components of multi-temporal data and supervised classifications by using the maximum likelihood method. Zha et al. [15] used the Normalized Difference Built-up Index (NDBI) to automate the mapping process of residential areas. Landsat TM (Thematic Mapper) satellite imagery was used. The results produced a reliable mapping solution with an accuracy of 92.6%.
Generally, producing thematic maps for information extraction related to LULC changes using remote sensing data has been based largely on performing digital image classification [16]. Generally, a variety of classification techniques of different levels of complexity and assumptions linked to their implementation in practice have been utilized to portray LULC changes. An overview of the available classification approaches applied to remote sensing imagery including their relative strengths and limitations was made available recently by Lu and Weng [17].
Yet, it is generally admitted that even until today accurate production of LULC thematic maps and their changes from EO observations remains an important challenge [18]. As a result, a large amount of effort has been directed towards developing more sophisticated algorithmic approaches, aiming to map LULC and its changes more accurately and robustly. This rapid development in image processing algorithms has also necessitated the need of their all-inclusive assessment with different types of EO imagery.
In this context, specifically the use of machine learning-based classifiers, such as Support Vector Machines (SVMs) [19], with satellite observations widely used in LULC-related applications, such as those from the Landsat sensor, comprises a topic of key importance and priority interest today. SVMs are a supervised classifier based on statistical learning theory. SVMs are one of the primary algorithms used for supervised non-parametric classification. SVMs belong to a group of theoretically superior machine learning algorithms [20]. Their use is becoming increasingly popular in remote sensing applications due to the several benefits of SVMs over other well-established classification methods [21]. SVMs have several advantages in comparison to other classifiers, as they can be applied requiring limited effort in their training and they do not make any assumption for the probability distribution of the training datasets used, as for example is the case for the maximum likelihood parametric classifier. Moreover, they are able to deal easily with high dimensionality datasets and have proven to effectively address the ill-posted problems providing high classification accuracies. SVMs have generally been implemented in many classification problems using different class numbers employed in classification and remote sensing data acquired at different spatial scales (e.g., SPOT, MODIS, Landsat TM/ETM+) producing reliable results [22][23][24][25].
Several studies have addressed the use of Landsat TM in LULC mapping and change detection at various geographical regions including the Mediterranean region. Many of these studies have been based on the implementation of image classification approaches, demonstrating the usefulness of different classification approaches [26][27][28][29][30]. For example, DiGirolamo [31] used two Landsat TM images for 1991 and TM+ for 2000 in an area near Atlanta (Gwinnett) in Georgia, and the differences were identified using three methods. Mancino et al. [32], via Landsat TM satellite data, in Basilicata in southern Italy, used vegetation markers to detect changes in land cover. The results of their study showed the current trend of natural expansion of forests occurring in drought-prone areas of the Mediterranean. The information on the land cover dynamics and of forest cover in particular shown in this study can be considered as a useful starting point for further analysis of spatial and temporal models of vegetation changes in degraded areas. Abd and Al-Najjar [33] studied the LULC changes of the wider Johor Bahru area of Malaysia using a maximum likelihood supervised classification with four classes of urban, water, vegetation and bare land for two series from 1995 and 2011, with Landsat TM and + ETM satellite data, respectively. Aclassification accuracy of 84.14% and 89.11%, respectively was reported. Detection of changes showed that urban areas increased by 3% and areas covered by vegetation decreased by 3%, while bare areas remained stable. In 2015, Giatanis et al. [34], studied 60 Years of Land Cover Change in the Marathon Area, Greece (next to the study area of this work) using a variety of data, such as panchromatic orthophotos, forest maps, vegetation maps etc., detecting a significant decrease in forest cover.
SVMs are a machine learning technique adopted in solving binary classification problems [35]. They project the points of the whole set of training in a space of more dimensions and find the superstructure that optimally separates the points of the two classes. The unknown points are classified according to the side of the superstructure in which they are located. The vectors that define the superstructure that separates the two classes are called support vectors. They have low computing costs, even in the case of non-linearity. Their peculiarity lies in the fact that they separate the data with the maximum possible distance, making the fewest possible errors. The SVMs classifier algorithm uses a kernel function, in the form of four different types of functions [36]. Arenas-Castro et al. [37] used SVMs over QuickBird images for mapping Iberian wild pear trees. In 2017, Iglesias et al. [38] applied SVMs in order to develop a tool to predict the influence of heartwood on wood density and pulp properties. Ramezan et al. [39] used this context in order to examine four sample selection methods by using data such as orthoimagery and LIDAR. In the same year, Xie et al. [40] used SVMs, among other classifiers, for the classification of land cover, forest, and tree species classes in China. More recently, SVMs have been used to for estimating the growing stock volume of Pinus massoniana plantations and measure the expansion of Eastern Redcedar into the deciduous woodlands within the Forest ( [41,42], respectively). However, to our knowledge, very few studies so far have investigated the potential importance of kernel function selection in the success of the SVMs classification [43,44]. In this context, our study explores the potential of a multi-temporal post-classification change detection scheme based on SVMs and TM imagery developed to obtain a LULC cartography and investigate the spatiotemporal dynamics throughout a period of 17 years (1993-2010) at the typical Mediterranean setting located in the Municipality of Mandra-Idyllia, Attica Region, Greece.

Study Area
The study area covers the area of the Municipality of Mandra-Idyllia, which includes the former  plantations and measure the expansion of Eastern Redcedar into the deciduous woodlands within the Forest ( [41,42], respectively). However, to our knowledge, very few studies so far have investigated the potential importance of kernel function selection in the success of the SVMs classification [43,44]. In this context, our study explores the potential of a multi-temporal postclassification change detection scheme based on SVMs and TM imagery developed to obtain a LULC cartography and investigate the spatiotemporal dynamics throughout a period of 17 years (1993-2010) at the typical Mediterranean setting located in the Municipality of Mandra-Idyllia, Attica Region, Greece.

Study Area
The study area covers the area of the Municipality of Mandra-Idyllia, which includes the former  The natural environment of Attica consists of a set of ecosystems dependent on each other that create the special conditions for the flora and fauna of the area. The Municipality of Mandra-Idyllia is part of the wider natural environment of the Attica region. The latter has high biodiversity, containing many endangered species of flora and fauna and including areas that have been designated as priority habitats by European legislation, in combination with the general local characteristics, and is called an "Attic Landscape". Regarding the morphology of the area, the main feature of the Municipality is the strong terrestrial slopes and the alternation of landscapes and the relief. The relief of the study area is characterized as mountainous in its western and northern part (600 to 800 m), as semi-mountainous to hilly in the zone of 100 to 600 m and as lowland in the coastal zone-southeast-as well as in the valleys of the rivers (0 to 100 m). Most of the study area is characterized by mild morphological inclinations of 2% to 17%. Morphological slopes greater than The natural environment of Attica consists of a set of ecosystems dependent on each other that create the special conditions for the flora and fauna of the area. The Municipality of Mandra-Idyllia is part of the wider natural environment of the Attica region. The latter has high biodiversity, containing many endangered species of flora and fauna and including areas that have been designated as priority habitats by European legislation, in combination with the general local characteristics, and is called an "Attic Landscape". Regarding the morphology of the area, the main feature of the Municipality is the strong terrestrial slopes and the alternation of landscapes and the relief. The relief of the study area is characterized as mountainous in its western and northern part (600 to 800 m), as semi-mountainous to hilly in the zone of 100 to 600 m and as lowland in the coastal zone-southeast-as well as in the valleys of the rivers (0 to 100 m). Most of the study area is characterized by mild morphological inclinations of 2% to 17%. Morphological slopes greater than 17% are observed in a significant part of the study area, on the slopes of the mountains. The slopes generally do not exceed 60%, except for some locations of local importance.
The main factors that contributed to the character of the existing vegetation in the wider region of Attica are long-term anthropogenic impact (logging, grazing, clearing, settlement creation, fires) and the climate [26,34]. In the wider region of Western Attica, there are three vegetation zones, which start from the sea and reach the tops of the mountains of Parnitha and Patera: the Mediterranean formations of the Eastern Mediterranean (Olea Ceratonion) start from the sea and cover the lowland and arid regions. This zone mainly occupies the areas with the lowest altitude. The Mediterranean formation of Aria (QuercionIlicis) is located above the previous zone, and occupies wetter areas. In this zone there is an uninterrupted spread of Aleppo Pine. The Mediterranean formation of the Cephalonian fir (Abietum Cephalonicae) occupies the peaks of the two mountains (Kitherona and Patera), while in Parnitha it forms a fairly remarkable forest. In this zone, where a short dry period and a higher amount of rainfall prevail, the main species and representative of the zone is the fir. Forests are the main form of land cover in the wider area (they occupy more than half of the total area-about 70%). The condition of the forests, where it has not been recently disturbed by fires, is at a good level. The grasslands in the wider area are smaller in extent. They are used as pastures, which is the case with most forest areas in Greece. In this area of Western Attica, grazing is literally out of control, as it is observed even in newly burned areas. Table 1 presents information concerning the spatial data (raster and vector) used to accomplish the current analysis.

Satellite Imagery
Freely distributed EO imagery was obtained from the Landsat-5 TM (Thematic Mapper) sensor dated 29-8-1993, 19-8-2001 and 12-8-2010 from the United Stated Geological Survey (USGS) archive (http://glovis.usgs.gov/). The images were obtained during the same month, ensuring consistency in seasonal vegetation. The selection of the imagery was based on the fulfillment of specific criteria, namely clear atmospheric conditions, high sun conditions, low water vapor, and temporal proximity to the fire event. The Landsat TM imagery was obtained already geometrically corrected and orthorectified. From these images, the thermal infrared band (band 6) was excluded. The bands of the image selected to be used in the study were 1-5 and 7. In addition to the TM imagery, the CORINE 2000 Land Cover (CLC) map (JRC-EEA, 2005) at a spatial resolution of 100 m for the site was obtained at no cost (from http://reports.eea.europa.eu/COR0-landcover/en).

Support Vector Machines (SVMs)
The SVMs proposed by [19,45,46] belongs to the category of supervised classification techniques, which work on the principle of statistical learning theory. During the process, the SVMs attempts to calculate the optimal hyperplane that separate the positive from the negative training points selected. For linear SVMs, the hyperplane with an offset b can be expressed as: where in the above equation x is a point lying on the hyperplane, w is a parameter determining the hyperplane orientation in space, and b represents the distance of the hyperplane from the origin. Each of the data points must fall on the proper side of the hyperplane, as follows: where the above 2 equations combined are expressed from the following inequality: During the training some data points are parallel to the optimum hyperplane designated as support vectors, and must satisfy the equation: Thus, the optimal hyperplane is generally achieved by minimizing → w, while satisfying Equation (4) from which it shown that this is happening when the margin is be equal to 2/ → w . However, in a different situation, in order to achieve maximum accuracy and higher dimensionality and to represent a complex hyperplane, the kernel functions K ( → χ ι · → w) are introduced by replacing the vector product in Equations (1) and (4). Another advantage with SVMs is the cost parameter C, which quantifies the penalty of misclassification errors. Apart from linear kernel, the other most popular kernel functions are polynomial, the Radial Basis Function (RBF) and the sigmoid kernels.

Methodology
This section provides details of the methodology implementation to satisfy the study objectives, whereas an overview of the approach as a whole is illustrated in Figure 2. For the data processing needs, the ENVI 5.1 remote sensing data processing software by EXELIS was used and for the mapping of the image processing products ArcGIS 10.1 by ESRI was used.

Landsat TM Pre-Processing
In deriving LULC maps from the obtained imagery, standard pre-processing procedures were subsequently applied. The Landsat TM imagery was imported into ENVI image processing software (ITT Visual Information Solutions SA) and the image georeferencing accuracy was initially checked with a reference map (1:10,000) that was available for the area. The positional accuracy was found to be within the level of the sensor pixel (RMS ~ 30 m) of the map, which was considered satisfactory for the purpose of the present study. Then, to assist a faster computational approach, a subset covering the study region was extracted from the image using the ENVI subset function. Before implementing the SVMs classifier, subsetting of the acquired imagery was performed for the region under investigation.

Classification
For the Landsat satellite data classification, nine different classes were used: Aleppo Pine Forest, Cephalonian Fir Forest, Forest Areas, Grassland Areas, Agricultural Areas, Burnt Areas, Bare Soils, Artificial Areas and Waters. The characteristics of Landsat images, and especially spatial resolution, determine the minimum mapping unit and, consequently, the mapping scale.
All six optical bands of the TM sensor were used. For training and testing of the SVMs model, as suggested by many researchers 10-30 p cases per class were selected, where p is the number of wavebands used [47][48][49]. Further, it is also observed that even with small training sets, SVMs gives reliable results [50]. In this study approximately 270 representative training pixels for nine classes were randomly selected from the Landsat TM image. Additionally, approximately 60 pixels for each class were also selected for validating the derived thematic maps. A CLC2000 map of the region was used for the selection of pure pixels. Before using the pixels for classification, class separability analysis was applied using the Jeffries-Matusita and the transformed divergence separability statistical measures (http://www.harrisgeospatial.com/portals/0/pdfs/envi/ENVI_User_Guide.pdf). The class separability test provides values in the range of 0 to 2.0, in which a value greater than 1.9 indicates a very good separability between the pixels chosen for classification. In the present study, spectral separability was performed by using six channels of the Landsat TM. The analysis indicated a class separability index of more than 1.58 for the six-band imagery, whereas for the validation points it was higher than 1.45. The lower separability value for the urban/industrial and the agricultural classes could be due to contamination of agriculture pixels by farm civil structures, which caused a

Landsat TM Pre-Processing
In deriving LULC maps from the obtained imagery, standard pre-processing procedures were subsequently applied. The Landsat TM imagery was imported into ENVI image processing software (ITT Visual Information Solutions SA) and the image georeferencing accuracy was initially checked with a reference map (1:10,000) that was available for the area. The positional accuracy was found to be within the level of the sensor pixel (RMS~30 m) of the map, which was considered satisfactory for the purpose of the present study. Then, to assist a faster computational approach, a subset covering the study region was extracted from the image using the ENVI subset function. Before implementing the SVMs classifier, subsetting of the acquired imagery was performed for the region under investigation.

Classification
For the Landsat satellite data classification, nine different classes were used: Aleppo Pine Forest, Cephalonian Fir Forest, Forest Areas, Grassland Areas, Agricultural Areas, Burnt Areas, Bare Soils, Artificial Areas and Waters. The characteristics of Landsat images, and especially spatial resolution, determine the minimum mapping unit and, consequently, the mapping scale.
All six optical bands of the TM sensor were used. For training and testing of the SVMs model, as suggested by many researchers 10-30 p cases per class were selected, where p is the number of wavebands used [47][48][49]. Further, it is also observed that even with small training sets, SVMs gives reliable results [50]. In this study approximately 270 representative training pixels for nine classes were randomly selected from the Landsat TM image. Additionally, approximately 60 pixels for each class were also selected for validating the derived thematic maps. A CLC2000 map of the region was used for the selection of pure pixels. Before using the pixels for classification, class separability analysis was applied using the Jeffries-Matusita and the transformed divergence separability statistical measures (http://www.harrisgeospatial.com/portals/0/pdfs/envi/ENVI_User_Guide.pdf). The class separability test provides values in the range of 0 to 2.0, in which a value greater than 1.9 indicates a very good separability between the pixels chosen for classification. In the present study, spectral separability was performed by using six channels of the Landsat TM. The analysis indicated a class separability index of more than 1.58 for the six-band imagery, whereas for the validation points it was higher than 1.45. The lower separability value for the urban/industrial and the agricultural classes could be due to contamination of agriculture pixels by farm civil structures, which caused a similar signature of the pixel. In addition to this, band ratio images were also used to visualize the area, so that pure pixels could be identified and selected to be included in training the SVMs.

In Situ Research and Classification
In order to determine the classes while applying the supervised classification method, an extensive in-situ data sampling was performed in places with different LULC. Landscapes of the representative classes that were similar to the defined classes were photographed. For field work, Google Earth images were used for sampling reasons. The major selected points that were photographed are shown below ( Figure 3). As mentioned, the methodology of the study was based on the supervised classification of Landsat TM images from 1993, 2001 and 2010, while the classification was based on the Corine Land Cover 2000 project and Anderson et al.'s [51] classification standards. According to Campell [52], the classification system to be followed in each study must be compatible with those used by other users and in previous applications. The classification of land cover forms of forest legislation as currently in force, the peculiarities of the configuration and diversity of the various forms of LULC constituted a new form of classification tailored to the needs of the study area. similar signature of the pixel. In addition to this, band ratio images were also used to visualize the area, so that pure pixels could be identified and selected to be included in training the SVMs.

In Situ Research and Classification
In order to determine the classes while applying the supervised classification method, an extensive in-situ data sampling was performed in places with different LULC. Landscapes of the representative classes that were similar to the defined classes were photographed. For field work, Google Earth images were used for sampling reasons. The major selected points that were photographed are shown below ( Figure 3). As mentioned, the methodology of the study was based on the supervised classification of Landsat TM images from 1993, 2001 and 2010, while the classification was based on the Corine Land Cover 2000 project and Anderson et al.'s [51] classification standards. According to Campell [52], the classification system to be followed in each study must be compatible with those used by other users and in previous applications. The classification of land cover forms of forest legislation as currently in force, the peculiarities of the configuration and diversity of the various forms of LULC constituted a new form of classification tailored to the needs of the study area. The following table (Table 2) lists the final classification, the classes' information, and the detailed description of these classes. In addition, the forest legislation as well as ecological criteria and the peculiarities of the study area were taken into account in the whole process. The conclusions The following table (Table 2) lists the final classification, the classes' information, and the detailed description of these classes. In addition, the forest legislation as well as ecological criteria and the peculiarities of the study area were taken into account in the whole process. The conclusions drawn from the process of unsupervised classification also played an important role in determining the final classification classes. Information classes were created for the classes with great diversity. Nine final classification classes were identified. Pure forests of Aleppo Pine that can yield forest products, as well as mixed stands with a dominant species of Aleppo pines.

Natural Renaissance (neophyte)
Dense pure clusters of Aleppo Pine in a new plantation, of existing mature forests or in a mixture with evergreen broadleaves.

Cephalonian Fir Forest
Pure forest of endemic Cephalonian Fir and mixed clusters with Cephalonian fir as the dominant species.

Forest
Sparse bushy forest areas of coniferous-broadleaf or mixed formations on degraded soils, native forest vegetation on fields and semi-natural areas.

Grassland
Areas with dry leafy vegetation on stony soils of uncultivated cultivation, as well as mixtures thereof with sporadic woody species.

Fallow
Abandoned fields with permanent or temporary abandonment of agricultural holdings.

Burnt Areas
Forest and rural areas after a recent fire.

Bare soil
Bare and unexploited areas with the revelation of the mother soil material, natural rocks, inactive quarries.

Artificial Areas
Areas with continuous or discontinuous urban fabric, industrialized areas, photovoltaic parks, road networks with asphalt coating, transport networks, waste disposal sites, etc.

Water
Sea, lakes and rivers.

SVMs Implementation
A number of studies with SVMs indicated that kernel selection is important for the performance of the SVMs classifier [44]. In order to select the best kernel functions in this study, four kernel functions were applied in the present study as given in Equations (6)-(9) below: Linear: Polynomial: Radial basis function: where γ is the gamma term in the kernel function for all kernels except linear, d is the polynomial degree term in the kernel function of the polynomial kernel and r is the bias term in the kernel function for the polynomial and sigmoid kernels. γ, d and r are user-defined parameters, as their correct definition significantly increases the SVMs accuracy solution.
For the best performance of the SVMs, optimization of each of the four kernels was provided, based on different combinations as recommended by users and checked according to classification accuracy as an objective function [53,54]. The kernel functions were optimized with respect to the penalty parameter, the pyramid levels and the classification probability threshold value.
As recommended, the penalty parameter was set to its maximum value (i.e., 100), which forces all pixels to converge to a class, with the pyramid parameter and classification probability threshold values of zero, so that the Landsat TM image gets processed at full spatial resolution and also restricts all image pixels to get one class label. The γ parameter was set to be equal to the inverse of the number of spectral bands (here equal to the value of 0.167 for six spectral bands of Landsat TM). For both the polynomial and sigmoid kernels, the bias was taken as equal to 1, with the degree of polynomial set to 2 for the polynomial kernel.

Accuracy Assessment
For accuracy assessment, well-known indices, such as error matrix (user's/producer's accuracy and omission/commission error), overall accuracy and kappa statistic [54,55], were used. Overall accuracy provides a summarized view of the overall classification performance in percentage (%). The kappa coefficient provides the value of deviation of reference data and the classifier (in this case the SVMs) versus the chance of agreement between the reference data and a random classifier. It varies from 0 and 1, with values greater than 0.80 indicating a good agreement.
The producer's accuracy indicates the efficiency of a classifier, while the user's accuracy expresses the probability that a pixel belongs to a given class and the classifier has labeled the pixel correctly into the same given class. The commission error (expressed as %) is derived by differentiating the user's accuracy (also expressed as %) from 100, whereas the omission error (%) is computed by differentiating the producer's accuracy (%) from 100.

Analysis of Band Ratio and Classification
The results of this study are considered satisfactory. In particular, for the classification needs of the study area after many tests and according to the relevant literature, a series of ratios for spectral bands was created in order to enhance the difference in reflectivity between the selected classification classes. This results in an increase in the size of the spectral separation (as a statistical indicator) between selected Regions of Interest (ROIs) for a specific input file. Research in the literature has shown that the spectral discrimination can be enhanced and refers mainly to the use of indices and ratios of band ratios. When these new variables are involved with image spectral bands in the image classification process, they significantly enhance the spectral resolution of selected ROIs [56,57]. Table 3 presents the initial variables and spectral band ratios used in this work.
These final images included all 14 variables from the initial 1, 2, 3, 4, 5 and 7 bands and the ratios of the Landsat TM's spectral bands (eight variables) for each image. These results mean that the performance of the distinction is very satisfactory and the stage of the supervised classification can be followed. The final pseudo-color images for the three reference years are shown in Figure 4.
The method of combining classes was used in the classes that had to be merged selectively, in the already sorted images. The basic classification classes were combined with their corresponding information while they were excluded from the final classification file in order to maintain the nine basic classes of the initial classification.
Subsequently, in order to produce final thematic maps for each year of study (1993,2001 and 2010), a method of clearing the isolated pixels was followed. Aggregation is a useful cleanup process after sorting. The three maps created for the three time periods of the study area are presented in Figure 5.
14 Clay soils b5/b7 These final images included all 14 variables from the initial 1, 2, 3, 4, 5 and 7 bands and the ratios of the Landsat TM's spectral bands (eight variables) for each image. These results mean that the performance of the distinction is very satisfactory and the stage of the supervised classification can be followed. The final pseudo-color images for the three reference years are shown in Figure 4. The method of combining classes was used in the classes that had to be merged selectively, in the already sorted images. The basic classification classes were combined with their corresponding information while they were excluded from the final classification file in order to maintain the nine basic classes of the initial classification.  Subsequently, in order to produce final thematic maps for each year of study (1993,2001 and 2010), a method of clearing the isolated pixels was followed. Aggregation is a useful cleanup process after sorting. The three maps created for the three time periods of the study area are presented in Figure 5.

Estimation of Classification Accuracy
The overall classification accuracy for images of 1993, 2001 and 2010 was 89.85%, 91.01% and 90.24%, respectively. The statistical coefficient (K) was 0.96, 0.89 and 0.99, respectively, which indicates an excellent agreement between classification and reference data. The differences between

Estimation of Classification Accuracy
The overall classification accuracy for images of 1993, 2001 and 2010 was 89.85%, 91.01% and 90.24%, respectively. The statistical coefficient (K) was 0.96, 0.89 and 0.99, respectively, which indicates an excellent agreement between classification and reference data. The differences between the K factor and the total accuracy are due to the fact that the exported errors in the table are handled differently. The results of the supervised classification with the classifier of the SVMs are presented below per year (Table 4):

Change Detection and Confusion Matrix Analysis
One of the primary challenges in the science of dynamic terrestrial ecosystems is making valued predictions regarding the changes that may occur in the next decades of the century and the physical relationship with human activity. With the contribution of remote sensing in this field, it is possible to observe and monitor the earth's environment as a whole and on a continuous scale. Using the ArcGIS 10 (ESRI, Redlands, CA, USA) software and ENVI 5.1 within ArcGIS 10 environment, the raster files created using the ENVI 5.1 software were transferred to the ArcGIS 10.1 environment in order to conduct the change detection assessment. For each period of the study, the thematic maps of land changes and the maps of the change of the type of LULC were created, and are shown in Figure 6.
A confusion matrix for year 2010 is provided as an example to understand the classifier performance ( Table 5). The columns of the table represent the ground truth classes, while rows represent the classes of the classified image to be assessed. Cells show the number of pixels for all the possible correlations between the ground truth and the classified image. Diagonal cells contain the number of correctly identified pixels. The results of 2010 indicate that in the case of fir, pine and burned areas, SVMs shows a good performance in correctly identifying the pixel followed by crop lands, brush lands, regeneration lands and bare soils. The least correctly identified pixels belong to the of class urban areas. In classified images, only the major classes are shown to reduce salt and pepper effects. Water bodies are not detected in the area. The next table (Table 5)  predictions regarding the changes that may occur in the next decades of the century and the physical relationship with human activity. With the contribution of remote sensing in this field, it is possible to observe and monitor the earth's environment as a whole and on a continuous scale. Using the ArcGIS 10 (ESRI, Redlands, California, CA, USA) software and ENVI 5.1 within ArcGIS 10 environment, the raster files created using the ENVI 5.1 software were transferred to the ArcGIS 10.1 environment in order to conduct the change detection assessment. For each period of the study, the thematic maps of land changes and the maps of the change of the type of LULC were created, and are shown in Figure 6.

Discussion
Our findings indicate that SVMs is a promising technique for land cover classification. Some other studies also indicated that SVMs is a powerful classification scheme. In the study by Srivastava et al. [54] the authors concluded that SVMs is performing well compared to other techniques in this field, and to get a best performance from SVMs, kernel optimization is very important. A study by Kavzoglu and Colkesen [25] also mentioned that SVMs' accuracy may show variation depending on the choice of the kernel function and its parameters. There results indicated that a radial-based function performed better than other kernel functions, which is also observed in this study. A similar study by Singh et al. in 2014 [58] also indicated that SVMs is powerful in LULC classification. They used a hybrid approach in combination with SVMs to provide the land cover of a mangrove ecosystem and revealed a better performance of the SVMs. Sukawattanaviji and Chen [59] demonstrated the use of SVMs on fused SAR and optical images and compared it against a Maximum Likelihood (ML) classifier; they also found that SVMs performs better than ML. SVMs was applied by Shao and Lunetta [60] for land-cover characterization using MODIS satellite data. They also compared the results with neural networks and classification and regression trees algorithms and concluded that SVMs outperforms both the algorithms. The selected classifier, SVMs, along with other non-parametric classifiers, such as Artificial Neural Networks (ANN), have been proven as being among the strongest classifiers [45,46]. Based on the thorough analysis of SVM's high performance in most of the cases, it was implemented here for land cover classification.
From the obtained data it is concluded that in 1993, forest areas were the dominant class, covering 41. The grassland areas in 1993 covered an extent of 45.810 km 2 (10.73%); in 2001 they covered 39.59 km 2 (9.27%) and finally in 2010 they appeared to occupy 4.52%, i.e., an area of 19.29 km 2 . This class, which has been observed mainly in forests and forest areas, including mountain pastures, is one of the forest ecosystems protected and managed by the Forest Service. The land use it represents can evolve depending on the dynamic natural evolution of ecosystems into forest ecosystems. However, these areas can be attributed to primary production, such as animal husbandry, which is traditionally occupied by the inhabitants of the area, but also to agriculture. From these results and from the reduction of primary production in the study area we can assess that the main reason for the decrease in the spread of these areas is the reduction of sheep grazing, which was a traditional occupation of the region's coastal populations for hundreds of years, during the examined time reference period of the study area.
Firs are a very important ecosystem. At the highest altitudes of Mount Kitheron, Cephalonian Fir trees are developed. The catastrophic fires of 2009 and 2010 did not harm the fir forests of the area, which were saved. Nevertheless, in these cases, wildfire assessment can be used together with the current methodology for retrieving information on a burnt area, by using remote sensing data [61][62][63][64][65]. The class of fir forests in the years 1993 and 2001 shows a similar coverage of about 1% (0.9% and 0.87, respectively) occupying an area of about 3.8 km 2 . In the year 2010 the coverage reached 1.25%, occupying an area of 5.33 km 2 . In other words, there was an increase of 0.35%. This is very important due to the lack of a good road network (low accessibility and low impact from human activity). Bare areas show a downward trend: in the year of 1993 they covered 2.09% and occupied an area of 8.91 km 2 , in 2001, 0.68% and an area of 2.89 km 2 and in 2010, 0.58% and an area 2.47 km 2 . Artificial areas throughout the research period remained stable and remained, for all three reference years 1993, 2001 and 2010, at 3%. These data were confirmed by the pilot project of the Municipality of Mandra-Idyllia in 2013, in which urban areas were found to consistently cover, from 1993 to 2013, 3% of the total area. This is due to a lack of incentives and investments for creating new businesses, which would boost the use of the industrial area of the Municipality of Mandra-Idyllia and thus attract new residents and increase the urbanization of the Municipality. These results are confirmed by the 2001 and 2011 population censuses (based on the published census results of the Hellenic Statistical Authority), which show a 12.37% decrease in the population, and by the 2006 statistics as described in the business plan of the Municipality of Mandra-Idyllia, according to which secondary production showed decreasing trends while the unemployment rate was very high. Additionally, the creation of infrastructure in the adjacent areas of the Municipalities of Aspropyrgos and Elefsina from 1997 onwards (Hospital, Olympic projects, etc.) further burdened the area. Besides that, in these areas there was an increase in the population as well. Agricultural land has been growing steadily over time. The coverage ranges from 20.22% and an area of 86.32 km 2 in 1993, to 21.60% and an area of 92.24 km 2 in 2001 and 24.84% and an area of 106.08 km 2 in 2010. The increase in agricultural land, which was higher after 2001, is due to the lack of industrial use in the region. The shift of the population to primary production was a necessity, as there was low supply in the secondary and tertiary production sectors, to which the period of economic crisis, which started in 2007, also contributed dynamically, as no new infrastructure was developed and development remained stagnant. The shift of the population to agricultural production is due to the large agricultural land that the Municipality has to offer for the productive process of the population in order to find financial resources and survival.
Overall analysis indicates that due to the expansion of agricultural uses over semi-natural, forest and grassland areas does not seem to take place for the production of goods for consumption, since the available statistics indicate a decline in agricultural activity. It is obvious that the strictest legal framework for the protection of forest ecosystems, which is becoming even stricter on reforestable areas after fires and their legal documented protection against goat herding, has been a strong incentive for the preservation and development of forest ecosystems. The region of the capital, which is huge in size and infrastructure classes, strongly creates the need for its redefinition and the creation of a new design. Expensive funding for the uncontrolled increase in size of the city of Athens has worked and continues to work to the detriment of projects of a purely developmental nature that could take place in the region (dams, irrigation projects, interregional highways, industrial areas etc.). The relationship between the creation of huge cities and infrastructures is interdependent, as the creation of infrastructure in large cities strengthens the tendencies of population concentration and economic activities in them. A typical example in the study area is the creation of the suburban railway, which left out the areas of the Municipality of Mandra-Idyllia, resulting in population and productive increase of the neighboring municipalities of Aspropyrgos and Elefsina.

Conclusions
This work provides evidence of how geoinformation technologies, and in particular remote sensing data and geospatial data analysis techniques along with SVMs can be coupled for assessing landscape changes through significant time periods. Specifically, in the presented research work the statistics of the intertemporal changes and the classification of the satellite images of the time series 1993, 2001 and 2010 provided useful conclusions about the most basic types of land use, which are agricultural use, forest use and grassland. By delving a little deeper into the changes of coverage types of interest, that is, in urban uses, in cultivated land (pastures and forests) and semi-natural areas, several interesting conclusions are drawn.
In general, the areas linked to the activities of the primary sector occupy a very large percentage of the total area of the Municipality of Mandra-Idyllia, a fact that is reduced throughout Greece. The character of these areas, most of which are linked to the activities of the primary sector and due to the fact that the Municipality of Mandra-Idyllia is the closest and largest Municipality near Athens, makes it more vulnerable to the economic consequences of urbanization. Although remote sensing techniques require advanced scientific data processing, they help in saving time. The results obtained from this process, including the statistical results and the resulting trends, may support decisions related to the environment and the economy or even address the need for a new strategic policy for the development of the Municipality in the future. It can also be a dynamic tool in the study of the environment in order to meet future needs for protection and sustainability. Furthermore, this research study can become an important source of information for any administrator who wishes to carry out a development project in the wider area of the municipality, considering it as a dynamically evolving but unused ecosystem. Finally, the main outcome and key conclusion of this research is that the use of remote sensing data and co-analysis with SVMs offers an efficient integrated methodological framework for addressing issues such as intertemporal LULC changes and similar environmental alterations. This change detection framework can become an optimum tool for making important and informed land management decisions.
Funding: This research received no external funding.