Land cover and land use classification performance of machine learning algorithms in a boreal landscape using Sentinel-2 data

ABSTRACT

In recent years, the data science and remote sensing communities have started to align due to user-friendly programming tools, access to high-end consumer computing power, and the availability of free satellite data. In particular, publicly available data from the European Space Agency’s Sentinel missions have been used in various remote sensing applications. However, there is a lack of studies that utilize these data to assess the performance of machine learning algorithms in complex boreal landscapes. In this article, I compare the classification performance of four non-parametric algorithms: support vector machines (SVM), random forests (RF), extreme gradient boosting (Xgboost), and deep learning (DL). The study area chosen is a complex mixed-use landscape in south-central Sweden with eight land-cover and land-use (LCLU) classes. The satellite imagery used for the classification were multi-temporal scenes from Sentinel-2 covering spring, summer, autumn and winter conditions. Using stratified random sampling, each LCLU class was allocated 1477 samples, which were divided into training (70%) and evaluation (30%) subsets. Accuracy was assessed through metrics derived from an error matrix, but primarily overall accuracy was used in allocating algorithm hierarchy. A two-proportion Z-test was used to compare the proportions of correctly classified pixels of the algorithms and a McNemar’s chi-square test was used to compare class-wise predictions. The results show that the highest overall accuracy was produced by support vector machines (0.758 ± 0.017), closely followed by extreme gradient boosting (0.751 ± 0.017), random forests (0.739 ± 0.018), and finally deep learning (0.733 ± 0.0023). The Z-test comparison of classifiers showed that a third of algorithm pairings were statistically different. On a class-wise basis, McNemar’s test results showed that 62% of class-wise predictions were significant from one another at the 5% level or less. Variable importance metrics show that nearly half of the top twenty Sentinel-2 bands belonged to the red edge (25%) and shortwave infrared (23%) portions of the electromagnetic spectrum, and were dominated by scenes from spring (38%) and summer (40%). The results are discussed within the scope of recent studies involving machine learning and Sentinel-2 data and key knowledge gaps identified. The article concludes with recommendations for future research.

Keywords:

• machine learning
• Sentinel-2
• land cover
• land use
• classification
• boreal

1. Introduction

In recent years, the data science and remote sensing communities have begun to align due to concurrent factors. First, popular competitions held by data science companies such as Kaggle (Google 2019) have demonstrated high classification accuracies using advanced machine learning algorithms. These competitions typically involve a sponsoring organization that posts its data on the Kaggle website for contestants to use as training input for their models and the winners are awarded cash prizes. Second, several new Earth-observing satellites, such as Sentinel-2 and Landsat-8, were recently launched and their data provided free to the public (Harris and Baumann 2015; Belward and Skøien 2015). Finally, consumer computing power has been dramatically increasing while its cost has been decreasing (Waldrop 2016). The combination of these factors has spurred the popularity of machine learning in the remote sensing and Earth-observation communities, particularly in the sub-field of land-cover and land-use (LCLU) classification.

One of the most widely used machine learning algorithms is random forests (RF) (Breiman 2001). The popularity of this algorithm is due to the fact it can be used for both classification and regression purposes, and thus can be used with categorical and continuous variables (Woznicki et al. 2019). Because of this flexibility, RF has been used in a wide range of Earth science applications including modeling forest cover (Betts et al. 2017), land-use (Araki, Shima, and Yamamoto 2018), land-cover (Nitze, Barrett, and Cawkwell 2015), and object-oriented mapping (Kavzoglu 2017). RF was compared with classification trees, which are also known as decision trees, by Rodriguez-Galiano et al. (2012) who found that RF produced a high accuracy of 92%, thereby outperforming classification trees. The higher accuracy of RF was attributed to its ensemble architecture in which several classification trees are trained on subsets of the training data.

Support vector machines (SVM) has been shown to outperform other classifiers due to its overall high capacity to generalize complex features (Shao and Lunetta 2012; Mountrakis, Jungho, and Ogole 2011). A land-cover classification study using Landsat-8 and involving six land-cover classes found that SVM was able to achieve a relatively high overall accuracy of 88% (Goodin, Anibas, and Bezymennyi 2015). Recently, Mansaray et al. (2019) analyzed the impact of training sample size on the overall accuracies of SVM and RF for mapping paddy rice in China in 2015 and 2016. They found that for 2015, SVM and RF achieved overall accuracies of 90.8% and 89.2%, respectively, using 10 satellite observations from Landsat-8 and Sentinel-1A. However, in 2016 SVM and RF achieved overall accuracies of 93.4% and 95.2%, respectively, using 14 satellite observations from Landsat-8, Sentinel-1A and Sentinel-2A.

Extreme gradient boosting (Xgboost) is a relatively new algorithm first described by Chen and Guestrin (2016). One of the earliest remote sensing applications of Xgboost was conducted by Georganos et al. (2018) using Bayesian parameter optimization on very-high-resolution WorldView-3 data. They found that Xgboost was able to outperform RF and SVM by 2–5% in larger sample sizes albeit with increased computational time. Man et al. (2018) compared five non-parametric classifiers using Landsat-8 data. They found that Xgboost slightly outperformed SVM by 0.3%. Xgboost was also found to slightly outperform RF by 0.2% in a recent six-class LCLU classification study using high-resolution data from RapidEye (Hirayama et al. 2019).

Deep learning (DL) is an abstract term that refers to a family of different algorithm architectures structured around neural networks. These architectures include multi-layer perceptrons, deep belief networks, stacked autoencoders, deep neural networks, and restricted Boltzmann machines, among others. Neural networks have been used in satellite image classification since at least the late 1980s and have been implemented in remote sensing software packages, usually with one or two hidden layers (Mas and Flores 2008). The number of hidden layers applied in satellite image classification remained low because data were expensive and computing power was inadequate. Considering the limitations at the time, some studies found no clear link between hidden layers and classification accuracy (Ardö, Pilesjö, and Skidmore 1997). This all changed at the turn of the 21^st century with the increased availability of Big Earth Observation Data and computing resources that merited the use of more (i.e. deeper) hidden layers and complex network architectures. Since 2015, DL has been used in a wide range of applications such as mapping land-cover (Li et al. 2016) and crops (Kussul et al. 2017; Zhong, Lina, and Zhou 2019), estimating crop yields (Kuwata and Shibasaki 2015), detecting oil palm trees (Li et al. 2017) and plant diseases (Mohanty, Hughes, and Marcel 2016) with accuracies often exceeding 90%.

