Multi-Temporal Sentinel-2 Data in Classiﬁcation of Mountain Vegetation

: The electromagnetic spectrum registered via satellite remote sensing methods became a popular data source that can enrich traditional methods of vegetation monitoring. The European Space Agency Sentinel-2 mission, thanks to its spatial (10–20 m) and spectral resolution (12 spectral bands registered in visible-, near-, and mid-infrared spectrum) and primarily its short revisit time (5 days), helps to provide reliable and accurate material for the identiﬁcation of mountain vegetation. Using the support vector machines (SVM) algorithm and reference data (botanical map of non-forest vegetation, ﬁeld survey data, and high spatial resolution images) it was possible to classify eight vegetation types of Giant Mountains: bogs and fens, deciduous shrub vegetation, forests, grasslands, heathlands, subalpine tall forbs, subalpine dwarf pine scrubs, and rock and scree vegetation. Additional variables such as principal component analysis (PCA) bands and selected vegetation indices were included in the best classiﬁed dataset. The results of the iterative classiﬁcation, repeated 100 times, were assessed as approximately 80% median overall accuracy (OA) based on multi-temporal datasets composed of images acquired through the vegetation growing season (from late spring to early autumn 2018), better than using a single-date scene (70%–72% OA). Additional variables did not signiﬁcantly improve the results, showing the importance of spectral and temporal information themselves. Our study conﬁrms the possibility of fully available data for the identiﬁcation of mountain vegetation for management purposes and protection within national parks.


Introduction
Mountain vegetation is particularly vulnerable to climate change, where the changes of the tree line and plant floor borders become visible [1]. The occurrence of species from various geographical regions in a relatively small area, often being glacial relics, endemics, or endangered species, makes their identification and monitoring extremely important for preserving natural wealth [2]. To achieve this, it is important to provide up-to-date vegetation maps of mountain protected areas.
Despite its high precision, field mapping requires a lot of time and work. In the case of high-mountain vegetation, the limited availability and shorter vegetation period compared to lowlands significantly affect the possibilities of field research. Due to rapid technological progress, remote sensing data, characterized by both greater objectivity and spatial coverage, are increasingly used [3]. The electromagnetic spectrum registered by remote sensing instruments, which create unique spectral characteristics of the analyzed objects, can support traditional methods of vegetation mapping by the use of image classification [3].
Recently, non-parametric classifiers are increasingly employed in vegetation classification [4][5][6][7][8] because of their more flexible approach to training data use than in parametric classifiers,

Study Area and Object of the Study
The study area, covering the highest parts of the Giant Mountains, was located at the Polish-Czech border, within the Krkonoše National Park and Krkonošský Národní Park, respectively; see Figure 1. There were five distinguished plant floors: foothills (up to 500 m above sea level a.s.l), lower (500-100 m a.s.l.) and upper montane zone (1000-1250 m a.s.l), subalpine (1250-1450 m a.s.l.), and alpine zone (above 1450 m a.s.l.). High mountain vegetation was the subject of our analysis (Figure 2), including the one located above the tree line, i.e., the areas located over 1250 m a.s.l. which were considered arctic-alpine tundra [7,47].
Remote Sens. 2020, 12, x FOR PEER REVIEW 3 of 29 accuracy obtained with these variables [30,40], but some authors also yielded a lower accuracy [6,41]. In mountain and upland areas, the use of digital terrain model (DTM) derivatives seems reasonable; however, predominantly using optical data requires an additional data source, like from airborne laser scanning (ALS) [42] or radar missions [31,43,44]. When using big datasets consisting of multiple scenes or/and additional bands, variable importance analysis is useful to select the best features affecting the accuracy, which was presented by many authors [29,30].
The main aim of this study was to assess the potential of Sentinel-2 multi-temporal data for vegetation types classification in the mountain ecosystem. Our assumption was to fully exploit plant phenology differences, allowing us to separate vegetation types between them, and as such, the multi-temporal classification approach was applied. Previous studies of the parts of this area employed different remote sensing images for vegetation classification [7,15,18,45,46]; however, no combined multi-temporal dataset was used for this purpose. In the context of the sensitivity of mountain ecosystems to climate and other changes, this approach seems to be particularly important, taking into account the variability of vegetation within the growing season, which may not be captured on a single image. We compared single-date and multi-temporal Sentinel-2 data, as well as the best combination of dates with calculated PCA bands and with vegetation indices to select the most optimal dataset. Based on reference polygons and prepared datasets, we performed iterative classification using the SVM algorithm and we assessed the obtained accuracies.

