Mapping Mediterranean Forest Plant Associations and Habitats with Functional Principal Component Analysis Using Landsat 8 NDVI Time Series

: The classiﬁcation of plant associations and their mapping play a key role in deﬁning habitat biodiversity management, monitoring, and conservation strategies. In this work we present a methodological framework to map Mediterranean forest plant associations and habitats that relies on the application of the Functional Principal Component Analysis (FPCA) to the remotely sensed Normalized Di ﬀ erence Vegetation Index (NDVI) time series. FPCA, considering the chronological order of the data, reduced the NDVI time series data complexity and provided (as FPCA scores) the main seasonal NDVI phenological variations of the forests. We performed a supervised classiﬁcation of the FPCA scores combined with topographic and lithological features of the study area to map the forest plant associations. The supervised mapping achieved an overall accuracy of 87.5%. The FPCA scores contributed to the global accuracy of the map much more than the topographic and lithological features. The results showed that (i) the main seasonal phenological variations (FPCA scores) are e ﬀ ective spatial predictors to obtain accurate plant associations and habitat maps; (ii) the FPCA is a suitable solution to simultaneously express the relationships between remotely sensed and ecological ﬁeld data, since it allows us to integrate these two di ﬀ erent perspectives about plant associations in a single graph. The proposed approach based on the FPCA is useful for forest habitat monitoring, as it can contribute to produce periodically detailed vegetation-based habitat maps that reﬂect the “current” status of vegetation and habitats, also supporting the study of plant associations.