1.1. Motivation and objectives

Some of the most important changes to the environment, such as urbanization, deforestation, and agricultural expansion, occur at the scale of landscapes and directly impact ecosystem processes (O’Neill et al. 1997; Belmaker et al. 2015). On the other hand, biotic interactions can physically alter landscapes and produce spatial patterns therein, a phenomenon termed ecosystem engineering (Hastings et al. 2007). Thus, it is vital to map the LCLU at the landscape scale in order to monitor and manage these changes. Classification using satellite data provides a crucial starting point for this endeavor.

Satellite data from the Sentinel-2 mission was chosen for this study for two main reasons (1) its relatively high 10 m spatial resolution, and (2) its radiometry includes three vegetation red edge bands. These two characteristics make the Sentinel-2 data appealing for LCLU mapping. However, there are presently no studies that utilize these data for assessing the performance of traditional and emergent machine learning algorithms for classification purposes and in complex boreal landscapes.

Thus, the objective of this study is to compare the classification performance of four popular machine learning algorithms over a boreal landscape in Sweden using Sentinel-2 data. The machine learning algorithms compared are two that have been widely used in the remote sensing community (SVM and RF) and two that are commonly used in the data science community (Xgboost and DL), and that are gaining popularity in remote sensing.

2. Study area

The study area (Figure 1) is a 10 km x 12 km mixed-use landscape located in the county of Uppsala in south-central Sweden. The mean annual temperature and precipitation are 5.6°C and 597 mm, respectively (Fick and Hijmans 2017). The elevation of the area is 40 meters above mean sea level (Tachikawa et al. 2011) and comprises stands of Norway spruce (Picea abies), Scots pine (Pinus sylvestris), and birch (Betula sp.), as well as extensive agriculture (Lundin et al. 1999). This area was specifically selected because of the weak annual amplitude of the vegetation greenness, which makes it difficult to capture seasonality (Jönsson et al. 2018) and thus poses a classification challenge, particularly for differentiation between vegetation types. The area also includes the Norunda research site that was established in 1994 for studies of greenhouse gas exchange, energy and water using the eddy covariance method. The research site is presently managed by Lund University and forms part of the Swedish contribution to the European research infrastructure Integrated Carbon Observation System.

Figure 1. a) Study area; b) location of the study area within the 33VXG Sentinel-2 tile; c) overview of the 33VXG tile’s coverage relative to Sweden and the Baltic region.

3. Data

3.1. Sentinel-2 data

The Sentinel-2 Multispectral Instrument (MSI) comprises two satellites that observe the Earth at 10 m, 20 m, and 60 m spatial resolutions (Drusch et al. 2012). The 10 m spatial resolution is the highest amongst freely available satellite products. Another unique aspect of the Sentinel-2 data is the presence of three red edge bands, which are able to capture the strong reflectance of vegetation in the near infrared portion of the electromagnetic spectrum (EMS).

The criteria for satellite imagery selection was that the scenes must contain little or no clouds and haze. The imagery must also be from different seasons to capture different plant phenological stages. As such, four Sentinel-2 scenes from 2017 and 2018 from tile 33VXG were included in the analysis. Three of the scenes were in 2017: May 04, July 06, and November 13, and one in 2018: January 27 (Figure 2). All the images were captured by Sentinel-2A except for the 2018 image, which was from Sentinel-2B. As shown in Table 1, the twin Sentinel satellites do not have identical band structures and the mean difference in the central wavelength of the ten bands used in this study is 3.3 ± 4.6 nm. The satellite imagery was downloaded from the Copernicus Open Access Hub (https://scihub.copernicus.eu/) on 29 May 2018. These were in Level-1C processing format, which means that they underwent geometric and radiometric correction but were not atmospherically corrected (Drusch et al. 2012). Atmospheric correction was performed using Sen2Cor (v2.5.5), which converts the top-of-atmosphere reflectance Level-1C data to a bottom-of-atmosphere (BOA) reflectance Level-2A product (Müller-Wilm 2018). BOA is also called “surface reflectance,” i.e. reflectance that would be measured at land surface, and is hereafter referred to as such. Ten bands that cover the red, blue, green, red edge, near- and short-wave infrared portions of the EMS were selected for inclusion into the classification procedure (Table 1). All the bands at 20 m spatial resolution were resampled to 10 m using bilinear interpolation to facilitate integration and consistency.

Table 1. Descriptions of the 10 Sentinel-2 bands used in this study. SR = Spatial Resolution, CW = Central Wavelength.

			S2A		S2B
Band	Spectral Region	SR(m)	CW(nm)	Bandwidth (nm)	CW(nm)	Bandwidth (nm)	CW Difference
2	Blue	10	496.6	98	492.1	98	4.5
3	Green	10	560.0	45	559	46	1
4	Red	10	664.5	38	665	39	0.5
8	Near Infrared	10	835.1	145	833	133	2.1
5	Red Edge	20*	703.9	19	703.8	20	0.1
6	Red Edge	20*	740.2	18	739.1	18	1.1
7	Red Edge	20*	782.5	28	779.7	28	2.8
8A	Near Infrared	20*	864.8	33	864	32	0.8
11	Shortwave Infrared	20*	1613.7	143	1610.4	141	3.3
12	Shortwave Infrared	20*	2202.4	242	2185.7	238	16.7

*Resampled to 10 m.

Figure 2. Natural color composites of the multi-temporal imagery used in this study. a) 4 May 2017 (Sentinel-2A); b) 6 July 2017 (Sentinel-2A); c) 13 November 2017 (Sentinel-2A); d) 27 January 2018 (Sentinel-2B).

Three spectral indices were derived from the resultant bands of each Sentinel-2 scene. The normalized difference vegetation index (NDVI) quantifies the ratio between energy absorbed by the vegetation canopy in the red portion of EMS and the energy reflected in the near infrared (NIR) (Rouse et al. 1973). The modified normalized difference water index (MNDWI) is a ratio index that maximizes the reflectance of water by using the green portion of the EMS while using the shortwave infrared (SWIR) portion to suppress the influence from artificial surfaces (Xu 2006). The normalized difference built-up index (NDBI) uses the SWIR band to identify artificial surfaces because of their strong reflectivity in that portion of the EMS while using the NIR to suppress the influence of vegetated surfaces (Zha, Gao, and Ni 2003). NDVI, MNDWI, and NDBI were included in the classification process in order help capture vegetation, water, and artificial surfaces, respectively. This resulted in each Sentinel-2 scene having 13 layers (10 bands + 3 indices) for a total of 52 layers.

