Comparison of Pixel-Based Classification Algorithms Using Landsat-8 OLI and Sentinel-2 MSI for Land Use/Land Cover Mapping in a Heterogeneous Landscape

1. Introduction

Land cover represents the characteristics of earth surface shaped by various natural agents or anthropogenic interventions. From earth-observation perspective, the term “land cover” defines the land types (i.e. vegetation, water bodies, crops, built-up areas, etc.) which can be detected from a distance. Land cover is a critical variable for earth surface studies since it can be changed over time [1]. On the contrary, the term “land use” refers to the way a particular land is used involving the associated economic purpose of this use [2]. Both concepts are interrelated. For example, a land cover type such as a forest may support a series of land uses (e.g. timber production, recreation, rangeland, etc.) while a land use such as agroforestry may include a series of land cover types (e.g. forests, plantations, annual crops, etc.). In this research, they are used complementary to depict all kind of existing land cover or land use types within the study area.

Land use/land cover information is essential for management and monitoring of natural resources, modeling, spatial planning, land administration and sound decision-making. Satellite-based classification provides land use/land cover spatial-explicit information and map generation at global, national or regional scales. However, medium-resolution (Landsat-like, 10–30 m) is more adequate to detect most human–nature interactions [3]. The opening of the Landsat archive in 2008 [4] and the launch of Sentinel-2 in 2015 provide optical multispectral imagery data at medium-high resolution. Free and open access policy on this imagery promoted the development of new products and applications across space and time, especially in the domain of land use/land cover mapping. This data policy combined with the increase of computing power and concurrent reduction of costs, has facilitated large area mapping and expanded the number of users worldwide [5].

Due to the complexity of land use/land cover characterization, several studies have been mainly focused to methods for mapping a single land cover type (i.e. forests, wetlands, forest fires, agriculture, urban areas or water). For example, [6,7] for forests, [8] for urban areas, [9,10] for croplands, [11,12] for wetlands, using either Landsat or Sentinel imagery. However, multiple-class characterization required for simultaneous and spatially exhaustive mapping [13]. Thus, effective and efficient methods are required for satellite imagery classification to provide meaningful information regarding all land use/land cover within a specific area.

A variety of classification approaches (unsupervised, supervised, parametric, non-parametric, object-oriented) has been developed and applied to derive land cover information with different degree of success. Per-pixel classification approaches remain the most popular in the analysis of satellite-derived imagery [14]. Here, we used supervised per-pixel classifications for multiple land use/land cover types mapping.

In supervised approaches, reference data are required to characterize the variability of land cover across space and time and serve as reference dataset for training and validating classification models. A suitable reference data is a fundamental requirement in supervised image classification [15]. We use existing authoritative geospatial datasets of higher accuracy as a pool for training and validation. The reference datasets spans forestlands, cultivated fields, discontinuous urban fabric, built-up areas, and wetland habitat types. The classification scheme of land cover classes is based on Copernicus Land Cover (CLC) nomenclature [16]. Based on CLC2018 land use/land cover distribution, a stratified random sampling scheme is deployed to train the classifier and access classification accuracy. Classification accuracy depends on the satellite imagery, the classification algorithm being used, and the nature of training data as well [17].

Four popular classifiers ML, RF, KNN and SVM selected, and their implementation in Erdas Imagine 2020 was used to run the experiments. Description and analysis on these classifiers can be found in the literature. For example, [18] for Bayesian classifiers; [19] for SVM; [20] for KNN; [21] for Random Forests.

Maximum likelihood method is included in our research due to its wide application and use in commercial image-processing software [22]. On the other hand, the above machine learning algorithms have gained great attention for classifying land use/land cover types in the last decade.

In recent evaluations, SVM and KNN, with the exception of Naïve Bayes (a Maximum Likelihood variant) performed similarly in per-pixel classification of 26 Landsat TM imagery 10kmx10km blocks [23]. In a peri-urban and rural with heterogeneous land cover area in Vietnam, SVM produced the highest overall accuracy (OA) using Sentinel-2 MSI, followed consecutively by RF and KNN. However, all three classifiers showed a similar and high OA (over 93.85%) when the training sample size was large enough (>750 pixels/class) [24].