Introduction
Phytosociology, today, plays key roles in European policies dedicated to the conservation of biodiversity, such as the 92/43/EEC Habitats Directive [1] and the Natura 2000 network. Plant associations and the other phytosociological levels, in fact, allow for the diagnosis of natural and seminatural habitats listed in Annex I of the Habitats Directive [2][3][4][5]. According to the Directive, each European member state is obliged to survey, monitor, and report about habitats, their distribution, and status of conservation every 6 years. The phytosociological maps (e.g., plant associations map) can adequately represent the distribution of the habitats and, if updated and repeated, would be helpful in evaluating and monitoring the status of conservation of the habitats [6][7][8]. These activities, to be performed every 6 years, require a lot of resources (i.e., field survey, photo-interpretation) that is constraint making the update of habitat maps particularly difficult due to costs and required time. Ichter et al., [9,10] reported that in Europe, after the adoption of the Habitats Directive in 1992, habitat project mapping based on phytosociological approach maps has increased, but very few maps have been updated and repeated. Discrimination of contiguous plant associations or upper syntaxa through the visual interpretation of remote images is based on subjective decisions and could be time and labor intensive. Moreover, the accuracy of these maps is often neglected or absent [10,11].
Remote sensing is an important and complementary source of information in field surveys that contributes to a better understanding of the diversity of natural and seminatural habitats, their spatial distribution, and their conservation status [12,13]. It provides multiple advantages over traditional mapping (i.e., photo-interpretation and field surveys), such as faster mapping production that facilitates repeatability even in inaccessible areas [11]. The increasing availability of remote sensing data free of charge of high spatial and temporal resolution (e.g., Sentinel-2 and Landsat-8) offers dense multi-temporal measures of greenness, such as the Normalized Difference Vegetation Index (NDVI) [14] times series, that are a useful proxy for the seasonal and inter-annual vegetation phenological changes [15]. These changes have proved useful in the discrimination of contiguous Mediterranean habitats [16][17][18] and to produce phenology-based mapping of the vegetation [19,20].
We hypothesize that Mediterranean forest plant associations have distinct phenological behavior during the year (NDVI time profile). Moreover, the application of the Functional Principal Component Analysis (FPCA) on remotely sensed NDVI enables the identification of the main seasonal phenological variations-FPCA scores-that are effective spatial predictors for plant associations mapping. Therefore, we present a novel methodological framework, as to the best of our knowledge, no previous study has used the FPCA for mapping forest plant associations supported by remotely sensed data to automatize the overall process. In this framework, our main objectives are: (i) identify the main seasonal variations (phenological) by applying the FPCA to the remotely sensed pixel-based NDVI time series; (ii) identify the main topographic and lithological gradients by applying the Principal Component Analysis (PCA) to the local terrain attributes (extracted from a Digital Elevation Model (DEM)) and the lithology; and finally, (iii) perform a supervised classification of the predictors to discriminate and map the Mediterranean forest plant associations and assess their importance in contributing to the final accuracy of the map. The methodology here presented aims to demonstrate the capability to generate high accuracy habitat maps based on phytosociological data with a reduced human effort due to automation of several processes.
average height of trees of 14.5 ±1.5m) are widespread. In the 1930s, the western sectors of Mount Conero were in a very degraded hydrogeological state and, for this reason, they were reforested [26].
As Cyclamino hederifolii-Quercetum ilicis and Cephalanthero longifoliaee-Quercetum ilicis belong to the habitat 9340, "Quercus ilex and Quercus rotundifolia forests" [5], and the area falls within the Natura 2000 network, the vegetation field plot surveys and the habitats mapping based on plant associations must be repeated every 6 years [1,27]. Nevertheless, because of the "traditional" mapping approach adopted and the complex spatial organization of the plant associations determined by the joint action of the topographic variability with biotic (e.g., dynamic vegetation processes) and anthropic factors, the production of vegetation maps has been irregular and discontinuous to date. The latest "traditional" plant association-based habitat map was carried out in the 2005 and was made available by the "Web Geographic Information System for Fauna, Flora, and Plant Landscape Data Management of Marche region" (http://sitbiodiversità.ambiente.marche.it/sitrem/) [7]. For all these reasons, the study area is suitable to test the proposed mapping approach. Figure 1. Study area. In green are the Mount Conero (Marche, Central Italy) forests investigated. The distribution is derived from the "Web Geographic Information System of Marche region" [7]. Red points are the reference data (field survey and photo-interpretation).

Remote-Sensed NDVI Times Series
The extent of the forest patches, according to O'Neill et al. [28], fit the spatial resolution of the Landsat 8 remote-sensed images (30 meters) here adopted. We looked for all cloud-free Landsat 8 Figure 1. Study area. In green are the Mount Conero (Marche, Central Italy) forests investigated. The distribution is derived from the "Web Geographic Information System of Marche region" [7]. Red points are the reference data (field survey and photo-interpretation).
As Cyclamino hederifolii-Quercetum ilicis and Cephalanthero longifoliaee-Quercetum ilicis belong to the habitat 9340, "Quercus ilex and Quercus rotundifolia forests" [5], and the area falls within the Natura 2000 network, the vegetation field plot surveys and the habitats mapping based on plant associations must be repeated every 6 years [1,27]. Nevertheless, because of the "traditional" mapping approach adopted and the complex spatial organization of the plant associations determined by the joint action of the topographic variability with biotic (e.g., dynamic vegetation processes) and anthropic factors, the production of vegetation maps has been irregular and discontinuous to date. The latest "traditional" plant association-based habitat map was carried out in the 2005 and was made available by the "Web Geographic Information System for Fauna, Flora, and Plant Landscape Data Management of Marche region" (http://sitbiodiversita.ambiente.marche.it/sitrem/) [7]. For all these reasons, the study area is suitable to test the proposed mapping approach.

Remote-Sensed NDVI Times Series
The extent of the forest patches, according to O'Neill et al. [28], fit the spatial resolution of the Landsat 8 remote-sensed images (30 m) here adopted. We looked for all cloud-free Landsat 8 Level-2 images (that had already been preprocessed with absolute radiometric correction) available for the period April 2013 (when Landsat 8 started)-April 2019. These six-year data are suitable for the 4th report ex Art. 17 (Table A2).
We co-registered the images, normalized the topographic effect [29] and then applied the NDVI index (by the R "RStoolbox" package [30]). Before analyzing the intra-annual variations (which represent the seasonality of the vegetation) we checked the inter-annual variation through a harmonic regression of the monthly median NDVI values of the 91 scenes collected (see Figure 4a). In addition to an obvious annual periodic variation, this analysis showed for the study area slight biennial and/or triennial variations. These slight variations allow consideration of the years of the period as useful replicas to describe the average intra-annual variations for the investigated period. Any land use changes that have occurred on small surfaces (e.g., coppicing, fire), which cannot be captured by harmonic regression, are detected by FPCA (see Section 2.2.2). Then, the NDVI surfaces, chronologically ordered according to the Day of Year (DoY), were collected in a raster stack. A preliminary analysis of the stack showed that the distribution of the acquisition dates of the scenes indicated by the DoY is irregular and not very dense (87 days out of 365). The same dates converted into weeks (weeks) are denser (40 weeks out of 53 total), and practically provide a complete coverage over the study period if converted into periods of 14 days (bi-week). In this last case, in fact, there are NDVI values in 25 bi-weeks out of the 26 totals (only bi-week 21 is missing). Therefore, after checking and removing the anomalous values (function clean.ts() of the R package "forecast"; [31,32]), we aggregate the DoY NDVI values to bi-weekly mean values and obtained a (one-year) bi-weekly pixel-based NDVI times series that represented the mean intra-annual seasonality variation of the forests.

Functional Principal Component Analysis of NDVI Time Series
FPCA is a tool of functional data analysis [33]. FPCA, explicitly accounts for the chronology of data (here bi-weeks) treating the whole NDVI curve as the statistical unit. Indeed, the classic PCA is not appropriate for time series analysis because bi-weeks would be considered to be independent vectors of values [34]. FPCA decomposes the space of functional data (e.g., NDVI time series) into low orthogonal functional principal components [35]. FPCA provides: (i) eigenvalues that account for the variation explained by each component; (ii) the vector of principal component scores (FPCA scores) that summarizes the similarity between the sampling units (here the pixel-based NDVI time series) and (iii) the eigenfunctions (which replace the eigenvectors of the classical PCA) that represent the principal "modes of variation" of the data [34]. To identify the main modes in which the forests vary during the year, we applied the FPCA to the mean bi-weekly pixel-based NDVI times series. We used the nonparametric method called "principal components analysis through the conditional expectation" (PACE) [36], implemented in the R "fdapace" package [37]. The pixel-based FPCA scores, which express the main modes in which forests vary during the year, were central to the subsequent analysis. They were used as input variables (predictors) for the supervised plant association mapping and for relationship analysis with lithological and topographic variables and plant associations. These scores are also useful to detect outliers (extreme FPCA scores indicate pixels with anomalous NDVI curves) that in some pixels (areas) could correspond to a change in the land use.

Topographic and Lithological Factors
Lithology and topography determine different environmental conditions that affect the forest plant association composition and distribution [25,38]. If combined with multi-temporal remote sensing data, Remote Sens. 2020, 12, 1132 5 of 22 they could improve the accuracy of habitat mapping [39,40]. Then, from a DEM with 30m resolution (the DEM has been derived from the NASA Shuttle Radar Topography Mission (SRTM) that ensures data 1 arc-seconds), we extracted a set of parameters by using the System for Automated Geoscientific Analyses (SAGA) GIS software [41]. These parameters are useful proxies of the main topo-climatic and hydrological variations [42] ( Table 1). The lithological information was derived from the lithological map of the Marche Region [43]. A standardized PCA (using the R package "RStoolbox"; [30]) of the terrain parameters and the binarized class lithology was applied. The pixel-based PCA scores, which express the main environmental (topographic and lithological) gradients were used as input variables (predictors) for the supervised plant association mapping.

Forest Plant Associations and Reference Data Collection
The collected reference data are based on field observations (phytosociological analysis) subsequently integrated with observations extracted through visual interpretation of high-resolution satellite images.
Phytosociological analysis was performed in 2 steps: in the first, phytosociological relevés were performed in 53 circular plots of 200 square meters (see Figure A1 showing the variation of vegetation in the study area). This plot dimension, suited to forest plant associations [45], allowed us to center the plot by a Real Time Kinematic (RTK) Global Positioning System (GPS) receiver in the pixel. At each plot, all vascular plant species and their Braun-Blanquet abundance scale values were recorded [46]. The sites x species abundance matrix was processed by hierarchical cluster analysis (the chord distance, and the Unweighted Pair Group Method with Arithmetic mean (UPGMA) linkage). The final number of clusters (plant associations) was chosen on expert-based comparison with the literature [23][24][25]47]. The Indicator Species Analysis [48] was performed to identify the indicator species of the plant associations (phi coefficient > 0.4, p-value = 0.05). In the second step, another 27 sites (pixels), by a field expert-based similarity comparison, were assigned to the plant associations identified in the first step.
Finally, in order to increase the number of reference plots and obtain a balanced space distribution for the reference data, we added reference plots assigned to the plant associations by visual interpretation of the Google Earth imagery (Google Earth images, 26/04/2018 and 31/08/2017 with 0.5m resolution, see Figure A2). The total number of reference plots was 175; 80 were collected by phytosociological analysis and 95 were photo-interpreted (Section 2.3 provides more details about reference plots). The final reference data classification was used as training data for the supervised plant association mapping and as for relationships analysis with the main phenological gradients identified by FPCA.

Relationships between NDVI Seasonality, Topography, and Forest Plant Associations
The first two FPCA scores (FPCA1 and FPCA2) were used to define the reduced phenological ordination space of all pixel-based NDVI time series. The correlations with terrain parameters (vectors) and lithologies (categorical) were tested (envfit function; vegan package in R statistical software [49]) and those that were significant (alpha level p < 0.05) were overlaid onto the ordination space. To evaluate and characterize the phenological behavior of the plant association, the reference data classification was overimposed, as a spiderplot, onto the ordination and used to fit the plant association NDVI temporal profiles (curve) (by the CreatePathPlot function of the fdapace package [37]).

Supervised Mapping of the Forest Plant Associations
As shown in the workflow diagram ( Figure 2), the supervised random forest classification (RF) [50] of the pixel-based FPCA scores (which represent the main modes of forest NDVI variation during the year) and the pixel-based PCA scores (which represent the topographic and lithological gradients) were used to discriminate and map the forest plant associations. The reference data classification was used as training data.
We resampled the pixel size of the predictors at 10m by a bilinear interpolation to ensure the maximum correspondence between the reference and predictors data. Then we extracted and assigned (by a bilinear interpolation using the four nearest raster cells) the predictor values to the reference data.
RF is a tree-based classification method that aggregates the results of many trees (ntree) to make a prediction for estimating the model error. The trees were created by subset of the training data ("in-bag samples") with the remaining "out-of-bag samples" to estimate the model error, known as the out-of-bag (OOB) error. In addition, there was a random subset of candidate predictor variables (mtry). RF makes it possible to evaluate the predictor importance, by selecting and ranking those predictors with the greatest ability to discriminate between the target classes [51].
It is known for the area that the Asparago acutifolii-Ostryetum association is rather limited, and thus we expected the training data to be unbalanced. Statistical classifiers such as RF can be biased if the proportions of training and validation data are unbalanced. Final classification can show an over-prediction of the majority classes and an under-prediction of the minority class. Over-sampling and under-sampling are used to produce more balanced datasets [52]. We incorporated the down-sampling, which makes the frequency of the majority class closer to those of the rarest class, into the RF. We finally set RF as the following: ntree = 1500; subset of training data = 4 * nmin, where 4 is the number of target groups. This means that the RF model for 1500 can take a balanced random sample without loss of information for the majority classes. Finally, we check mtry from 2 to the square root of the number of input variables, as usually set [51]. As shown in the workflow diagram (Figure 2), the supervised random forest classification (RF) [50] of the pixel-based FPCA scores (which represent the main modes of forest NDVI variation during the year) and the pixel-based PCA scores (which represent the topographic and lithological gradients) were used to discriminate and map the forest plant associations. The reference data classification was used as training data. To calibrate the model and provide a robust estimate of the accuracy by limiting the bias, we repeated the 10-fold cross-validation 20 times. A mean overall accuracy (OA) index, Kappa index (and their respective Standard Deviations (SD)), and a cross-validated confusion matrix (representing the error distribution for class among the 20 repeats) were calculated. The importance of the predictors was evaluated and ranked by the normalized RF OOB Mean Decrease Accuracy index. The higher the values, the greater the relevance of the predictor. We run a RF model only with phenological variables (FPCA scores) to evaluate how topography and lithology contribute to the final OA (see Section 3.3). Furthermore, the map resulting from the application of the RF model obtained for all the pixels was visually interpreted and compared with the previous traditional maps (these analyses were performed using the packages "raster" [53] and "caret" [54]).

Forest Vegetation Field Survey and Plant Associations
In the first phase of field phytosociological analysis, the abundance values of 77 plant species were recorded in the 53 phytosociological relevés. The hierarchical cluster analysis of the relevés x species abundance data matrix provided the dendrogram with four principal clusters ( Figure 3).
predictors was evaluated and ranked by the normalized RF OOB Mean Decrease Accuracy index. The higher the values, the greater the relevance of the predictor. We run a RF model only with phenological variables (FPCA scores) to evaluate how topography and lithology contribute to the final OA (see sub-section 3.3). Furthermore, the map resulting from the application of the RF model obtained for all the pixels was visually interpreted and compared with the previous traditional maps (these analyses were performed using the packages "raster" [53] and "caret" [54]).

Forest Vegetation Field Survey and Plant Associations
In the first phase of field phytosociological analysis, the abundance values of 77 plant species were recorded in the 53 phytosociological relevés. The hierarchical cluster analysis of the relevés x species abundance data matrix provided the dendrogram with four principal clusters ( Figure 3). According to Biondi et al. [23][24][25]47], Cluster 1 corresponds to the Cephalanthero longifoliae-Quercetum ilicis association, Cluster 2 to the Cyclamino hederifolii-Quercetum ilicis, Cluster 3 to the Asparago acutifolii-Ostryetum carpinifoliae, while Cluster 4 is the plantation of Pinus halepensis and P. pinea. The frequency species composition and the indicator species of plant associations are summarized in Table A1. Through the subsequent field analyzes (step 2) and the visual interpretation of the high-resolution images ( Figure A2), we obtained the training data of the plant associations. The training dataset, as expected, is quite unbalanced where the Asparago acutifolii-Ostryetum carpinifoliae association is the minority class (Table 2). Table 2. Reference data used for training. Number of plot references for plant associations acquired by field phytosociological analysis (in parenthesis the number of plots analyzed by phytosociological relevés) and visual interpretation of high-resolution images.