4. Methods

4.1. Study design and sample selection

Training data are crucial components in supervised learning and most machine learning algorithms require a large number of training data samples. However, the delineation and acquisition of reference data from satellite imagery can be a daunting task (Chi, Feng, and Bruzzone 2008). In-situ collection of representative training data is difficult over large, remote or distant areas (Inglada et al. 2017) and the use of timely high-resolution satellite or aerial imagery for this purpose is not always feasible and often cost-prohibitive. The use of old maps as ancillary data for classifying past land-cover is not an uncommon practice in the remote sensing community. For example, Tran, Tran, and Kervyn (2015) used thematic maps from the 1970s as training input to classify Landsat-1 data from 1973 over a study area in Vietnam. A novel approach being applied in a growing number of recent studies is that training data are acquired from extant high-quality land-cover maps (Wessels et al. 2016; Zhang and Roy 2017; Hermosilla et al. 2018). The use of an existing maps as training data for LCLU classification can introduce errors inherent in the previous classification. However, land-cover maps with reasonably high overall accuracy can produce a large number of training samples with increased efficiency allowing for a wide range of feature representation (Hermosilla et al. 2018). Further, it is important to randomly check the accuracy of the training samples collected from these land-cover maps in order to minimize misclassification, for example, to ensure that a sample labeled as “forest” actually falls within a forest.