Study Area and Object of the Study
The study area, covering the highest parts of the Giant Mountains, was located at the Polish-Czech border, within the Krkonoše National Park and Krkonošský Národní Park, respectively; see Figure 1. There were five distinguished plant floors: foothills (up to 500 m above sea level a.s.l), lower (500-100 m a.s.l.) and upper montane zone (1000-1250 m a.s.l), subalpine (1250-1450 m a.s.l.), and alpine zone (above 1450 m a.s.l.). High mountain vegetation was the subject of our analysis (Figure 2), including the one located above the tree line, i.e., the areas located over 1250 m a.s.l. which were considered arctic-alpine tundra [7,47]. The largest areas of the entire Giant Mountains were covered by a mosaic of dwarf pine and grasslands (1500 ha on the Polish side and 3200 ha on the Czech side), playing an important ecological role by protecting the lower forests from snow avalanches [47]. Dwarf pine thickets were dominated by Pinus mugo species. Between their patches, there were shrub communities, consisting The largest areas of the entire Giant Mountains were covered by a mosaic of dwarf pine and grasslands (1500 ha on the Polish side and 3200 ha on the Czech side), playing an important ecological role by protecting the lower forests from snow avalanches [47]. Dwarf pine thickets were dominated by Pinus mugo species. Between their patches, there were shrub communities, consisting of bilberry Remote Sens. 2020, 12, 2696 4 of 24 with a dominance of Calluna vulgaris, Vaccinium myrtillus, and co-occurring Calamagrostis villosa and Deschampsia flexuosa grasses. Within the postglacial cirques, very rich subalpine deciduous shrubs were developed, where there was a regional scale unique community of endemic Salix lapponum willow complex [48]. Bottoms of cirques with fertile soils were covered by herbs from Adenostyles alliariae and Athyrietum distentifoli alliances. The significantly angled slopes above the valleys and postglacial cirques were dominated by species-rich grasslands, where the most popular were Calamagrostis villosa species. A unique feature of the Giant Mountains landscape were the subalpine-subarctic mountain peat bogs. The transition between subalpine and alpine zones was largely dominated by floristically poor communities of grasslands with the domination of Nardus stricta. In the alpine zone, the mountain grasslands are dominated by Juncus trifidus species. Rocks and screes in the upper parts of the Giant Mountains were an excellent habitat for species-rich vegetation and epilytic lichens from Umbilicarion and Rhizocarpion alliances.
Both sides of the Giant Mountains were affected by various factors, including agricultural activities, such as past deforestation and grazing, tourism once the parkland was classified as a protected national park, and sporadic avalanches and debris flows [7,15]. Additionally, in the 1980s and 1990s, an ecological disaster occurred due to strong winds and air pollution [49]. Other vegetation disturbances were caused by pests. The dynamic of vegetation changes were also noted due to expansive species like Molinia caerulea or Calamagrostis villosa encroachment and spreading [50]. All of these factors require a constant need to monitor the vegetation of this valuable mountain area.
Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 29 of bilberry with a dominance of Calluna vulgaris, Vaccinium myrtillus, and co-occurring Calamagrostis villosa and Deschampsia flexuosa grasses. Within the postglacial cirques, very rich subalpine deciduous shrubs were developed, where there was a regional scale unique community of endemic Salix lapponum willow complex [48]. Bottoms of cirques with fertile soils were covered by herbs from Adenostyles alliariae and Athyrietum distentifoli alliances. The significantly angled slopes above the valleys and postglacial cirques were dominated by species-rich grasslands, where the most popular were Calamagrostis villosa species. A unique feature of the Giant Mountains landscape were the subalpine-subarctic mountain peat bogs. The transition between subalpine and alpine zones was largely dominated by floristically poor communities of grasslands with the domination of Nardus stricta. In the alpine zone, the mountain grasslands are dominated by Juncus trifidus species. Rocks and screes in the upper parts of the Giant Mountains were an excellent habitat for species-rich vegetation and epilytic lichens from Umbilicarion and Rhizocarpion alliances. Both sides of the Giant Mountains were affected by various factors, including agricultural activities, such as past deforestation and grazing, tourism once the parkland was classified as a protected national park, and sporadic avalanches and debris flows [7,15]. Additionally, in the 1980s and 1990s, an ecological disaster occurred due to strong winds and air pollution [49]. Other vegetation disturbances were caused by pests. The dynamic of vegetation changes were also noted due to expansive species like Molinia caerulea or Calamagrostis villosa encroachment and spreading [50]. All of these factors require a constant need to monitor the vegetation of this valuable mountain area.

Sentinel-2 Satellite Data
Satellite images with no/minimal cloud coverage from 31 May, 7 and 27 August, and 18 September 2018 were downloaded from the Copernicus Open Access Hub browser. As mountain vegetation physiognomy is variable throughout the growing season (i.e., in late summer/early autumn, alpine and subalpine grasslands begin to discolor, making them easier to distinguish from, e.g., subalpine tall forbs), our assumption was to select data from different phenological phases to capture the vegetative phase, blooming, and senescence of the plants; this data can be reliably collected via the high temporal resolution of Sentinel-2 sensor. However, because study area was located in the mountains, the cloud cover was a significant limitation, which made it impossible to include a given image in the analyses; as such, all available images from the period from March to November, where clouds did not cover more than 10% were taken into account-June and July images were not included, because there were no images available that matched the criteria.