Plant Associations Training Data
Phytosociological Analysis

Visual Interpretation of Images Total
Cephalanthero longifoliae-Quercetum ilicis 29 (22) Figure 4a shows the inter-annual neglectable variations that confirms an unchanged land use or the absence of extensive changes. Figure 4b shows the six-year period NDVI variations combined in a unique year according to the DoY and finally. Figure 4c shows the 7292 bi-weekly pixel-based NDVI time series obtained and analyzed through FPCA. Seven orthogonal FPCA components, which reflected the main contrasting modes of NDVI forest variation during the year, were identified. The first two FPCA components explain 90.7% of the total variation. FPCA components, in practice, represent the amount of deviation from the overall mean of the NDVI time series. The FPCA1 component accounted for 68.0% of the total variation. It opposes and orders the pixel-based bi-weekly NDVI times series that, in the autumn-winter period (months of January, February, and March, bi-weeks 1-8; months of October, November, and December, bi-weeks [18][19][20][21][22][23][24][25][26], have low NDVI values, to those with high values (Figure 5a,c). The FPCA2 component accounted for 22.7% of the total variation. It opposes and orders the pixel-based NDVI time series that, in the spring-summer period (from March to September; bi-weeks 9-19), have quite constant NDVI values to those that increase up to early August (18th bi-week) (Figure 5b,c).  the positive part of FPCA1 and FPCA2. Its NDVI profile shows a low NDVI seasonality. NDVI values slightly increase in the summer period up to the 17th bi-week (end of August), where they reach the maximum value (about 0.85); then, the NDVI values decrease slightly until the end of the year.

Main Seasonal NDVI Variation of the Forests
The Pinus halepensis and P. pinea plantation (purple spiderplot and NDVI profile in Figure 5c) has NDVI values ranging from 0.75 to 0.8. Unlike the other plant associations, it has the NDVI autumn-winter values greater than the spring and summer ones. This plantation is thus allocated in the positive part of FPCA1 and the negative part of FPCA2.

Relationships between NDVI Seasonality, Topography, and Forest Plant Associations
The first two FPCA scores (FPCA1 and FPCA2) allowed us to define the reduced phenological ordination space, where the significant correlated topographic and lithological variables were overlaid. FPCA1 is positively correlated with ELEV, TPI, SOL and the lithology (L1) while negatively with L4. FPCA2 is positively correlated with WIN, SLO, while negatively with SOL (Figure 5c). The plant associations, superimposed as spiderplots, show a distinct behavior during the year, clearly represented by the respective NDVI fitted temporal profiles (Figure 5c).
The Asparago acutifolii-Ostryetum carpinifoliae (green spiderplot and NDVI profile in Figure 5c), which has the lowest NDVI values in the autumn-winter period and high NDVI values in spring-summer period, is allocated in the negative part of FPCA1 and positive of FPCA2. The NDVI profile shows a great NDVI seasonality. The 5th and 6th bi-weeks (the beginning of March) are the beginning of the season where a rapid change of the slope in the fitted NDVI values occurs, quickly reaching their maximum (0.82) at the 13th bi-week (end of June). Then, the values gradually begin to decrease until the end of the year.
The Cephalanthero longifoliae-Quercetum ilicis (orange spiderplot and NDVI profile in Figure 5c) has NDVI values about 0.75 (ranging from 0.6 to 0.8) in the autumn-winter period and the highest values in the spring-summer period. Indeed, it is allocated in the middle part of FPCA1 (both in the negative and positive parts), and in the positive part of FPCA2. The NDVI profile shows that this plant association, on average, has less NDVI seasonality than the previous one. In the fifth and the sixth bi-weeks, the NDVI values start to increase, reaching a maximum value of 0.85 at the 15th bi-week (end of July). Then the values begin to gradually decrease up to 0.75.
The Cyclamino hederifolii-Quercetum ilicis (red spiderplot and NDVI profile in Figure 5c) has high NDVI values in the autumn-winter period and in the spring-summer period. It is allocated in the positive part of FPCA1 and FPCA2. Its NDVI profile shows a low NDVI seasonality. NDVI values slightly increase in the summer period up to the 17th bi-week (end of August), where they reach the maximum value (about 0.85); then, the NDVI values decrease slightly until the end of the year.
The Pinus halepensis and P. pinea plantation (purple spiderplot and NDVI profile in Figure 5c) has NDVI values ranging from 0.75 to 0.8. Unlike the other plant associations, it has the NDVI autumn-winter values greater than the spring and summer ones. This plantation is thus allocated in the positive part of FPCA1 and the negative part of FPCA2.

Mapping Performance and Methodological Considerations
The results suggest that the Mediterranean forest plant association has distinct remotely sensed phenological behaviors and that the main seasonal phenological variations, useful to discriminate the contiguous Mediterranean habitats [17,19,20], extracted from the NDVI Landsat 8 time series using FPCA, are efficacy predictors for mapping the forest plant associations. Thus, the supervised random forest classification, similarly to Zhu and Liu [55], revealed that the main remotely sensed phenological seasonal variations (expressed by the first four pixel-based FPCA scores) contributed to the high OA (87.5%) of the map, much more than the lithological and topographic features.
This mapping methodology (based on the FPCA scores), if compared to traditional maps produced so far (e.g., [7,24]), shows several different advantages and useful potential implications regarding the Habitats Directive. This mapping process is time saving, (it required 150 hours for 650 hectares), furthermore, it facilitates the repetition of plant association mapping (every 6 years) as established by the Habitats Directive [1] and the Italian guidelines regarding monitoring systems [27]. It is also an effective tool useful for monitoring and detecting the changes over time of forest habitats, due to climate change and/or changes in management. For example, attention should be paid to the forests of Cephalanthero-Quercetum ilicis (habitat 9340), in which the dominance of the holm oak was favored by the coppicing [24]. Today, this forest is no longer managed but left to its free evolution (passive management) and the natural inter-specific competition could favor deciduous species such  The pattern distribution of the plant associations is really close to the last phytosociological map made in the traditional way. Cephalanthero longifoliae-Quercetum ilicis and Asparago acutifolii-Ostryetum carpinifoliae are distributed mainly on the north-facing slopes, exposed to the cold winter winds, on the "scaglia rossa" formation, landslide deposits, and slope deposits. The Cyclamino hederifolii-Quercetum ilicis and Pinus sp. plantations, on the contrary, are spread on the sunny south and west-facing slopes, protected from the cold winter winds on "scaglia rossa" formation and high topographic position index. The Cephalanthero longifoliae-Quercetum ilicis also occurs, in linear forms, along the impluvium line of the south-facing slopes ( Figure 6).
The whole mapping process took about 150 hours: 110 for field analysis (90 for the first step and 20 for the second), and 40 for the other laboratory analyses (visual interpretation of images, Landsat 8 images collection and pre-processing, DEM analysis, FPCA, and RF supervised classification).

Mapping Performance and Methodological Considerations
The results suggest that the Mediterranean forest plant association has distinct remotely sensed phenological behaviors and that the main seasonal phenological variations, useful to discriminate the contiguous Mediterranean habitats [17,19,20], extracted from the NDVI Landsat 8 time series using FPCA, are efficacy predictors for mapping the forest plant associations. Thus, the supervised random forest classification, similarly to Zhu and Liu [55], revealed that the main remotely sensed phenological seasonal variations (expressed by the first four pixel-based FPCA scores) contributed to the high OA (87.5%) of the map, much more than the lithological and topographic features.
This mapping methodology (based on the FPCA scores), if compared to traditional maps produced so far (e.g., [7,24]), shows several different advantages and useful potential implications regarding the Habitats Directive. This mapping process is time saving, (it required 150 hours for 650 hectares), furthermore, it facilitates the repetition of plant association mapping (every 6 years) as established by the Habitats Directive [1] and the Italian guidelines regarding monitoring systems [27]. It is also an effective tool useful for monitoring and detecting the changes over time of forest habitats, due to climate change and/or changes in management. For example, attention should be paid to the forests of Cephalanthero-Quercetum ilicis (habitat 9340), in which the dominance of the holm oak was favored by the coppicing [24]. Today, this forest is no longer managed but left to its free evolution (passive management) and the natural inter-specific competition could favor deciduous species such as Ostrya carpinifolia and Acer obtusatum in the overstory layer at the expense of Quercus ilex, with the consequent loss of the habitat 9340. Floristic variations at the plot level, as well the variations in the remotely sensed NDVI curves, could be important warning signals for changes and loss of habitats. However, this mapping provides important information to identify areas where it is possible to increase the presence of the habitat. For example, the reforested coniferous areas do not show the capacity to renew themselves. Indeed, in some sectors, in the understory layer, there is a spontaneous growth of species such as Quercus ilex, Juniperus oxycedrus, and Pistacia lentiscus which are indicators of the Cyclamino hederifolii-Quercetum ilicis association (habitat 9340-see Table A1). In these sectors, the phenological response of forested areas tends to be similar to those of Cyclamino hederifolii-Quercetum ilicis (see Figure 5). Therefore, the current mapping would allow us to highlight these areas where appropriate management measures could favor the conversion towards forests of Cyclamino hederifolii-Quercetum ilicis, with a consequent gain of the habitat 9340. A further advantage, compared to traditional mapping, is that this method provides the estimate of uncertainty, never provided in previous maps (e.g., [7,24]). The uncertainty of results can be used to improve the mapping process (e.g., to reinforce the number of permanent plots) and is important for the end users to evaluate the distribution of uncertainty over the adopted classes, suggesting how habitat classes could be misclassified with other ones (e.g., Cephalanthero longifoliae-Quercetum ilicis - Table 3). Besides the statistical data highlighted in the confusion matrix (Table 3), it is important to evaluate the quality and relevance of the errors. For example, in the Habitats Directive's perspective, the misclassifications between the Cephalanthero longifoliae-Quercetum ilicis and the Cyclamino hederifolii-Quercetum ilicis could be "acceptable" since both belong to the same habitat (9340). On the contrary, the misclassification of one of the two plant associations mentioned above with Asparago acutifolii-Ostryetum carpinifoliae or with Pinus sp. plantation should be considered to be an error. However, the high user accuracy we obtained (with an average of 90%, all plant associations have values above 80%) makes the map reliable for multiple applications in the context of the Habitats Directive (i.e., establishing conservation measures or the assessment of Plans and Projects significantly affecting Natura 2000 sites as stated in Art. 6).
The approach here proposed, mainly based on FPCA, is a valid tool to mitigate one of the main problems affecting Remote Sensing and GIS for Habitat Quality Monitoring. The key problem is the difficult comparability and direct relation between sensor-based and field monitoring data [56]. The FPCA scores, in fact, used in subsequent analyzes (e.g., correlation) have allowed us to test the relationships between the main remotely sensed phenological gradients with the topographic and floristic features and to express these relationships, similarly to the well-known (to phytosociologists and ecologists) PCA ordering technique, in a reduced phenological space ( Figure 5). This reduced (phenological) space is therefore an effective graphical tool which summarizes the relationships between phytosociological information detected in the field, the environmental information derived from a DEM, and the main remotely sensed phenological gradients, and facilitates the comparison and ecological interpretations of different perspectives of the same habitat. FPCA1 is the autumn-winter NDVI phenological gradient that, increasingly, puts in order the Asparago acutifolii-Ostryetum caprinifoliae, Cephalanthero longifoliae-Quercetum ilicis and Cyclamino hederifolii-Quercetum ilicis, and is positively correlated with the aridity and heat gradient due to lithological and topographic features. FPCA2 is the spring-summer NDVI phenological gradient that discriminates the coniferous needleleaf reforested areas from the other broad-leaved plant associations ( Figure 5).
Even the temporal NDVI profiles are an effective tool, which directly expresses the relationship between the remote-sensed phenology and the vegetation classified in the field, as similarly demonstrated in [18,57] for physiognomic types. The temporal NDVI curves of Figure 5 suggest that the NDVI time profiles can be a functional and diagnostic signature of the Mediterranean forests, not only at the physiognomic level but also at the plant association level. The comparison and the ecological interpretation of remotely sensed and field monitoring data could benefit of the plant association NDVI curves empowered by the FPCA phenological plot ( Figure 5).
Considering the increasing availability of remotely sensed data with a high spatial and temporal resolution (e.g., Sentinel 2) [56], we believe that FPCA could promote the integration between classical phytosociological analyses and remotely sensed phenological data. For example, in addition to the methodology proposed here, we hypothesize that a preliminary unsupervised classification (e.g., k-means clustering) of the FPCA scores, which could identify distinctly homogeneous areas for phenological characteristics (called pheno-cluster in Bajocco et al., [57,58]) will contribute to define efficient phytosociological sampling strategies. Furthermore, FPCA scores and NDVI profiles could be functional vegetation attributes useful for validating new (syntaxonomic) vegetation classifications [59].

Some Considerations of the Limits of This Study and Future Works
Although this mapping process is accurate and promising, with positive implications for the Habitats Directive, it should be considered that the study area is limited to 650 hectares, even if with a wide variety from the phytosociological point of view. We are extending the analysis to other central Apennine areas and the preliminary results confirm the capability of the above-discussed method to map five forest categories, two shrublands, and a grassland at the level of plant associations. Furthermore, it should be mentioned that the study excluded fragmented, residual, and linear forests, because they are not suited to the spatial resolution of Landsat-8 images (30 m), but are relevant to be mapped because, in some cases, e.g., if dominated by Quercus pubescens s.l., they are attributable to the habitat 91AA* [5]. The Sentinel-2 images, with a spatial resolution of 10 meters, were not used because of the limited availability of the second level images (available from May 2017 not covering the reference period 2013-2018). Considering the higher temporal and spatial resolution of this set of images, the proposed methodology should be evaluated to also understand the differences over time. Finally, the phenological response was anomalous in some cases, such as in the areas of the northern steep cliff, where winter images are fully shaded [60] (despite the topographic normalization carried out in the pre-processing activities of images). These anomalies are visible, for example, in Figure 5 (the elongated branches of the spiderplot are the pixels investigated in the field in the full shaded situations). However, problems related to topography and low sun elevation angle, can be avoided by using other types of images (e.g., PlanetScope by Planet). Finally, the harmonic regression (that shows neglectable inter-annual variation-see Figure 4a) and the FPCA (that not detected outliers in addition to those of the shaded pixels due to geo-morphology-see Figure 5c) confirmed the absence of meaningful land use changes (e.g., deforestation, fire . . . ). Then, the phenological responses of plant associations have not been "distorted" or "confused" by possible land use changes. However, it could be interesting and useful in the future to analyze (at the pixel level) at the same time both inter-annual and intra-annual variations by applying the FPCA.

Conclusions
The problem of mapping forest plant associations and habitats has a relevant impact on the decision taken by policy makers to monitor and preserve our environments. The generation of habitat-vegetation maps is today mainly performed by photo-interpretation, which is a time-and cost-consuming method. This work presented a new methodological framework based on FPCA to map the forest plant associations. In this work, NDVI time series from Landsat-8 were used. The obtained results show that it is possible to identify the main intra-annual seasonal variations by using Landsat-8 images and that the Mediterranean forest plant associations have different phenological behaviors. The supervised mapping of the plant associations achieved a global accuracy of 87.5%. The main seasonal (NDVI) phenological variations (FPCA components), extracted by FPCAs from the NDVI time series, contributed to the global accuracy of the map much more than the topographic and lithological variables. This aspect opens-up a new way of creating more detailed (temporal and spatial) maps. The FPCA was important to map the plant associations considering the phenology, making explicit reference to the seasonality of the vegetation without giving up the chronological order of the data. It should also be considered that the FPCA scores, which embody the chronological order of changes, can be used in subsequent analysis (e.g., ordination, clustering, correlation) enabling the comparison and interpretation of results from different points of view of the same habitat. The FPCA scores have an enormous potentiality to develop integrated analysis between the traditional field plot phytosociological and vegetational data and the remotely sensed data. An unsupervised classification of the FPCA scores could be useful to steer the phytosociological sampling, especially in partially explored and investigated areas. The NDVI time profiles and the FPCA scores could be important attributes useful in the validation processes of the vegetation classification. The NDVI annual average profile and the seasonal variations identified through FPCA showed to be very effective functional descriptors of plant associations. They could be considered to be real Plant Functional Types (PFTs) complementary to species-based approaches in plant community ecology, syntaxonomy, and biogeography. Indeed, they act as a sensible and pragmatic tool to show, monitor and preview environmental changes [61]. The obtained results demonstrated that this kind of approach: (i) is cost saving, effective, accurate and could be run periodically to output maps that reflect the "current" status of vegetation and habitats also reducing costs and risks of delaying the reporting requested by the Habitats Directive (92/43/EEC) [19]; (ii) the FPCA is a suitable and elegant solution to simultaneously express the relationships between remotely sensed and ecological field data. The approach could be applied to other vegetation typologies also supporting policy makers with updated maps to monitor and preserve biodiversity with a reduced effort in terms of required resources.
Future works will be steered to assess other vegetational indexes derived from data acquired with satellite platforms such as Sentinel-2, as well the comparison with other machine learning algorithms other than RF. Moreover, it could be interesting to analyze at the pixel level both inter-annual and intra-annual variations by applying the FPCA. For these analyzes, however, it would be necessary to collect a greater number of scenes by integrating Landsat 7, Landsat 8, and Sentinel 2.
Author Contributions: Conceptualization S.P.; Reference Data Collection and Management S.C., S.P. and G.Q.; Data analysis and programming, S.P. and G.Q.; Methodology, S.P., S.C. and A.M.; Writing S.P., S.C. and A.M.; editing, S.P. and A.M. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.    (e,f) Asparago acutifolii-Ostryetum carpinifoliae. (g,h) Pinus halepensis and P. pinea plantation. Pictures (a,c,e) allow distinguishing of Cyclamino hederifolii-Quercetum ilicis from Cephalanthero longifoliae-Quercetum ilicis and Asparago acutifolii-Ostryetum carpinifoliae. In (a) dark green, dominant, corresponds to Quercus ilex while flowering ash (Fraxinus ornus), the only deciduous species of this community, appears white because of the abundance of flowers. In (c) the abundance of flowering ash (white) and other deciduous species (light green) such as Acer obtusatum, Ostrya carpinifolia, Sorbus sp. alternating with Quercus ilex (dark green), is very evident. In (e) holm oak is absent, indeed only deciduous tall trees occur such as Ostrya carpinifolia, Acer obtusatum, and Fraxinus ornus; (g,h) are useful to recognize the occurrence of Pinus sp. plantations. These latter pictures clearly show the unique features of Pinus community of a vegetative grows in late summer. In (h) the occurrence of light green is evident, while in (g) it completely lacks.