The study design (Figure 3) involves initial training sample selection with orthophotos then increasing the sample size using a high-quality LCLU map of the study area. A total of 100 training samples were initially collected for each land-cover class in the study area using random sampling from 25 cm orthophotos of the study area captured on the 2^nd and 3^rd of July 2015 by the Swedish Cadastral and Land Registration Authority (https://zeus.slu.se). These samples were then compared to RGB and false color composites of the multi-temporal Sentinel-2 scenes to discern any visible change that took place between 2015 and 2017/18 such as clearcutting or infrastructure development. The sample labels were modified in the event that a change was detected. In the next step, each of the samples was compared to the 2018 National Land Cover Database of Sweden (Nationella marktäckedata, NMD), produced by the Swedish Environmental Protection Agency with a minimum mapping unit of 0.01 hectares. Again, sample labels were modified where appropriate as described in the next paragraph. The NMD dataset was the result of an integrated processing chain that involves several ancillary datasets from Swedish Cadastral and Land Registration Authority (Lantmäteriet), Statistics Sweden (Statistiska centralbyrån), road and railway network data, 25 cm orthophotos, a 2 m digital elevation model, and airborne lidar. The NMD is part of Sweden’s national geo-strategy for 2016–2020 and is available from https://www.naturvardsverket.se/.

Figure 3. Flowchart of the methods.

Fifteen of the 25 NMD classes occurred in the study area and were aggregated to eight classes (Table 2). These are: (1) Deciduous forest: comprising birch, aspen, alder, beech, oak, elm and ash; (2) Coniferous forest: comprising Scots pine and Norway spruce; (3) Water: lakes, rivers and canals; (4) Artificial: urban areas, construction sites, and roads; (5) Wetlands: saturated land including marshes and bogs; (6) Agriculture; (7) Clear cuts within and outside wetlands; and (8) Open land with and without vegetation. The remaining seven classes were discarded due to negligible presence in the study area. A quantitative sample-by-sample inspection resulted in a low (1.75%) mismatch between the NMD dataset and the 800 samples (100 samples per class) collected using the orthophotos, and a qualitative visual inspection showed a good match between the NMD classes and the orthophotos (Figure 4). In this way, the reference data can reliably replicate on-ground conditions and meet the good practice guidance for accuracy assessment (Olofsson et al. 2014).

Table 2. Total number of pixels in each sample and the number of pixels selected for the training and evaluation. The sum of training and evaluation pixels represents 1% of all pixels in the study area. Also shown are the original Swedish National Land Cover Database (NMD) classes that were aggregated to form each of the eight land cover types in this study. Forest classes are all outside of wetlands.

LCLU Type	Class ID	Original NMD Class Description	Pixels (Total)	Samples (Training/Evaluation)
Deciduous	Class 1	1.1.5 Birch, Aspen, and Alder forest	180,827	8271/3544
		1.1.6 Beech, Oak, Elm, and Ash forest
Coniferous	Class 2	1.1.1 Scots pine forest	411,356	8271/3544
		1.1.2 Norway spruce forest
		1.1.3 Mixed coniferous forest
Water	Class 3	6.1 Water	45,222	8271/3544
Artificial	Class 4	5.1 Developed land, building	52,839	8271/3544
		5.2 Developed land, no building/road
		5.3 Developed land, road
Wetland	Class 5	2 Wetland	44,835	8271/3544
Agriculture	Class 6	3 Agriculture	261,391	8271/3544
Clear Cut	Class 7	1.1.8 Clear cut outside wetlands	125,476	8271/3544
		1.2.8 Clear cut within wetlands
Open Land	Class 8	4.1 Open land without vegetation	56,497	8271/3544
		4.2 Open land with vegetation

Next, the number of samples was increased using the NMD dataset due to the data-driven framework of machine learning models (Brink, Richards, and Fetherolf 2017). Stratified random sampling was applied in order to produce 1477 samples per class. This number was selected because it represents 1% of all the pixels in the study area when all the class samples are summed (1477 x 8 = 11,816). This amount of samples was chosen because it is large enough to be adequately dispersed across the study area (Figure 5) and does not exhaust available computing power. The 11,816 samples were then split into two portions: a training dataset comprising 70% (8271) of the samples and an evaluation dataset comprising the remaining 30% (3544). Each LCLU class was assigned the same number of training (1034) and evaluation (443) samples (Table 2).

Figure 4. Comparison between the 25 cm orthophotos from 2015 (a-labeled panels) and the 2018 Swedish National Land Cover Database (NMD) map (b-labeled panels) for different parts of the study area.

Figure 5. Overview of the training (yellow dots) and evaluation (red dots) samples in selected parts of the study area: (a) The full study area; (b) red box: within the vicinity of the Norunda research station; (c) turquoise box: agricultural fields of Örbyhus in the northeast of the study area; (d) yellow box: the town of Björklinge in the southeast..

4.2. Machine learning algorithms

This section provides descriptions of the algorithms used in the classification. Readers seeking deeper understanding of the theoretical background of a particular algorithm should consult the reference provided at the end of each description.

4.2.1. Support vector machines

SVM was first described in Cortes and Vapnik (1995) based on the work of Vapnik (1982) and is a supervised learning technique commonly used in a range of remote sensing applications. The SVM algorithm finds the optimum minimization, i.e. decision boundary, of ambiguous classifier outputs in a problem space. This decision boundary is referred to as the hyperplane and it distinguishes the classification problem into a predefined set of classes that are consistent with the training data. The algorithm undergoes an iterative process of finding the optimum hyperplane boundary in an n-dimensional classification space to distinguish patterns in the training data then apply the same configuration to a separate evaluation dataset. The dimensions in this context are the number of spectral bands and the vectors are the individual pixels in a multiband composite (Mountrakis, Jungho, and Ogole 2011). There are different kernels through which the hyperplane boundary can be defined. Here, a radial basis function kernel was used because some of the Sentinel-2 bands are not linearly separable. A detailed mathematical description of this algorithm can be found in Cortes and Vapnik (1995).

4.2.2. Random forests

RF is an ensemble learning algorithm based on the idea that a combination of bootstrap aggregated classifiers perform better than a single classifier (Breiman 2001). The bootstrap component means that each individual tree is parameterized using a randomly sampled set of observations with replacement from the training data (Hastie, Tibshirani, and Friedman 2009). This helps to de-correlate the trees thereby reducing multicollinearity. The proportion of observations that are not used for this purpose are included in the evaluation and are referred to “out-of-bag” samples. Several of these decision tree models are created on different groupings of the input variables and the resultant output is the unweighted majority vote of each class that is averaged across all trees. A detailed mathematical description of RF is provided in Breiman (2001).

4.2.3. Extreme gradient boosting

The classical gradient boosting machine (GBM) builds an additive model of shallow decision trees that are weak learners and then generalizes them by optimizing an arbitrarily defined loss function to make stronger predictions (Friedman 2001). Xgboost is a relatively new implementation of the GBM that simultaneously optimizes the loss function while building the additive model (Chen and Guestrin 2016). The novelty of Xgboost lies in the fact that it comprises an objective function, which combines the loss function and a regularization term that controls model complexity. This enables parallel calculations and the maintenance of optimal computational speed. The softmax multiclass classification objective function was used in this study. Softmax is a function that normalizes each class into a probability distribution with an interval of (0, 1) that sums to 1. A detailed mathematical description of Xgboost is provided in Chen and Guestrin (2016).

4.2.4. Deep learning

The type of DL architecture implemented in this study is a multilayered feed-forward deep neural network (DNN) with error back-propagation. DNN links several functions joined together in hierarchically structured neural networks that are typically deeper than three hidden layers. The information flows from the input data through the activation function, error is calculated and propagated back to the earlier layers, and finally fed to the output at the conclusion of the predefined iteration (Candel and Erin 2018; Goodfellow, Bengio, and Courville 2016). In this way, every sequential hidden layer merges values in the previous layer and subsequently learns to form more abstract representations until the output layer. The algorithm exploits the activation function by learning the values of parameters that result in the best approximation (Beysolow 2017). These are then taken by an output function that computes class probabilities. The hyperbolic tangent (tanh) activation function with a softmax output classification function were used in this study. Tanh is a rescaled logistic sigmoid function, i.e. $f(x)={\left(1+e^{-x}\right)}^{-1}$, that provides zero-centered outputs and allows model parameters to be more frequently updated in feed-forward neural networks. A detailed description of the DL method applied in this study is provided in Goodfellow, Bengio, and Courville (2016) and in Candel and Erin (2018).

4.3. Model training and hyper-parameter optimization

Hyper-parameter optimization (also called tuning) was part of the model training process whereby optimal hyper-parameters were selected for the algorithms. Models were trained by optimizing hyper-parameters using a repeated k-fold cross-validation. In the case of SVM, RF and Xgboost, optimization was performed using randomized sampling of all hyper-parameter combinations up to a specified number of iterations (Kuhn and Johnson 2013b; Bergstra and Bengio 2012). This is called the tune length and it defined the total number permutations that were evaluated. The training dataset was randomly split into k sets of equal size, of which one set was retained for validation and the remaining k–1 samples were used for training. The entire procedure was then repeated using a predefined set of repeats in order to reduce model variance (Kim 2009; Kuhn and Johnson 2013a). A 10-fold cross-validation with 5 repeats was chosen as a tradeoff between lowering the variance, ensuring a robust model, and a reasonable computational time.

In case of DL, optimization was performed using a random grid search. Here, the algorithm trained a model for all possible combinations of the hyper-parameters in the grid and selected the best one. Once the tune length was reached, the performance metrics from the cross-validation were used to select the best model parameters (Table 3). SVM, RF, and Xgboost were run using the caret package (Kuhn 2008) in the R statistical software environment version 3.4.2 (R Core Team 2017). DL was performed using the H2O package, also in R (Cook 2016). All processing was performed on an 8-core 3.60 GHz Xeon server with 48 GB of RAM running Windows 10 64-bit.

Table 3. Results of the hyper-parameter optimization process showing the final values used in each model. The tune length was set to 1000 iterations. SVM = Support Vector Machines, RF = Random Forest, XGB = Extreme Gradient Boosting, DL = Deep Learning. Algorithms were run in R using the caret package (Kuhn 2008) and the H2O package (Cook 2016).

Model	Hyper-parameter Value	Definition
SVM	sigma = 0.00425, C = 64, Kernel = “rbf”	rbf = radial basis function, a multidimensional Gaussian distribution function describing the distance between an input vector and a pre-defined center vector. sigma = weight of the RBF kernel. C = cost of misclassification.
RF	mtry = 7, ntree = 1000	mtry = number of variables randomly sampled as candidates at each split. ntree = number of trees.
XGB	nrounds = 500 max_depth = 10 eta = 0.05 gamma = 0.3 nodesize = 1	nrounds = maximum number of iterations. max_depth = maximum depth of a tree. eta = learning rate by which to shrink the feature weights. gamma = minimum relative improvement in squared error reduction in order for a split to happen. nodesize = minimum number of rows to assign to the terminal nodes.
DL	activation = “tanh,” hidden_layers = 6, neurons_per_layer = 200, epochs = 500	activation = activation function. hidden_layers = number of hidden layers. neurons_per_layer = size of each hidden layer. epochs = number of times to iterate.

4.4. Accuracy assessment and area estimation

The accuracy of each algorithm was assessed using a number of metrics derived from an error matrix. These include Overall Accuracy (OA), Producer’s Accuracy (PA), and User’s Accuracy (UA). Error matrices were enhanced by providing unbiased estimations of the proportional area of cells within the matrix as per Olofsson et al. (2014). The area proportions of the mapped classes were included in the output because each class required its own estimation weight. These proportions were necessary to estimate OA and PA in order to include differences in sampling between classes. Conversely, UA was computed from within a given LCLU class and was thus retrieved directly from the error matrix (Olofsson et al. 2013).

McNemar’s chi-squared (χ²) test (McNemar 1947) was used to statistically compare error matrices by testing for the marginal homogeneity between two classifiers. Marginal homogeneity refers to the equality (i.e. lack of statistically significant difference) in the overall distributions of row or column variables predicted by one algorithm compared to another. It is a simple yet powerful method to compare class-wise predictions between algorithms. The test is parametric, has a low type I error and consists of a straightforward formulation (Dietterich 1998; de Leeuw et al. 2006). Additionally, a two-proportion Z-test (Lachin 1981) was used to compare the proportions of correctly classified pixels (PCCP) from two algorithms at a time. This test produced a two-tailed probability value that tests the null hypothesis of no difference between PCCP of each algorithm pair. The square of the Z-statistic produced by the test followed a χ² distribution with one degree of freedom (Wallis 2013). For both methods, a χ² value of greater than 3.84 indicated statistically significant difference at the 5% level.

The conventional method of computing areas of mapped classes involves multiplying the area of a pixel by the total number of pixels in a class. This method does not account for classification errors (Czaplewski 1992) and introduces bias into the resultant LCLU map. Therefore, in order to provide unbiased area estimates of the output classes, error-adjusted stratified estimation of the error matrix was conducted following Olofsson et al. (2013). Each area estimate was accompanied by a confidence interval (95%) that quantified uncertainty. For a thorough explanation of the method used to acquire these unbiased accuracies and quantify uncertainty in area estimates, readers are strongly encouraged to consult the good practices guide by Olofsson et al. (2014) and the recommendations described in Olofsson et al. (2013).

5. Results and discussion

5.1. Classifier comparison

The OA of the four algorithms were relatively close to one another as shown in Table 4. These accuracies were derived from the error matrix shown in Figure 6. The highest OA was produced by SVM (0.758 ± 0.017), closely followed by Xgboost (0.751 ± 0.017) and RF (0.739 ± 0.018), and finally DL (0.733 ± 0.0023). All four algorithms produced similar maps that were visually appealing (i.e. coherent classes and minimal speckling) and represented the area fairly well (Figure 7). The adjusted areas of the classes within a 95% confidence interval were similar for SVM, Xgboost and RF. However, relative to the other algorithms, DL overestimated coniferous forests and open land, and underestimated the remaining classes except water (Figure 8). All eight LCLU classes had the same number of training and evaluation samples so that the OA is not biased towards classes with more training samples (He and Garcia 2008). One drawback of balanced training samples in conventional accuracy assessment is that the proportions of land-cover classes are not taken into consideration. However, area proportions were taken into account in this study using by producing unbiased OA estimates following Olofsson et al. (2014), which could have led to lower OA. It is possible to obtain higher OA values using an imbalanced class distribution where the accuracy of certain well-represented classes increases OA. But, this focus on OA ignores the performance of individual classes and is particularly disadvantageous to those that are underrepresented in the study area (Maxwell, Warner, and Fang 2018). Therefore, the choice of sampling design should be dependent upon whether the objective of the study is to obtain the highest possible OA regardless of class distribution or equally represent all LCLU classes.

Table 4. Area-adjusted performance metrics for each machine learning algorithm and LCLU class derived from the error matrices in Figure 6. The proportional weights of each class that were used to create the unbiased error estimate for each class. SVM = Support Vector Machines, RF = Random Forest, XGB = Extreme Gradient Boosting, DL = Deep Learning, OA = Overall Accuracy.

		Producer’s Accuracy	User’s Accuracy	Class Weight	Area in km^2	Number of Pixels
SVM OA = 0.758 ± 0.017	Deciduous	0.63	0.69	0.1412	17	166,886
	Coniferous	0.94	0.74	0.3216	38	380,212
	Water	0.93	0.99	0.0381	5	45,047
	Artificial	0.45	0.72	0.0589	7	69,645
	Wetland	0.69	0.90	0.0466	6	55,134
	Agriculture	0.93	0.85	0.1938	23	229,077
	Clear cut	0.73	0.71	0.1280	15	151,311
	Open land	0.49	0.64	0.0718	8	84,900
RF OA = 0.739 ± 0.018	Deciduous	0.63	0.64	0.1461	17	172,733
	Coniferous	0.93	0.74	0.3163	37	373,906
	Water	0.92	0.98	0.0381	5	45,090
	Artificial	0.40	0.72	0.0535	6	63,212
	Wetland	0.67	0.85	0.0486	6	57,419
	Agriculture	0.93	0.86	0.1954	23	230,974
	Clear cut	0.67	0.65	0.1256	15	148,453
	Open land	0.47	0.59	0.0765	9	90,425
XGB OA = 0.751 ± 0.017	Deciduous	0.64	0.68	0.1427	17	168,754
	Coniferous	0.94	0.75	0.3175	38	375,365
	Water	0.93	0.98	0.0388	5	45,872
	Artificial	0.45	0.72	0.0586	7	69,301
	Wetland	0.65	0.89	0.0470	6	55,573
	Agriculture	0.93	0.85	0.1957	23	231,346
	Clear cut	0.70	0.67	0.1247	15	147,476
	Open land	0.46	0.59	0.0749	9	88,525
DL OA = 0.733 ± 0.0023	Deciduous	0.58	0.63	0.1334	16	157,741
	Coniferous	0.95	0.75	0.3548	42	419,436
	Water	0.93	0.98	0.0384	5	45,381
	Artificial	0.36	0.71	0.0424	5	50,078
	Wetland	0.64	0.88	0.0451	5	53,309
	Agriculture	0.90	0.78	0.1884	22	222,755
	Clear cut	0.70	0.66	0.1173	14	138,674
	Open land	0.44	0.64	0.0802	9	94,838

Figure 6. Error matrices showing correct and incorrect cross-tabulations of the evaluation samples by each machine learning algorithm. The LCLU classes are numbered 1 through 8, where 1 = Deciduous, 2 = Coniferous, 3 = Water, 4 = Artificial, 5 = Wetland, 6 = Agriculture, 7 = Clear cut, 8 = Open land. The error matrices form the basis for the calculation of the performance metrics shown in Table 4.

Figure 7. Final classified maps of the study area for (a) Support Vector Machines; (b) Random Forests; (c) Extreme Gradient Boosting; (d) Deep Learning.

Figure 8. Adjusted area estimates of each LCLU class resulting from image classification using the four algorithms tested. The error bars denote the 95% confidence interval. SVM = Support Vector Machines, RF = Random Forest, XGB = Extreme Gradient Boosting, DL = Deep Learning.

The lowest PA across all four algorithms was for open land and artificial surfaces (Table 4), which highlights the difficulty distinguishing these classes from the rest. Open land was composed primarily of grasses and herbs, and was visually similar to clear cuts, particularly those with regrown shrubs and herbaceous cover. Thus, it remained difficult to classify due to the similar spectral properties of the vegetation found in it (Buck et al. 2015). On the other hand, artificial surfaces in the study area consisted of suburban features such as small towns and villages interspersed with vegetation such as trees and grass. The inclusion of land surface temperature (Abdi 2019) could have enhanced the detection of man-made features but there is presently no straightforward way to quantify emissivity because the MSI lacks a thermal band.

SVM produced the highest accuracy for open land (PA = 0.49) among the four classifiers, and tied with Xgboost for artificial surfaces (PA = 0.45). The ability of SVM to capture these classes was probably because it is based on a relatively small number of complex decision boundaries (Cortes and Vapnik 1995). In cases such as open land and clear cuts, where the data were linearly non-separable, the feature vectors were projected with a nonlinear vector mapping function to a higher dimension feature space. This facilitates the creation of a decision boundary that seems nonlinear in the original feature space. The added computation intensity of this projection was offset, to a certain extent, by kernel functions that enable a simplified representation of the data (Mountrakis, Jungho, and Ogole 2011). Indeed, Khatami, Mountrakis, and Stehman (2016) found that SVM was the most efficient algorithm for most applications and outperformed several classifier families, including RF, neural networks, and decision trees in direct comparison.

A somewhat surprising result was that DL produced the lowest accuracy for open land and for artificial surfaces (PA = 0.44 and 0.36, respectively), and was overall the poorest performing algorithm. Use of the tanh activation function was the probable reason for the underperformance of DL in this study. The tanh function easily saturates (i.e. the so-called “vanishing gradient problem”) and slows the training procedure when applied to DNNs (Rakitianskaia and Engelbrecht 2015). In the case of high-dimensional data such as multi-temporal satellite imagery, this saturation was possibly dependent on feature complexity. The tanh activation function has a conventional sigmoid curve that is centered at zero and is restricted to the range of (−1, 1). It returns a near-zero slope when the input values are large (Shi et al. 2018). This suggests that it was incapable of the separating the complex multiclass features in the study area. The DL classes with the highest PA were those that can be relatively easily distinguished from multi-temporal data such as coniferous forests, water or agriculture (Table 4). The lowest performing DL classes are spectrally complex and thus had considerable variation in the relative intensity of pixel reflectance and the tanh activation function performed poorly in distinguishing between them. In retrospect, the rectified linear units (ReLU) activation function may have provided a better output due to its ability to hold the incremental gradient descent, i.e.$f(x)=max\left[0,x\right]$, because it is a non-saturating nonlinearity function (Nair and Hinton 2010).

McNemar’s test results (Table 5) showed that 42% of the class-wise predictions between paired algorithms exhibited significance at the 1% level (P ≤ 0.01); 20% at the 5% level (P ≤ 0.05); and 11% at the 10% level (P ≤ 0.10). Additionally, 27% of class-wise algorithm pairings were statistically non-significant. Most of the non-significant pairings were for wetland, except between Xgboost and SVM (P ≤ 0.05) and between Xgboost and DL (P ≤ 0.10). This means that apart from these two pairings, accuracies for wetland between the other classifiers were statistically insignificant, i.e. there is no real difference between them. Both coniferous forests and clear cuts were consistent at the 1% level across all classifier pairings, suggesting that the accuracies achieved for these classes were highly significant.

Table 5. McNemar’s chi-squared (χ²) test with associated probability value (P). SVM = Support Vector Machines, RF = Random Forest, XGB = Extreme Gradient Boosting, DL = Deep Learning. *** = P ≤ 0.01, ** = P ≤ 0.05, * = P ≤ 0.1, NS = Not Significant.

	McNemar’s Test	Deciduous	Coniferous	Water	Artificial	Wetland	Agriculture	Clear Cut	Open Land
SVM v. RF	χ^2	11.589	42.654	3.184	6.298	0.985	0.041	27.876	7.557
	P	***	***	*	**	NS	NS	***	***
XGB v. SVM	χ^2	18.178	45.797	2.132	3.838	4.621	0.071	31.045	4.951
	P	***	***	NS	**	**	NS	***	**
DL v. SVM	χ^2	5.320	36.445	3.184	1.098	1.071	3.923	52.115	59.602
	P	**	***	*	NS	NS	**	***	***
RF v. XGB	χ^2	8.983	41.150	1.841	12.445	2.756	0.018	28.019	3.71
	P	**	***	NS	***	NS	NS	***	**
DL v. RF	χ^2	1.227	32.345	2.750	6.892	0.386	4.498	47.036	52.998
	P	NS	***	NS	***	NS	**	***	***
XGB v. DL	χ^2	3.432	35.003	1.841	4.331	2.984	4.188	51.336	46.066
	P	*	***	NS	*	*	**	***	***

The Z-test results (Table 6) showed that the PCCP of two algorithm pairings were statistically significant thereby rejecting the null hypothesis. These are SVM and RF (χ² = 4.832), and SVM and DL (χ² = 4.701). This indicates that the PCCP produced by SVM is markedly different from those produced by RF and DL despite the fact their OA are fairly close to one another (Table 4). For the remaining four algorithm pairings, the null hypothesis was not rejected, meaning their accuracies are similar. Although the PCCP produced by SVM and Xgboost are not statistically different from one another (χ² = 1.369), the latter failed to reject the null hypothesis when compared to other algorithms (Table 6).

Table 6. A Two Proportion Z-test to compare the proportions of correctly classified pixels with its associated probability value (P). SVM = Support Vector Machines, RF = Random Forest, XGB = Extreme Gradient Boosting, DL = Deep Learning. Algorithm pairs that exhibit a statistically significant (P ≤ 0.05) difference in the proportion of correctly classified pixels are in bold.

Algorithm Pair	Two Proportion Z-test
SVM v. RF	χ^2 = 4.832
	P = 0.027
XGB v. SVM	χ^2 = 1.369
	P = 0.241
SVM v. DL	χ^2 = 4.701
	P = 0.030
RF v. XGB	χ^2 = 0.996
	P = 0.318
DL v. RF	χ^2 = 0
	P = 1
XGB v. DL	χ^2 = 0.937
	P = 0.333

A possible limitation of the classification process was the choice of user-defined parameters for each algorithm (Table 3). The hyper-parameter optimization applied here was based on iterative tuning and a random grid search without empirically examining their optimum values. This was done in order to compare the algorithms in a fair manner and not introduce a priori knowledge. Thus, the accuracies reported here may not represent the maximum possible that could be attained. That said, the removal of open land, which was one of the confusion classes, increased the overall accuracy for all algorithms by an average of 15% with Xgboost having 0.7% higher accuracy than SVM (not shown). A complex landscape invariably includes classes that are spectrally similar to one another, such as open land and clear cuts in this study, but that have different land-use categorizations. Therefore, an algorithm’s ability to generalize training samples within such classes will ultimately be reduced. This leads to a non-separation of these classes in a feature space causing overlap between class distributions and a lower overall classification accuracy (Jansen and Gregorio 2002; Smits, Dellepiane, and Schowengerdt 1999).

5.2. Variable importance metrics

All four algorithms provide variable importance metrics as part of their output (Supplementary Figure 1–4). The twenty highest scoring layers across the three algorithms displayed two generalizable patterns of importance. First, nearly half belonged to the red edge (25%) and shortwave infrared (SWIR) (23%) bands. Second, they were dominated by scenes from May (38%) and July (40%). None of the spectral indices were ranked highly in this study. In the case of NDVI, it was probably due to the red edge bands providing most, if not all, of the information required to account for the vegetation signal. Additionally, NDVI is known to saturate at high values (Box, Holben, and Kalb 1989), which further degrades its ability to capture high-biomass classes such as forests. The absence of the other two indices could be explained by inability to separate urban areas and barren land in the case of NDBI (Zha, Gao, and Ni 2003), and between water and shadowed surfaces in the case of MNDWI (Feyisa et al. 2014).

A distinguishing feature of the Sentinel-2 satellites is the presence of red edge bands that provide an added value previously only available in commercial satellites such as RapidEye. The red edge region is so-called because leaves reflect strongly between 680 nm and 780 nm and is dependent on chlorophyll concentration (Horler, Dockray, and Barber 1983; Gates et al. 1965). The advantage of the red edge bands in LCLU classification is evident when Sentinel-2 data is compared with other satellites such as Landsat-8 that do not capture this portion of the EMS. For example, red edge Band 5 (705 nm) was found to be important for mapping crop types, particularly cereals and legumes (Forkuor et al. 2018), red edge Band 6 (740 nm) was useful for differentiating crop types from grassland at the sub-pixel level (Radoux et al. 2016), and red edge band 7 (783 nm) was used for classifying edge pixels in boreal forests comprising Scots pine, Norway spruce and birch (Zerega 2018).

The high importance of the SWIR bands could be due to the dominance of forest classes in the study area (Eklundh, Harrie, and Kuusk 2001). Sentinel-2’s SWIR bands are centered at 1610 nm (Band 11) and 2190 nm (Band 12) and were thus able to detect variability in water content between different tree species (Lukeš et al. 2013). The combination of red edge bands with SWIR has been shown to enhance the modeling of boreal forest canopy cover and leaf area index, producing root-mean-square error values 1.6–7.2% lower than Landsat-8 (Korhonen, Hadi, and Rautiainen 2017). Sentinel-2 Band 11 was also found to correlate well (r = −0.68, p = 0.01) with Mediterranean forest growing stock volume (Chrysafis et al. 2017), although in a Swedish forest Zerega (2018) found a lower correlation for Band 11 (r = −0.26) and higher correlations (−0.51 and −0.52) for red edge bands 6 and 7, respectively. Recent work by Persson, Lindberg, and Reese (2018) to classify Norway spruce, Scots pine, hybrid larch, birch and pedunculate oak found that red edge bands captured in May and July and SWIR bands captured in April and October had the highest separation of the tree species.

A possible limitation of the variable selection process was the fact that different algorithms calculate variable importance differently. For example, the DL implementation in H2O used a method by Gedeon (1997) that accounts for each input variable’s connecting weights to the first two input layers as a basis for computing its importance. Both RF and Xgboost computed variable importance in similar ways that were based on the mean decrease of a predefined loss function such as mean squared error (Breiman 2001; Chen and Guestrin 2016). Despite these differences, the results (Supplementary Figure 1–4) illustrate the consistency with which red edge and SWIR bands from May and July were selected by the algorithms.

5.3. Mono-temporal vs. multi-temporal Sentinel-2 data

The rationale for incorporating data from four seasons was to help distinguish vegetation classes that had different phenological cycles (e.g. coniferous, deciduous, crops). However, the importance of the scenes from May and July was clear in this study, which was probably due to the dominance of forest area (Figure 8), and the increased leaf expansion in the boreal zone that usually takes place in late spring – early summer (Tang et al. 2016). Increase in the light-saturated mean daily rate of photosynthesis also takes place around that time as the seasonal light intensity increases with the approaching boreal summer (Sukenik, Bennett, and Falkowski 1987; Letts et al. 2008). Canopy maturation takes place around July, which is also an important time for the maturation of certain crops, such as wheat, that require relatively dry sunny conditions (Brown 2013; Tang et al. 2016).

The use of multi-temporal data may have inadvertently increased the level of noise in the classification process as increased data dimensionality leads to higher redundancy (Nitze, Barrett, and Cawkwell 2015). Indeed, bootstrapped accuracy provided by RF showed a steep decline with the addition of new variables after reaching peak accuracy with just six bands. This may suggest that mono-temporal (i.e. a single scene) Sentinel-2 surface reflectance data, captured within an optimal date (Vuolo et al. 2018), may offer an alternative to multi-temporal time series. This could be due to the relatively high radiometric resolution of Sentinel-2 data, particularly the presence of three red edge and two near infrared bands that enable the satellites to sense important biophysical parameters. Furthermore, differences in radiometry between Sentinel-2A and Sentinel-2B were unlikely to have impacted the results as the mean difference in the central wavelength of the bands used was 3.3 ± 4.6 nm and there was considerable overlap between the bandwidths of the two satellites (Table 1).

Early work on Sentinel-2 data by Immitzer, Vuolo, and Atzberger (2016) found that mono-temporal Sentinel-2 red edge and SWIR were important for mapping both tree species and crop types. However, the timing of the imagery used (13 August for trees and 30 August for crops) caused confusion within these classes. The dates of these acquisitions were not optimal because most crops were approaching senescence while images acquired at the end of spring would have better differentiated tree species (Immitzer, Vuolo, and Atzberger 2016; Persson, Lindberg, and Reese 2018). The authors thus concluded that they expect multi-seasonal information contained within time-series of Sentinel-2 data to provide improved classification results. Persson, Lindberg, and Reese (2018) used Sentinel-2 images from four months (April, May, July and October) in 2017 to classify tree species in Sweden. They found that although the highest OA of 88.2% was obtained when all bands from all four dates were included, the highest accuracy provided by a single image (80.5%) was the one from May.

It is not implausible that a single, optimally-timed, atmospherically-corrected Sentinel-2 scene could provide equal or better results than a multi-temporal stack for detecting spectral differences in vegetation due to the accumulation of noise that is inherent in such data (Feilhauer et al. 2013). Of course, this is highly dependent on the aims that a study hopes to achieve. For example, Prishchepov et al. (2012) found that in order to map land abandonment with classification accuracies greater than 80%, it was necessary to use at least three multi-temporal images within a single year before and after abandonment. Phenological information contained within multi-temporal data does allow for better discrimination between spectrally similar and phenologically distinct vegetation. However, it should be used with caution in managed land due to anthropogenic changes that could be confused for phenological change (Ghioca-Robrecht, Johnston, and Tulbure 2008).

6. Conclusions and recommendations

The Sentinel-2 Multispectral Instrument is unique among presently operating Earth-observation satellites due to its three red edge bands that can capture plant chlorophyll content and its medium-high spatial resolution of 10 m. This study represents a first assessment of traditional and emergent machine learning algorithms to classify land-cover and land-use using multi-temporal Sentinel-2 data over a complex boreal landscape. Two of the machine learning algorithms, support vector machines and random forests, are widely used in the remote sensing community. The other two, extreme gradient boosting and deep learning, are commonly used in the data science community but are gaining popularity in remote sensing.

The four tested algorithms produced similar overall accuracies ranging between 0.733 to 0.758. The Z-test comparison of classifier accuracies showed that a third of algorithm pairings were statistically distinct from one another, while McNemar’s test results showed that 62% of class-wise predictions were significant. The highest classification accuracy produced by support vector machines is due to the algorithm’s relatively small number of complex decision boundaries. The lowest performance produced by deep learning is probably due to a large number classes and the saturation of the hyperbolic tangent activation function used in the study. Finally, the variable importance metrics show that nearly half of the top twenty bands belonged to the red edge and shortwave infrared bands and were dominated by scenes from May and July.

Studies that compare machine learning algorithms should do so in a consistent manner in order to not introduce bias. The algorithms being compared should undergo equally robust hyper-parameter selection. Bias can be introduced in an experiment when some algorithms are executed with hyper-parameter values from a priori knowledge while those of other algorithms are calibrated using default or random values. Random iteration across a defined number of parameter combinations, as done in this study, can be used to eliminate a priori knowledge but at the cost of the algorithms not reaching optimum accuracies. Modifications to this approach include (1) significantly increasing the number of random iterations, which will invariably increase computing requirement, or (2) individually assessing optimum hyper-parameter thresholds for each algorithm.

Spectral vegetation indices are often used in land-cover and land-use classification to help distinguish between vegetation types. However, a quarter of the top twenty important bands in this study were red edge whereas none of the spectral indices were ranked highly. This suggests that the presence of the red edge bands in the Sentinel-2 satellites might render the use of vegetation indices obsolete in boreal landscapes. Thus, there is need for more studies that assess the efficacy of Sentinel-2 red edge bands relative to vegetation indices for the purposes of land-cover and land-use classification, particularly in regions dominated by forests.

The dominance of the imagery from May and July in the band scoring metrics of this study raises important questions about the use of mono-temporal vs. multi-temporal imagery in land-cover and land-use classification. It is presently unknown whether the performance of an optimally-timed mono-temporal Sentinel-2 scene is equivalent to, or better than, a multi-temporal stack in classifying land-cover and land-use over complex landscapes. More studies are needed to assess the added value of the Sentinel-2 red edge bands in both mono-temporal and multi-temporal experiments involving different vegetation classes.

Acknowledgements

Open access funding for this article was provided by Lund University. I would like to sincerely thank all three reviewers for their time and effort in providing detailed and constructive comments. Funding to undertake this study was provided by the Department of Physical Geography and Ecosystem Science, Lund University.

Disclosure statement

No potential conflict of interest was reported by the author.

Supplementary material

Supplemental data for this article can be accessed here.