Sentinel-2 Satellite Data
Satellite images with no/minimal cloud coverage from 31 May, 7 and 27 August, and 18 September 2018 were downloaded from the Copernicus Open Access Hub browser. As mountain vegetation physiognomy is variable throughout the growing season (i.e., in late summer/early autumn, alpine and subalpine grasslands begin to discolor, making them easier to distinguish from, e.g., subalpine tall forbs), our assumption was to select data from different phenological phases to capture the vegetative phase, blooming, and senescence of the plants; this data can be reliably collected via the high temporal resolution of Sentinel-2 sensor. However, because study area was located in the mountains, the cloud cover was a significant limitation, which made it impossible to include a given image in the analyses; as such, all available images from the period from March to November, where clouds did not cover more than 10% were taken into account-June and July images were not included, because there were no images available that matched the criteria. The data with 2A processing level were selected, and the quality of atmospheric correction was assessed by the use of 9 field-collected spectra for radiometrically stable areas (asphalt and concrete). Based on the root mean square error (RMSE), we calculated a consistency of 6.5% for the satellite data. Reflectance 20 m bands were resampled to 10 m in SNAP (ESA Sentinel Application Platform v2.0.2, http://step.esa.int, Brockmann Consult, Skywatch, Sensar, and C-S). To avoid errors caused by the absorption by water, the three 60 m atmospheric bands 1 (coastal aerosol), 9 (water vapor), and 10 (Cirrus) were excluded from the analysis, leaving the bands most commonly used in land applications [6,29-31].

Additional Variables Calculation
To improve the method, we decided to use only variables calculated based on one data type (Sentinel-2 images). Our assumption was to add these variables separately to the multi-temporal dataset with the highest accuracy, to check whether these variables improve the accuracy (see Section 2.2.2). The first type of variables were the PCA bands derived from a statistical method of linear data transformation, which determines the new main axis of the coordinate system along with the largest possible data variance by projecting variable values in a multidimensional space [35]. New, uncorrelated variables called principal components were then calculated using the ENVI 5.3 software (Harris Geospatial Solutions, Broomfield, CO, USA). Based on the correlation table of variables that make up the best dataset and then on eigenvalues, the first PCA bands were selected as the most informative. The second type of variables were commonly used vegetation indices that determine the condition of vegetation or canopy water content [52]. Based on availability and the literature, we selected 18 vegetation indices from the spectral resolution of Sentinel-2 data and calculated them using ENVI software (the list of indices presents Table A1 in Appendix A).

Multi-Temporal Datasets Creation
To investigate the effect of combining data from different growing seasons on the obtained accuracies, we stacked the images from four dates and additional variables into 13 combinations (Table 1). We tested each image separately (datasets 1-4) and combined them using two, three, and four different dates (datasets 5-15). Additionally, for the best result, we added calculated vegetation indices and PCA bands (datasets 16-17). Table 1. Datasets prepared for classification process (the_best indicates the best-classified dataset assessed using overall accuracy (OA). X-the total number of bands in the best dataset; Y-the number of dates that make up the best dataset; Z-the number of principal component analysis (PCA) bands selected as the most informative). Due to clouds and topographical shadows on individual images, a mask was needed that would allow for the correct interpretation of the result on the final map. The cloud and shadow mask was created based on spectral values from 2 (blue) and 8 (NIR) bands, respectively, and the water vapor map (WVM) product was used in the mask development process as one of the auxiliary materials. The spectral similarity of deep shadows and water surfaces resulted in lakes also being masked. Additionally, the mask defining the range of the study, corresponding to areas above 1250 m a.s.l., was created based on a digital terrain model with a spatial resolution of 1 m derived from georeferenced point cloud from ALS data, filtered and classified in LAStools software (Rapidlasso, GmbH, Gilching, Germany). All masks were developed and applied to the images in ENVI 5.3 software (Harris Geospatial Solutions, Broomfield, CO, USA) using the Build Mask and Apply Mask functions.

Reference Data
On-ground vegetation type polygons were acquired using GPS Trimble GeoXT receiver (Trimble Inc., Sunnyvale, CA, USA) from 20-30 August 2013 and 30-31 August 2014 field campaigns in the Polish part of the Giant Mountains. As a reference legend, we used the botanical non-forest vegetation map created by Wojtuń andŻołnierz [47], which contains two levels of vegetation organization-communities and types. It was important to check the consistency between the remote sensing and reference data. We decided to use a vegetation type unit because some of the vegetation communities took up very small patches (e.g., 9 m 2 ), which would make it impossible to correctly designate them as training data on the Sentinel-2 pixel. Groups of vegetation types forming communities were more representative in this context. To fit the spatial size of Sentinel-2 data, we sampled field polygons registered in Projected Coordinate System Universal Transverse Mercator_(UTM, zone_33N) in patches equal to or greater than 20 × 20 m per vegetation type, based on the satellite data grid. We also updated on-ground information by checking the consistency between high-resolution Airborne Prism EXperiment (APEX) image (3.12 m, used bands: 640, 547, and 471 nm) from 2012 and PlanetScope image (3 m, used bands: 590-670, 500-590, and 455-515 nm) from 2018, both acquired in September. The visual analysis allowed us to conclude that there were no changes in the polygons within the study area, which, with the adopted 10 m resolution, could indicate that resources are outdated. Based on them, homogeneous patches were identified visually and translated into polygons representing vegetation types (Table 2).

Classification with Iterative Accuracy Assessment
All stacked datasets were classified using the SVM algorithm with the 'e1071' package [53] of R software [54]. We employed this algorithm due to the non-parametric character allowing for flexibility of training data used, high classification accuracy, and limited classification errors confirmed by many authors [9,12,13,15]. The most commonly used are linear and radial kernel functions, expressed as, respectively [55]: Remote Sens. 2020, 12, 2696 where x i x j are the feature vectors, and γ is the gamma parameter.
To select the best SVM parameters, including kernel function, γ, and cost of the penalty (C), tuning of parameters was performed using the 'tune.svm' function [53].
To ensure the contribution of each observation in the classification, we decided to use a 100-times iterative procedure of classification and validation of the results following Ghosh et al. [56]. We split the reference polygons into 60:40 for training and validation using a stratified random sampling approach, which allowed us to assess the variability within classes.
Quantitative assessment was based on an error matrix with the OA, as well as producer (PA), user (UA) [57], and F1 [58] accuracies for individual classes. After performing 100 repetitions of classification, the median of each measure was calculated for all datasets, and the distribution of results was presented as a boxplot. The resulting map was developed based on a dominant value calculated from 100 iterations of the best classified dataset.
Additionally, to assess the importance of individual bands from corresponding terms of data acquisition, from the entire ABCD dataset, we calculated variable importance using a receiver operating characteristics (ROC) curve analysis conducted for each vegetation type class using the "caret" package [59], which was based on an SVM model and integrated bootstrapping was repeated 100 times-this was the same procedure we performed for the best dataset with added vegetation indices and PCA bands. The results of these analyses, limited to the 20 most important variables from each dataset, can be seen in Figures A2-A4 in Appendix A.

Selection of the best Parameters and Dataset
The automatic manipulation of algorithm parameters provided the best configuration possible. Table 3 presents the accuracies obtained for four single-date scenes using default and optimized parameters of the SVM. Table 3 shows that the highest OA was achieved when the radial function was applied. The rest of the datasets proceeded using only optimized parameters, increasing the efficiency of the classification process. The classification of vegetation types was successfully derived from the Sentinel-2 time series, which, as mentioned prior, was better than using single-date data ( Figure 3). The OA of the classification for each of the multi-temporal datasets, obtained from the 100-times iterative procedure of accuracy assessment, was at least two p.p. (percentage points) higher than the results obtained for single-date datasets and ranged from 76.3% to 79.5%. The highest accuracy was achieved for the ABC dataset (the confidence interval (95%) for the classification was (0.7789, 0.8098)). Based on this, PCA bands and vegetation indices were included in this dataset. The correlation table of variables calculated for PCA bands selection showed the lowest correlations (less than 0.35) observed for bands 5, 6, 7, and 8 from A; 6, 7, and 8 from B; and 6, 7, and 8 from dataset C ( Table A2 in  of the eigenvalues allowed for the selection of the first 10 bands of PCA transformation to be the most informative for further processing, as they contained 99.4% of the total variance ( Figure A1 in Appendix A); as such, two combinations with additional variables consisted of 40 bands (30 spectral bands + 10 PCA bands) and 84 bands (30 spectral bands + 18 indices calculated for the images acquired on three dates). The combinations comprised of PCA bands and vegetation indices resulted in a lower OA: 77.1% and 79.2%, respectively ( Figure 3). classification for each of the multi-temporal datasets, obtained from the 100-times iterative procedure of accuracy assessment, was at least two p.p. (percentage points) higher than the results obtained for single-date datasets and ranged from 76.3% to 79.5%. The highest accuracy was achieved for the ABC dataset (the confidence interval (95%) for the classification was (0.7789, 0.8098)). Based on this, PCA bands and vegetation indices were included in this dataset. The correlation table of variables calculated for PCA bands selection showed the lowest correlations (less than 0.35) observed for bands 5, 6, 7, and 8 from A; 6, 7, and 8 from B; and 6, 7, and 8 from dataset C ( Table A2 in the Appendix A). The analysis of the eigenvalues allowed for the selection of the first 10 bands of PCA transformation to be the most informative for further processing, as they contained 99.4% of the total variance ( Figure A1 in Appendix A); as such, two combinations with additional variables consisted of 40 bands (30 spectral bands + 10 PCA bands) and 84 bands (30 spectral bands + 18 indices calculated for the images acquired on three dates). The combinations comprised of PCA bands and vegetation indices resulted in a lower OA: 77.1% and 79.2%, respectively ( Figure 3).

Vegetation Types Classification Results
The best-classified types of vegetation turned out to be forest and subalpine dwarf pine scrubs. Both classes reached a median of PA, UA, and F1 accuracies of around 90%-i.e., 94%, 89% ( Figure  4), and 90% ( Figure 5), respectively. Due to the high internal homogeneity, these classes also represent the greatest stability, which is confirmed by the small width of the distribution of achieved accuracy. Satisfactory results, above 60% of both accuracies and an F1 score above 70%, were also obtained for classes whose spectral characteristics constituted a specific mixture of signals, i.e., rock and scree vegetation and those characterized by unique physiognomy during subsequent stages of the growing season, i.e., grasslands. Deciduous shrubs were the worst classified type of vegetation,

Vegetation Types Classification Results
The best-classified types of vegetation turned out to be forest and subalpine dwarf pine scrubs. Both classes reached a median of PA, UA, and F1 accuracies of around 90%-i.e., 94%, 89% (Figure 4), and 90% ( Figure 5), respectively. Due to the high internal homogeneity, these classes also represent the greatest stability, which is confirmed by the small width of the distribution of achieved accuracy. Satisfactory results, above 60% of both accuracies and an F1 score above 70%, were also obtained for classes whose spectral characteristics constituted a specific mixture of signals, i.e., rock and scree vegetation and those characterized by unique physiognomy during subsequent stages of the growing season, i.e., grasslands. Deciduous shrubs were the worst classified type of vegetation, for which the median for both accuracies was less than 50% and F1 did not exceed 60%. For some classes, PA was lower than UA, because it was more difficult to fit the result with the validation data, e.g., for openwork and heterogeneous deciduous shrub vegetation class. For homogeneous classes like subalpine dwarf pine scrubs, this fit was simpler, and PA exceeded that of UA.
The analysis of the error matrix shows that the most frequently confusing class was deciduous shrub vegetation, which, in approximately 30% of cases, was classified as heathlands and in approximately 20% as subalpine tall forbs (Table 4). Physiognomically similar classes were confused quite often too-about 21% of the subalpine tall forbs' areas were confused with grasslands and nearly 14% with heathlands.
The visual assessment of the final map confirms the accuracy obtained during the quantitative evaluation ( Figure 6). The map of high-mountain vegetation created in the process of classification is largely similar to the reference map of non-forest vegetation [47]. The main differences concern the classes that obtained the lowest accuracy and which, because they are composed of heterogeneous complexes of communities with dominant species that are also species accompanying the other Remote Sens. 2020, 12, 2696 9 of 24 separations, have been mixed-i.e., the sites of heathlands, subalpine tall forbs, and deciduous shrub vegetation.
Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 29 for which the median for both accuracies was less than 50% and F1 did not exceed 60%. For some classes, PA was lower than UA, because it was more difficult to fit the result with the validation data, e.g., for openwork and heterogeneous deciduous shrub vegetation class. For homogeneous classes like subalpine dwarf pine scrubs, this fit was simpler, and PA exceeded that of UA.  The analysis of the error matrix shows that the most frequently confusing class was deciduous shrub vegetation, which, in approximately 30% of cases, was classified as heathlands and in approximately 20% as subalpine tall forbs (Table 4). Physiognomically similar classes were confused quite often too-about 21% of the subalpine tall forbs' areas were confused with grasslands and nearly 14% with heathlands. Table 4. Error matrix generated for the ABC dataset (based on iteration with OA closest to the median; SDPS-subalpine dwarf pine scrubs; F-forest; G-grasslands; NV-non-vegetation; for which the median for both accuracies was less than 50% and F1 did not exceed 60%. For some classes, PA was lower than UA, because it was more difficult to fit the result with the validation data, e.g., for openwork and heterogeneous deciduous shrub vegetation class. For homogeneous classes like subalpine dwarf pine scrubs, this fit was simpler, and PA exceeded that of UA.  The analysis of the error matrix shows that the most frequently confusing class was deciduous shrub vegetation, which, in approximately 30% of cases, was classified as heathlands and in approximately 20% as subalpine tall forbs (Table 4). Physiognomically similar classes were confused quite often too-about 21% of the subalpine tall forbs' areas were confused with grasslands and nearly 14% with heathlands. Table 4. Error matrix generated for the ABC dataset (based on iteration with OA closest to the median; SDPS-subalpine dwarf pine scrubs; F-forest; G-grasslands; NV-non-vegetation; Figure 5. F1 accuracy obtained during the 100-times iterative procedure of accuracy assessment for the best dataset.  G  1  0  311  1  11  2  0  40  25  391  NV  1  11  31  132  1  3  1  0  14  194  BF  4  0  27  3  118  0  0  7  3  162  RSV  0  0  3  6  0  127  2  3  2  143  DSV  2  9  0  1  0  1  27  1  24  65  STF  0  1  14  9  16  4  12  115  41  212  H  7  8  28  2  16  3  18  26  158  266  848  356  420  176  196  168  60  192  272  2416 Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 29

Discussion
Previous studies highlighted the advantages of single-date imagery in their vegetation classification for this mountain area [7,18,46]; however, our study goes further and demonstrates the advantage of using the freely available Sentinel-2 multi-temporal data for accurate vegetation mapping. We divided the discussion section into three parts: the first is devoted to mountain vegetation classification with special attention to the Giant Mountains study area (Section 4.1); in the second, we describe the use of multi-temporal data in classification (Section 4.2); in the third, we discuss the sense of including additional variables in the classification (Section 4.3).

Mountain Vegetation Classification
Image remote sensing has a great deal of potential for identifying types of mountain vegetation due to a wide range of available resolutions-spectral, spatial, and temporal, a properly generalized legend, and the selection of a classification algorithm adequate to the character of the data. The results obtained in this study (OA equal to 79.5%) show the usefulness of multispectral satellite data for the identification of types of high mountain vegetation. Other studies, which used comparable data, corresponding categories, and classification algorithms produced similar results, thus, confirming our results. Suchá et al. [45] classified eight vegetation types above the tree line in the Krkonoše Mts. National Park using Landsat-8 data (spatial resolution 30 m, seven spectral bands) and three per-pixel algorithm classifiers (ML, SVMs, and NNs) and obtained the best OA for ML classifier (78.3%). The study in [18] performed a similar investigation to ours, where eight types of vegetation were classified using simulated Sentinel-2 data and SVM classifier instead, which resulted in 81.9% OA for the dataset consisting of six bands as an effect of PCA and 78.3% OA for the total set of bands. Analogous studies with multispectral Sentinel-2 data were provided by Kupková et al. [7] where SVM, MLC, and NN algorithms were used to classify eight vegetation classes in Eastern tundra of the Giant Mountains, which yielded lower OAs than our results-71.0% and 79.5% based on SVM, respectively.
The complexity of the mountain vegetation means that the classification at a higher level of detail than the vegetation type requires the involvement of sensors registering images in many spectral bands, while also maintaining high spatial resolution. For this purpose, aerial hyperspectral data, such as AVIRIS (Airborne Visible/Infra-Red Imaging Spectrometer; 224 spectral bands), DAIS 7915 (Digital Airborne Imaging Spectrometer, 79 bands), AISA Dual (Airborne Imaging Spectrometer; 494 bands), or APEX (288 bands) data were employed [7,15,19]. These attributes allow to classify heterogeneous mountain vegetation at the community level and obtain high accuracy (74%-84% OA), despite its complicated structure; however, multispectral Sentinel-2 data also has the advantage of spectral resolution due to existence of SWIR, NIR, and red-edge bands, which was confirmed by a study of mountain vegetation [10] and other studies [28,30,32]. In our case, SWIR (11 and 12) and NIR (8a) bands were the most important in the classification of the entire ABCD dataset, which confirms the 20 first variables ( Figure A2 in Appendix A). SWIR1 was the first most important variable for deciduous shrub vegetation; SWIR2 for grasslands, bogs and fens, heathlands, and subalpine tall forbs; and NIR (8a band) for forests and rock and scree vegetation. For subalpine dwarf pine scrub classification, the most important was NIR (7 band). Depending on the class, different bands were placed in the 20 most important bands in classification; however, for all classes except subalpine dwarf pine scrubs, SWIR bands from four different dates always occurred. To improve the classification model all used spectral bands could be correlated with each other to select only uncorrelated ones for further study. However, in most of the articles that discuss the use of SVM for vegetation classification, all available Sentinel bands (except the so-called "atmospheric" bands) are used [7,60,61].
The complex character of the classified vegetation types causes divergent results. Large-area forests and subalpine dwarf pine scrubs growing in homogeneous patches turned out to be the most identifiable classes, reaching medians of PA and UA above the 90% (forest: PA-95.0%; UA-96.5%; subalpine dwarf pine scrubs: PA-95.2%; UA-90.0%). The specific texture of the subalpine dwarf pine scrubs and the characteristic spectral reflection of mosaics of forest-forming species makes them one of the best-classified types of alpine vegetation. Similar results, where PA and UA accuracies fluctuated around 90%, were obtained by other authors classifying subalpine dwarf pine scrubs and forests on Landsat-8 [45], Sentinel-2 [7,18], the Environmental Mapping and Analysis Program (EnMAP) [18,46], APEX [7,15,46], and AISA Dual [7] data. Among the best-classifying types of vegetation, rock and scree vegetation was also frequently mentioned, which is well distinguishable even by sensors with a spatial resolution of no more than 10 m, i.e., Sentinel-2: PA-91.5%, UA-87.0% [18] and PA-92.7%, UA-95.0% [7] or EnMAP: PA-97.5%, UA-96.3% [18]. The clue is the fact of mixing the signal from plants and rocks, which results in a high reflectance in NIR and SWIR bands. In our work, it was the third best-classified type of vegetation, achieving a PA and UA of 83.0% and 88.1%, respectively. On the other hand, classes with heterogeneous community complexes proved to be the worst-classified types of vegetation, with an accuracy not exceeding 60%, whose dominant species are often also accompanying species of other types-i.e., heathlands (PA-60.3%, UA-55.2%) and subalpine tall forbs (PA-54.9%, UA-52.8%) and those that form too small clusters to be well distinguishable by the sensor-i.e., deciduous shrub vegetation (PA-25.6%, UA-41.0%). Subalpine tall forbs and heathlands, due to their high spectral similarity, practically prevent proper separation on multispectral data-in Landsat-8, the OA often does not exceed 50% [45], and in the case of Sentinel-2, the OA often does not exceed 60% [7,18]. Higher values are obtained for classifications based on hyperspectral data, where PA and UA reach a 70% of OA, i.e., EnMAP: PA-79.2%, UA-67.9% (subalpine tall forbs), and PA-44.8%, UA-52.0% (heathlands) [46], or AISA Dual: PA-85.1%, UA-85.8% (subalpine tall forbs), and PA-81.6%, UA-83.8% (heathlands) [7]. Table 2 shows that the datasets for classification were not balanced. It was also noted by other authors as a problem affecting the accuracy of machine learning [16,17]. Our results show this effect in the most and least numerous classes as subalpine dwarf pine scrubs and deciduous shrub vegetation, respectively; however, for forests and rock and scree vegetation mentioned above also occupied by smaller areas (closer to the least than the most numerous) the accuracies were some of the highest. Balancing the dataset, particularly in the natural vegetation area, is not as straightforward as confirmed by other authors [6, 29,30]; however, Thanh Noi and Kappas [17] proved that for SVM, in comparison to random forest (RF) and k-nearest neighbor (k-NN), the difference between various sample sizes was insignificant, which can support our method choice and obtained results, and we conclude that this effect was not very pronounced.

Multi-Temporal Classification
The main aspect working in favor of Sentinel-2 data, in addition to its open-access character, is the high temporal resolution enabling the generation of multi-temporal compositions. Taking into account the images in which successive stages of vegetation were captured, it is possible to obtain much better classification results compared to the results based on single-date data. The results obtained in this study confirm that each of the analyzed multi-temporal compositions, two, three, or four dates, resulted in a higher value of OA concerning the classification on single-date data-the highest value for single-date data was 74.2%, whereas the lowest value for multi-temporal composition was 76.3%. A similar tendency, regardless of the sensor and used algorithm, can also be observed for other objects of the study, i.e., land cover [62], tree species [28][29][30][31], swamp [25], and grassy [27,33] and shrubby [6] communities (Table 5).
Selecting the optimal composition that generates the highest accuracy is complex. The quality of the composition is influenced by both the number of images that comprise it and the date of registration. It was shown that at some stage, the addition of further images does not cause a further increase in accuracy, and the stabilization of the result depends largely on the dates taken into account [27,38]. In our study, the highest OA was obtained for a composition consisting of three images (79.5%) out of four available (78.5%). Similarly, in the study where shrub communities were classified, despite the access to four images (12% OA), a composition using only two of them resulted in a higher accuracy (68% OA) [6]. In the case of classification of tree species, where authors had access to 18 images, the composition with five images obtained the highest result-92.1% and 92.4% OA, respectively [29].   The key assumption in multi-temporal classification is that the vegetation varies between different terms of data acquisition within a year when using inter-seasonal data. When there is not enough spectral information to divide spectrally similar groups of vegetation, even with SWIR or NIR regions registered by Sentinel-2, the use of temporal information as additional variables demonstrates the advantage of the approach proposed in this work. The registration dates of the images constituting the multi-temporal composition are the key issue determining the quality of the obtained results. A satisfactory classification result largely depends on the characteristics of the studied vegetation and its phenological cycle. By capturing those periods of the year in which key stages of a plant life cycle are observable, it is possible to generate compositions that produce the best results. In most cases, compositions including the contrast of spring and autumn are considered to be the most informative because it is the time of intensified discoloration associated with flowering and senescence of vegetation [6, 28,29]. Slightly less often, especially in the case of grassy vegetation, early and late summer are also indicated, i.e., the period of dynamic growth [27,63]. In the case of high-mountain vegetation, which is the subject of this study, the composition generating the highest result (79.5%) was composed of spring (31 May) and late summer (07 August, 27 August) images. In this case, grasslands that discolor as a result of drying are the most important indicator for the separation of different types of vegetation [7,15,45]. Additional variable importance analysis performed for the ABCD dataset revealed that, for all classes, the most numerous important features located at the first 20 places were from A, B, and C datasets ( Figure A2 in Appendix A). The spectral band from spring was indicated as the most important for subalpine dwarf pine scrubs, grasslands, heathlands, and deciduous shrub vegetation; the band from the B dataset was the best in the classification of forests and rock and scree vegetation; and band from C dataset was the best in the classification for bogs and fens and subalpine tall forbs.

Additional Variables
Apart from different factors, e.g., sensor, legend generalization, algorithm, or the number of images analyzed, the impact on the obtained classification result also has additional processing. It includes those that reduce spectral space (transformations), and also those that increase the information capacity in the form of vegetation indices added as new bands. The transformation was performed to extract key information that differentiates the analyzed classes. Although in many studies this procedure resulted in an increase of the OA-AISA Eagle II from 72.8% to 82.1% (minimum noise fraction transformation; MNF) [9]; AISA Dual from 74.2% to 84.3%, APEX from 77.7% to 82.6% (PCA) [7]; Sentinel-2 from 78.3% to 81.9% (PCA) [18]-in our work, the PCA transformation reduced the final result by slightly more than 2 pp. (from 79.5% to 77.1%); however, this is not an exception, because there are also studies in which transformations led to lower accuracies-APEX from 82.7% to 81.0% (PCA) [15]; EnMAP from 82.9% to 56.3% (PCA) [18]; Sentinel-2 from 72.0% to 51.0% (PCA) [6]. Calculation of the most important variables in the classification of the ABC_PCA dataset allowed us to confirm this conclusion because only two/three PCA bands were placed in 20 top variables for each vegetation type classification ( Figure A3 in Appendix A). In the analysis of individual classes, in the first place, PCA bands were noted for rock and scree vegetation, bogs and fens, and subalpine tall forbs. Overall, from the 10 calculated PCA bands only 4 occurred in these 20 most important variables, and the most frequently occurring PC1 and PC2 bands.
When including vegetation indices as an input to the classification of multispectral data to increase its quality as, e.g., in forest species classification improving the model performance by around five percentage points [30], in our study it led to a decrease in the OA from 79.5% to 79.2%; however, in the variable importance analysis of the ABC_IND dataset, we noted that in most cases, they were located in the top-20 variables, and for each class, at the first place vegetation index was located ( Figure A4 in Appendix A). The most frequently occurring were EVI and VARI, based on VIS (visible) + NIR and only VIS spectral bands, respectively, but band 11 from SWIR was similarly frequent. For five classes, the first most important variable was ReNDVI (forest, rock and scree vegetation, grasslands, heathlands, and deciduous shrub vegetation) supporting the importance of red-edge in combination with the NIR band. Many publications described strategically selected spectral ranges of Sentinel-2 data [6, 28,41,64]. This multispectral satellite, registering the electromagnetic spectrum in 13 spectral ranges, has several narrow (less than 20 nm) bands. As a consequence, the only spectral bands themselves can generate higher results than sets enriched with secondary information from, e.g., vegetation indices. In a study where grassland species were the authors only analyzed the classification of spectral bands, they were able to reach 90.4% OA, whereas with NDVI, they were able to reach 88.6% [41]. Another research, describing shrubs classification, reported 72.0% OA for only spectral bands and 59.0% for a dataset consisting of spectral bands, NDVI, and PCA [6].
Based on our results and additional analysis of important variables in classification (Figures A2-A4 in Appendix A), this study can be investigated further to create the most optimal models with only uncorrelated spectral bands after correlation analysis and with the best variables after feature selection. However, as our study aim stated, we wanted to assess the potential of Sentinel-2 multi-temporal data for mountain vegetation types classification by analyzing plant phenology differences through the growing season; hence, here we decided to use whole datasets to recommend the best dates for creating a multi-temporal composition.

Conclusions
The presented problem of mountain vegetation types mapping in the Giant Mountains allowed us to determine the usefulness of Sentinel-2 multi-temporal satellite data. Analysis of the obtained results led to the following conclusions: 1.
Sentinel-2 multispectral data allow us to classify high-mountain vegetation at a satisfactory level of accuracy, assuming the right level of generalization of the legend, the selection of a classification algorithm adequate to the character of the data, and the use of the advantages associated with high temporal resolution-classification based on multi-temporal compositions allows achieving better results compared to the results generated based on single-date data. Contrary to high-accuracy hyperspectral data not fully available at this moment even for single-date collection and limited for use in local-scale analysis, Sentinel-2 data can be assessed as more applicable.

2.
The quality of the temporal composition, in addition to the number of images, is primarily due to the date of acquisition-compositions containing contrasting spring and autumn, i.e., the time of intensified discoloration associated with flowering and senescence vegetation, were considered to be the most informative. Lower OA of a single image does not exclude it as a valuable component of the multi-temporal composition, as after adding an image from late August gave better accuracies than the two preceding images from the beginning of August and the end of May.

3.
The additional variables (vegetation indices and PCA transformation bands) tested on the best-classified dataset did not contribute to the increase in OA, which suggests that in the case of the classification of multi-temporal Sentinel-2 data, the most important variables for a satisfactory result are the images themselves (number and dates of acquisition), not their additional processing; however, the inclusion of vegetation indices can be investigated more deeply, taking into account the most influential indices for particular vegetation types classification to build the models based on only the most informative features.
The map of mountain vegetation types in the Giant Mountains developed based on Sentinel-2 data is an objective source of information that can support monitoring works, especially because the high temporal resolution, ensuring access to constantly supplemented data resources, enables its continuous updating. Funding: This research received no external funding. The APC was funded by the Faculty of Geography and Regional Studies, University of Warsaw.

Acknowledgments:
The authors wish to thank EUFAR, DLR and VITO for APEX data acquisition in 2012. We are also grateful to Lidia Przewoźnik, Bronisław Wojtuń, and Bogdan Zagajewski for their support in field data collection and to both national parks authorities for giving the access to provide this research. We are also grateful to Reviewers and Editor, who allowed us to improve the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.     Figure A3. The 20 most important variables in classification of each vegetation type using ABC_PCA dataset.