In this study, we test the above mentioned classifiers to derive land use/land cover information. We explore classification performance under six different scenarios. The study investigates the performance of the above classifiers using Landsat-8 OLI and Sentinel-2 MSI across a heterogeneous Mediterranean watershed, based on the same available land cover reference, training and validation data per sensor. In case of Landsat, spectral indices (EVI and NDMI) were included. These indexes have been reported in the literature that improve land use/land cover classifications accuracies [25]. We evaluate classification performance for an area with complex landscape and investigate how single date and seasonal optical multispectral data impact land use/land cover classification accuracy.

2. Study Area

The test site is Mygdonia basin, which is located in Northern Greece. Its watershed covers an area of 190,285ha. It lies East of Thessaloniki city (40˚40΄56,49N and 23˚18΄21,15E, WGS84) at a distance of 41,6 km (Figure 1). The watershed is surrounded by mountains in the North (Mount Krousia) and in the South (Mount Cholomontas), by hills in the West and by Rentina Gorge and Kerdylia mount in the East. Its elevation ranges from 35m to 1,129m. At the center of the watershed, there are two lakes (Koronia and Volvi). The watershed is drained through seasonal and intermittent streams at these lakes.

The two lakes along with their surrounding wetlands are listed as a Wetland of International Importance by the Ramsar Convention since 1975 (GR005: 16,388ha). Along with the valley of Rentina Gorge, they have been designated as Special Conservation Zones within the Natura2000 network (GR1220001 and GR1220003: 28,734.90ha) in 2017. These protected sites constitute a unique complex of interconnected natural ecosystems of lakes, seasonal streams, channels, riparian forests, shrubs, wet meadows and fields.

At the southernmost end of the watershed, on the slopes of mount Cholomontas, there is a portion of another protected area (GR1270001). It has an area of 15,651.14ha dominated by beech, oak and pine forests.

Non-irrigated arable lands are distributed across the watershed up to the productive forests in the south. Intensively cultivated lands, mainly irrigated, surround the wetland ecosystem. According to CLC2018, 49,28% of the watershed is under agricultural use while forestlands occupy 42,04%, water 5,46%, discontinuous urban fabric 1,31%, wetlands 1,17%, and developed areas only 0,33%. Most of the land under agricultural use is used as cropland (93,770 ha) while the area of perennial crops such as fruit and olive tree plantations and vineyards as well, account for only 0,65%. Irrigated lands cover 14,27% while non-irrigated 50,97% of croplands. Approximately 30,98% of forestlands are broadleaf forests, 2,73% are pine forests, 12,04% mixed forests, 36,80% shrubs and 10,10% transitional woodlands (Figure 2).

The climate is considered temperate (Csa-Mediterranean mainland) with warm and dry summers and cool winters. The mean annual temperature is 22.6°C in summer and about 4°C in winter. The mean annual rainfall is 593mm according to the records of the last 40 years.

3. Materials and Methods

3.1. Satellite Imagery

Landsat-8 Operational Land Imager (OLI) surface reflectance (C2L2) data were obtained from the United States Geological Survey website [26]. Two scenes (path/row: 183/032 and 184/032) required to cover entirely the study area. Following a search, cloud free (<10%) scenes were carefully selected for summer and winter seasons. Acquisitions dates were 01 July 2018 and 22 June 2018 for dry season and 28 January 2020 and 17 February 2019 for winter season. A mosaic contained six bands (blue, green, red, near infrared (NIR), shortwave infrared (SWIR 1, SWIR 2) was created at the study area limits

Sentinel-2 (L2a) MSI imagery downloaded from the Sentinels Scientific Data Hub [27]. Each product consists of 100x100 sq. km orthorectified granules or tiles. Four cloud free (<10%) granules required to cover entirely the study area, sensed in summer 2018, were selected (Table 1).

The 13 spectral bands of Sentinel-2a span from the visible to SWIR spectrum, at 10m, 20m and 60m spatial resolutions. The bands at 60m spatial resolution are dedicated primarily for detecting atmospheric features. Therefore, they have been excluded from the analysis [28]. A mosaic of ten bands (2-8, 8a, 11 and 12) was created at the watershed limits. Nearest neighbor interpolation was employed to downscale the spatial resolution of 20m bands at 10m. This process has been shown to perform very satisfactory compared to other approaches [29]. Both Landsat 8 OLI and Sentinel-2a image scenes are spatially registered to Universal Transverse Mercator (UTM)/World Geodetic System 1984 (WGS84) projection.

3.2. Land Cover Reference Data

A series of existing land cover reference data were retrieved from existing national databases (Table 2).

Land use/land cover types in 2018 for the entire watershed were obtained from European Copernicus Program (Corine Land Cover product-CLC 2018). In Europe, CORINE Land Cover (CLC) provides harmonized and comprehensive maps of land cover and land use change at European level [30 Buttner, 2014]. The program was established by the European Commission (EU) in 1990 for facilitating policy making at European level. The most recent CLC2018 comprises of 44 thematic classes at the third level with a minimum mapping unit (MMU) of 25 Ha for areal features, and 5 Ha for changes, respectively. It is an excellent tool for strategic analysis and planning at European level. However, CLC’s thematic content comprises a mixture of land cover and land use classes. In addition, its MMU serves well the needs of the European Union but is not suited for national or local detailed land use/land cover mapping [31].

Information regarding plantations and vineyards either irrigated or non-irrigated retrieved from the Land Parcel Identification System (LPIS). However, these data refer only to parcels for which there are individual claims for subsidies made by farmers and receive European Union Aid [32]. Therefore, they do not represent the entire number of cultivated fields within the entire watershed.

Information on habitat types acquired through the national large-scale Natura2000 database. We retrieved spatial-explicit information on the habitat types and vegetation species dominated the wetland.

Forestlands retrieved from the Forest Map national program. Forest Map is a very-high-resolution diagram at the scale of 1:5,000, depicting forests and non-forests, according to the current legislative framework of Greece [33]. Furthermore, we obtained available forest management plans, from the Hellenic Forest Service to retrieve information on (co)dominant forest species at the stand level and land use types within managed forests. However, forestlands within the available plans cover only a portion of forests equal to 26,857ha (28% of forestlands, according to the Forest Map). Forestlands outside plans are mainly unmanaged of different structure and crown cover, distributed across the watershed and comprised of degraded broadleaf forests (mixed or not), evergreen shrubs, and reforested pine forests.

3.3. Sampling Design

A stratified random sampling design is adopted. Copernicus CLC18 product used as the basis for sampling units distribution across all identified classes. Based on these classes, the sample size estimated to be equal to 2,356 for a required thematic accuracy of 85%. The samples distributed randomly, proportional to CLC2018 class area (Table 3).

3.4. Training Data

Based on the above sample distribution, sampling plots were defined at the pixel spatial resolution (30x30m) of Landsat imagery. Each random point was located at the center of the respective pixel using a gridded fishnet on Landsat-8. Each plot was divided into 3x3 pixels to coincide with Sentinel-2a spatial resolution (10m).

Land cover reference data were processed to generate the following seven thematic datasets with the highest accuracy: Based on the Forest Map, we excluded non-forest areas and created a dataset exclusively for forests. In areas, where the forest dataset overlapped with forest management plans, we extracted information on forests (brooadleaf, needleaf or shrubs) at the stand level, based on the dominant species. In wetlands, we excluded forestlands based on the above forest dataset. Then, we create a natural habitats (inland marshes, shrubs, wet meadows and high reeds) dataset, excluding all other land use/land cover types, based on their unique Natura2000 database 4-digit codes. Discontinued urban fabric (small towns and villages) areas extracted from urban zones provided by the Forest Map. A dataset regarding roads and built-up areas generated by processing the cadastral database. Plantation trees (olive, fruit and forest ones) and vineyards extracted from the LPIS database. The last generated dataset consist of crops either irrigated or non-irrigated.

In addition, Google Earth high-resolution (2019) imagery was used for visual interpretation of each plot, based on physiognomic attributes (color, shape, size, pattern and texture). This orthoimagery was the closest existing one to satellite imagery acquisition dates.

Based on the above interpretation, cross-referenced by each of the thematic datasets, each plot was assigned a land use/land cover unique type. Thus, a consistent large database was generated for selection of training data.

3.4. Classification

Four classification methods were applied, one parametric (ML) and three machine learning classifiers, KNN, RF and SVM. All procedures in this study were implemented using the Erdas Imagine 2020 commercial software.

We tested the utility of single-dated (summer 2018) and combination of summer-winter spectral bands of Landsat-8 OLI and Sentinel-2a MSI as data input, developing six different scenarios. Two of them refer to the use of spectral indices (EVI and NDMI) with single-date and seasonal Landsat-8 OLI spectral bands. EVI is sensitive to vegetation intra-annual variations while NDMI is sensitive to moisture content. They both used for different types of vegetation and irrigated fields discrimination. We acknowledged that many other different combinations of spectral and temporal features or approaches could be used. We decided to limit our research to the aforementioned features in our analysis.

In the initial phase, we tested numerous iterations of classifications with different combinations and number of land use/land cover classes on both types of imagery. However, classification performance was unacceptable (<85%) for over ten classes in both types of imagery. Land cover classification accuracy is affected by the number of classes identified. Overall classification accuracy decreases by increasing the number of classes [34].

Therefore, we adopted a classification scheme of 9 classes using Sentinel-2a and 8 classes using Landsat-8 OLI (Table 4). Rare and small-sized classes either grouped to form new classes or integrated to existing large enough classes. For example, vineyards, fruit and olive trees integrated to arable non-irrigated lands. Mineral extraction sites, discontinued urban fabrics, and industrial/commercial classes form a new class entitled as “Artificial Surfaces”. Classes such as, land principally occupied by agriculture, complex cultivated patterns, were deleted. They are generic land use/land cover types and include several other land types.

During the process, we selected 3,535 training data (polygons) for Landsat-8 OLI and 2,753 for Sentinel-2a MSI classification. We defined a set of training polygons by random sampling 70% of the points selected for class validation. This data is then used to train supervised classification algorithms. The remaining 30% of samples was used for classification validation(Table 4).

Training polygons were manually generated at random locations of sample plots taking in account that cover types should be spectrally homogeneous. For this reason, in many cases, we forced to generate training polygons away from sample locations. We avoided long and thin training polygons. Small polygons tend to be prone to edge effect. Moreover, we selected more training polygons in areas where land cover reference data was missing or in highly heterogeneous areas, in order to increase classification accuracy. The generation of training data in areas where land cover reference data are missing proved to be an issue. Their selection was based on our expert knowledge of the study area in relation to spectral data.

3.5. Accuracy Assessment

All classifiers were tested on the entire watershed based on the same training and validation data per sensor. We evaluated classification performance using the Overall Accuracy (OA), Producers’ Accuracy (PA), Users’ Accuracy (UA) and Kappa Coefficient (KC). For accuracy assessment, we selected 1,260 validation points for Landsat-8 OLI and 782 points for Sentinel-2A MSI.

3.5.1. Landsat-8 Scenarios

In scenario 1, we used single date Landsat-8 imagery acquired in summer. The ML classifier produced a slight higher overall accuracy (91.67%) comparing to machine learning classifiers (Table 5). In terms of class accuracy, the best results for ML achieved for water bodies, non-irrigated arable land and shrubs. However, needleleaf forest and high reeds have the lowest user accuracy and K-coefficient. Needleleaf forest is confused with broadleaf forest in mixed-forest areas. Moreover, unsuccessful reforested pine forests are confused with tall sclerophylous vegetation (evergreen shrubs).

Very dense high reeds in wetlands are confused with broadleaf forest in the hillsides showing similar spectral signature. The low producer accuracy reported in high reeds class is also an issue for all machine learning classifiers.

In scenario 2, we used single date Landsat-8 OLI combined with EVI and NDMI in-dices acquired in summer as well. Again, the ML classifier produced the highest overall accuracy (90.95%) over machine learning classifiers (Table 6). The contribution of indices in the classification process is not sufficient. The performance of all machine classifiers is similar and close to ML overall accuracy.

In scenario 3, we used multi-dated Landsat-8 OLI (summer and winter). For the first time, all classifiers have an overall accuracy slightly over 90% (Table 7). In this scenario, KNN produced the highest overall accuracy (91.90%) followed by ML, SVM and RF (Table 7). Accuracy in all classes is improved except those of needleleaf forests and high reeds.

In scenario 4, we used multi-dated Landsat-8 OLI (summer and winter) combined with the respective EVI and NDMI indices. Overall accuracy of all classifications improved but show similar results. The KNN classifier produced the highest overall accuracy (93.02%) (Table 8).

3.5.2. Sentinel-2A

In scenario 1, we used single-dated Sentinel-2a imagery. RF classifier produced slightly higher OA (93.86%) compared to ML and KNN (Table 9) which have similar results. SVM classifier has the lowest overall accuracy. The best results achieved for water bodies, broadleaved forests, needleleaf forests and non-irrigated arable lands. The lowest class producers’ accuracy is observed for roads in SVM and high reeds in KNN. However, high reeds class accuracy was low in SVM as well.

In scenario 2, we used multi-dated Sentinel-2a imagery (summer and winter). The ML classifier produced the highest overall accuracy (96.68%). Machine learning classifiers produced lower but similar results (Table 10). We observed that roads class is confused with artificial surfaces class, especially within urban areas. In addition, the class of shrubs is confused partially with non-irrigated lands, whereas small-sized fields are surrounded by shrubs.

4. Results

Table 11 shows the obtained overall accuracy per classifier and sensor in each developed scenario. In reference to Landsat-8 OLI classification, KNN was the best classifier in scenario 4, achieving the highest OA=93.02%. In this case, KC reached the highest value (0.9227). Under the same scenario, the ML classifier reached the second highest OA=92.46%, followed by SVM with OA=92.30%, and RF with OA=91.67%.

In scenario 4, all classifiers achieved the highest OA among all previous scenarios with minor variations (<1%) due to sufficient training data. Thus, we concluded that the use of multi-dated multispectral seasonal Lansdat-8 OLI data combined with spectral indices increases the performance and overall accuracy of classification for land use/land cover mapping. However, the contribution of spectral indices (EVI and NDMI) in classification performance was not significant (+1%) in all scenarios. The resulted classified maps is presented in Figure 3.

In reference to Sentinel-2a MSI classification, the ML classifier has the highest OA=96.68% under scenario 2, when intra-annual seasonal multispectral data is used (Figure 4). It is observed that the OA of ML is higher than the OA of RF (OA=93.86%) which ranked first in scenario 1. In scenario 2, the SVM classifier produced the highest accuracy (OA=92.84%), followed by KNN (OA=92.71%) and RF (OA=91.82%). We observed a major variation amongst ML and machine learning classifiers performance. The OA of ML is higher (>4%) than machine learning classifiers.

5. Conclusions

The aim of this work was to analyze the performance and accuracy of different classification classifiers (ML KNN, RF and SVM) and evaluate Landsat-8 OLI and Sentinel-2A MSI imagery for the identification and mapping of land use/land cover types in a highly heterogeneous Mediterranean site.

Based on the results, the achieved overall classification accuracies for both satellite imagery was acceptable (>85%) and the performance of selected machine learning classifiers was quite similar and statistically not significant in all scenarios. The best performing classifier is ML using seasonal bi-temporal Sentinel-2a imagery (OA 96.68%). This OA is higher than the best performing classifier KNN (93.02%) using seasonal bi-temporal Landsat-8 OLI imagery combined with spectral indices EVI and NDMI. Following a visual assessment of the respective classified map by the ML classifier, we realize that it overestimates artificial surfaces (Figure 4). Artificial surfaces are confused with high reflectance bare soils or open areas with no vegetation cover across the watershed. However, this finding cannot be supported by the respective confusion matrix (Table 12).

As far the performance of machine learning methods, we believe that more training data are required in case of Sentinel-2a classification. Machine learning methods require enough training samples to make optimum decisions [23]. However, the high spatial variability and spatial structure of the study area (small-sized area and sparsely land cover classes) affects the selection of proper training data.

In terms of class user’s accuracy, the lowest accuracy observed for High reeds (80%) followed by Needleleaf forests (88,37), according to confusion matrix in scenario 2 (Table 12). High reeds are commonly confused with shrubs in the wetlands, whereas they form mixed associations with shrubs (Tamarix sp.). In areas of unsuccessful reforestation with pine trees, needleleaf forests are confused with evergreen shrubs.

In reference to Landsat-8 (scenario 4), the lowest class users’ accuracy observed for Broadleaf forests (84,73%). This can be explained by the existing mixed formations of low in height broadleaf trees with shrubs, which span the watershed. PA’s lowest accuracy 80,33% observed for Needleleaf Forests and 82,76% for High Reeds.

We would recommend the classifier KNN using Landsat-8 imagery for land use/land cover mapping. However, we would be prudent with the application of the ML classifier using Sentinel-2 imagery. Machine learning algorithms are stable and produce similar results in all scenarios.

This paper presents a methodology for testing different classifiers for land use/land cover mapping of a high heterogeneous and complex landscape. It also includes a process for processing available geospatial databases, and a multi-source training data preparation. Integration of intra-annual temporal-spectral data into classification produces land use/land cover maps of high accuracy. This study represents an important step toward multiple-class land use/land cover mapping using spectral-temporal Landsat-8 or Sentinel-2 features by providing a quantitative assessment on classification accuracy. Our work contributes to the evaluation of classification algorithms for updating Copernicus Land Cover product. It documents that major classes at the 3rd level of Copernicus nomenclature, such as urban fabric, roads, irrigated and non-irrigated lands, broadleaf or needleleaf forests, shrubs, water bodies, large enough streams and wetlands vegetation can be classified with high accuracy based on seasonal multispectral data.

We believe that more research is required in the domain of land use/land cover mapping. Future research should be oriented towards the development of novel methods by integrating ancillary geospatial data or by integrating time-series spectral-temporal data into a classification model for land use/land cover mapping.

Table 13. Confusion matrix for scenario 4 (Landsat-8 KNN classification).

	Artificial surfaces	Broadleaf forest	Needleleaf forest	Non-irrigated arable land	Permanent freshwater lakes	Permanently irrigated land	High reeds	Shrubs	Total	PA (%)	UA (%)	KC
Artificial surfaces	88	0	0	6	0	0	0	0	94	95.65	93.62	0.9311
Broadleaf forest	0	172	7	0	0	2	4	18	203	91.01	84.73	0.8203
Needleleaf forest	0	2	49	0	0	0	0	0	51	80.33	96.08	0.9588
Non-irrigated arable land	3	1	0	350	0	0	0	8	362	95.11	96.69	0.9532
Permanent freshwater lakes	0	0	0	0	56	0	0	0	56	100.0	100.0	1.0000
Permanently irrigated land	0	2	0	0	0	78	0	2	82	95.12	95.12	0.9478
High reeds	0	0	0	0	0	0	24	0	24	82.76	100.0	1.0000
Shrubs	1	12	5	12	0	2	1	355	388	92.69	91.49	0.8778
Total	92	189	61	368	56	82	29	383	1260	93.02		0.9109

Table 13. Confusion matrix for scenario 4 (Landsat-8 KNN classification).

	Artificial surfaces	Broadleaf forest	Needleleaf forest	Non-irrigated arable land	Permanent freshwater lakes	Permanently irrigated land	High reeds	Shrubs	Total	PA (%)	UA (%)	KC
Artificial surfaces	88	0	0	6	0	0	0	0	94	95.65	93.62	0.9311
Broadleaf forest	0	172	7	0	0	2	4	18	203	91.01	84.73	0.8203
Needleleaf forest	0	2	49	0	0	0	0	0	51	80.33	96.08	0.9588
Non-irrigated arable land	3	1	0	350	0	0	0	8	362	95.11	96.69	0.9532
Permanent freshwater lakes	0	0	0	0	56	0	0	0	56	100.0	100.0	1.0000
Permanently irrigated land	0	2	0	0	0	78	0	2	82	95.12	95.12	0.9478
High reeds	0	0	0	0	0	0	24	0	24	82.76	100.0	1.0000
Shrubs	1	12	5	12	0	2	1	355	388	92.69	91.49	0.8778
Total	92	189	61	368	56	82	29	383	1260	93.02		0.9109