Recursive Feature Elimination and Random Forest Classiﬁcation of Natura 2000 Grasslands in Lowland River Valleys of Poland Based on Airborne Hyperspectral and LiDAR Data Fusion

: The use of hyperspectral (HS) and LiDAR acquisitions has a great potential to enhance mapping and monitoring practices of endangered grasslands habitats, beyond conventional botanical ﬁeld surveys. In this study we assess the potentiality of recursive feature elimination (RFE) in combination with random forest (RF) classiﬁcation in extracting the main HS and LiDAR features needed to map selected Natura 2000 grasslands along Polish lowland river valleys, in particular alluvial meadows 6440, lowland hay meadows 6510, and xeric and calcareous grasslands 6120. We developed an automated RFE-RF system capable to combine the potentials of both techniques and applied it to multiple acquisitions. Several LiDAR-based products and di ﬀ erent spectral indices (SI) were computed and used as input in the system, with the aim of shedding light on the best-to-use features. Results showed a remarkable increase in classiﬁcation accuracy when LiDAR and SI products are added to the HS dataset, strengthening in particular the importance of employing LiDAR in combination with HS. Using only the 24 optimal features selection generalized over the three study areas, strongly linked to the highly heterogeneous characteristics of the habitats and landscapes investigated, it was possible to achieve rather high classiﬁcation results (K around 0.7–0.77 and habitats F1 accuracy around 0.8–0.85), indicating that the selected Natura 2000 meadows and dry grasslands habitats can be automatically mapped by airborne HS and LiDAR data. Similar approaches might be considered for future monitoring activities in the context of habitats protection and conservation.


Introduction
Protecting and monitoring natural habitats is essential for the mitigation of biodiversity decline, caused by the negative effects of increased human activity [1] and the natural adaptation to climate change. Back in 1992, the European Union established the Habitats Directive [2], a cornerstone European nature conservation policy, developing the wide Natura 2000 ecological network of protected areas with the aim of preserving the most valuable sites in the European landscape. Within the Natura 2000 habitats, grasslands are one of the most biodiverse habitats in Europe [3], delivering essential The recent enhancements in the spectral and spatial resolution of sensors have enabled the development of new opportunities for a wider application of hyperspectral imaging (HS) in environmental monitoring [32]. Characterized by both high spatial and spectral resolution, they typically record the reflected object's light across hundreds of narrow spectral bands, within the visible, near infrared (NIR), mid infrared (MIR), and short-wave infrared (SWIR) parts of the electromagnetic spectrum [33,34], resulting in a 3-dimensional hyperspectral data-cube [35]. More and more studies have recently focused on exploiting this wealth of spectral information for natural habitats mapping [19,36]. For example, few works have been found in the current literature focusing on the exploitation of HS potentials for heathland habitats mapping and conservation, using different techniques such as, spectral unmixing [37], machine learning classification [38], or object-based image analysis [39,40].
In recent years, it has also been proven and recognized that light detection and ranging (LiDAR) systems, also referred to as airborne laser scanning (ALS), are a good alternative for vegetation mapping and characterization [41][42][43][44][45]. Some studies have shown how the combination of ALS with HS data can enhance the accuracy of habitats mapping [46][47][48][49], emphasizing the importance of simultaneously exploiting both the spectral and topographical information when targeting the automatic identification of such types of natural habitats [50]. However, the effective fusion of the two complementary data sources is a challenging task. The high dimensionality of the data and the limited availability of ground-truth samples, especially when it comes to natural habitats mapping and characterization, pose major challenges in data handling and processing [51,52]. The first one, also referred to as the curse of dimensionality [53], can be addressed by using dimensionality reduction methods [54], supposed to remove redundant data without losing original and relevant information. Among these techniques, the feature-selection based approaches, whose objective is to find a subset of the relevant data from the original dataset, are reported to be simpler, easier to implement, and very efficient at the same time [55][56][57]. The future-selection methods are also reported to be effective methods in handling the fusion of diverse data sources, such as HS and LiDAR data, as they allow to consider all available information simultaneously in a decision-making system in which the different sensor data can be selected [58]. In the last years, several investigators have carried out research in this specific field, combining the feature-selection capabilities with machine learning (ML) classification [52,[59][60][61], for a broad range of applications [62][63][64]. Among these techniques, it has been highlighted by several authors that recursive feature elimination (RFE) combined with random forest (RF) classification could provide unbiased and stable results in different application fields [65][66][67][68]. However, to our knowledge, this method was not investigated yet for mapping meadows and dry grasslands.
The main objective of this study is thus to assess the performances of RFE-RF for integrating airborne hyperspectral and LiDAR data in the context of mapping selected Natura 2000 habitats found along Polish river valleys; in particular, alluvial meadows (6440), lowland hay meadows (6510), and xeric and calcareous grasslands (6120). With this work, we aim at identifying the main features required to automatically classify these habitats and at the same time emphasizing the added value of the fusion between such complementary and diverse datasets.
Airborne flight campaigns with HS and LiDAR sensors were conducted simultaneously to botanical field-surveys three times in the growing seasons (Spring/Summer/Autumn), along three significative river valleys in Poland, with a good distribution of lowland meadows and dry grasslands habitats, as described in Section 3.1. With the aim of combining the RFE method with RF classification, an RFE-RF automated classification system was developed and is described in Section 3.4. Three different experiments were designed in line with the objectives of this paper, as described in Section 3.5. In the first one, we aimed at identifying the relevant spectral features for meadows and dry grasslands mapping, by applying the RFE on the hyperspectral dataset alone and comparing results with a standard, well-known dimensionality reduction technique, the minimum noise fraction transformation (MNF) [69][70][71][72]. The most effective method for hyperspectral dimensionality reduction was then retained for the second experiment, in which several spectral-and topographical-based metrics were computed and added to the RFE-RF classification system. Using the RFE technique, we derived for each case study, a list of optimal features selection which might represent the best-to-use features for the investigated classification problem. In the third experiment, as a sort of verification test, only the best-to-use selected features were used to run the final classification maps. In this way, we could attempt to generalize a best-to-use features selection to be calculated from ALS and HS data respectively, for the classification of the selected Natura 2000 grasslands (6440, 6510 and 6120) in the analyzed river valleys.

Study Areas
The research study areas are located in chosen lowland river valleys in Poland (part of the North European lowland), on sections characterized by the presence of well-developed semi-natural, open Natura 2000 natural habitats (in particular habitats 6440, 6510 and 6120).
The BN study area is approximately 16 km 2 (Figure 1), covering the estuary of the Biebrza and Narew rivers and its floodplain, in the northeast of the country. On the floodplain there is a varied fluvial microrelief, visible on the digital terrain models (DTM) computed from the ALS flights ( Figure 2). The characteristic elements are: flat-bottomed depressions, narrow longitudinal depressions of overgrowing old river beds, flat elevations, small longitudinal elevations of old river bars, natural levees, as well as remnants of higher terrace levels (with traces of aeolian process) [73]. The floodplain is dominated by mineral and organic-mineral deposits, and also other organic materials, such as mud. Higher terrace levels are made of sandy and sandy-silty formations. The development of natural habitats in both valleys is influenced by the system of embankments and dikes extending both along the Biebrza (limiting the floodplain from the west) and Narew rivers (the embankment located south of the riverbed).
Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 35 derived for each case study, a list of optimal features selection which might represent the best-to-use features for the investigated classification problem. In the third experiment, as a sort of verification test, only the best-to-use selected features were used to run the final classification maps. In this way, we could attempt to generalize a best-to-use features selection to be calculated from ALS and HS data respectively, for the classification of the selected Natura 2000 grasslands (6440, 6510 and 6120) in the analyzed river valleys.

Study Areas
The research study areas are located in chosen lowland river valleys in Poland (part of the North European lowland), on sections characterized by the presence of well-developed semi-natural, open Natura 2000 natural habitats (in particular habitats 6440, 6510 and 6120).
The BN study area is approximately 16 km 2 (Figure 1), covering the estuary of the Biebrza and Narew rivers and its floodplain, in the northeast of the country. On the floodplain there is a varied fluvial microrelief, visible on the digital terrain models (DTM) computed from the ALS flights ( Figure  2). The characteristic elements are: flat-bottomed depressions, narrow longitudinal depressions of overgrowing old river beds, flat elevations, small longitudinal elevations of old river bars, natural levees, as well as remnants of higher terrace levels (with traces of aeolian process) [73]. The floodplain is dominated by mineral and organic-mineral deposits, and also other organic materials, such as mud. Higher terrace levels are made of sandy and sandy-silty formations. The development of natural habitats in both valleys is influenced by the system of embankments and dikes extending both along the Biebrza (limiting the floodplain from the west) and Narew rivers (the embankment located south of the riverbed).   (BG2, BN, BG1). The spatial distribution of reference polygons is also displayed (yellow dots).
The selected study areas have in common a highly variable topographical micro-relief (as described above and highlighted in Figure 2), a typical condition for the development of the investigated Natura 2000 habitats. In such highly variable topographical context, the extraction of LiDAR-based metrics could be very meaningful and is expected to bring a significant contribution to the quality of the classification results.  The study areas in the Bug river valley represent two sections of the river (BG1, BG2, Figure 1), typical for the middle and lower course, with a total area of about 100 km 2 . Both sections of the valley are characterized by a highly diverse fluvial microrelief ( Figure 2), with fragments of flat flood-basin and upper terraces, longitudinal, narrow depressions of oxbow lakes, natural levees, extensive flat-bottomed depressions on the floodplain, as well as longitudinal and narrow raised berms or crests above the floodplain surface [74]. The soil cover is definitely dominated by alluvial soils developed on sandy and sandy-silty formations, locally showing traces of aeolian processes.
The selected study areas have in common a highly variable topographical micro-relief (as described above and highlighted in Figure 2), a typical condition for the development of the investigated Natura 2000 habitats. In such highly variable topographical context, the extraction of LiDAR-based metrics could be very meaningful and is expected to bring a significant contribution to the quality of the classification results. Alluvial meadows of river valleys of the Cnidion dubii are semi-natural, fertile communities, usually mown twice, in Poland associated with the floodplains of main rivers. They arise in places where flooding or inundation occurs periodically [5]. Thus, they are adapted to changing moisture conditions, and are usually spotted in shallow depressions, on fertile, soils, on the slopes of alluvial forms or on slightly elevated flat parts of the plain built of sandy or silty material. They show high variability in species composition [75,76]. In addition to the characteristic and distinctive plant species (e.g., Cnidium dubium, Allium angulosum, Scutellaria hastifolia, Viola stagnina, and Gratiola officinalis), it may occur species that require higher moisture content, as well as those that are typically associated with dry habitats. The main current threats include no periodic flooding, too intensive agricultural use, or abandonment of landuse [77,78]. Habitat occurs in both chosen case studies for this research. In the Biebrza and Narew valleys (BN), habitat 6440 is commonly associated with flat-bottomed, extensive and shallow depressions ( Figure 3). In the sections of the Bug valley (BG1 and BG2), it occurs less frequently and it is related mostly to narrow depressions of former oxbows or longitudinal and narrow raised berms or crests.

Habitat 6510: Lowland Hay Meadows
Lowland hay meadows are semi-natural communities where hay is most often obtained twice a year. They originated on fertile mineral soils in various morphological conditions and because of long-term not intensive agricultural use [5,79]. The habitat is composed of plant communities which arise in river valleys on upper terraces-beyond the extends of floods and inundation [80]. The lowland hay meadows are considered as one of the richest in plant species, hosting numerous fauna (mainly invertebrates) and being highly variable in terms of physiognomy [79], but at the same time at high risk of transformation into arable fields [79]. The most often observed characteristic plant species of the habitat are: Galium mollugo, Campanula patula, Knautia arvensis, Rumex thyrsiflorus. Because of the lack of typically shaped upper terraces, as well as the light sandy alluvial deposits, patches of habitat in the analyzed section of the BN valley are rare and scattered. They are found on the few flat elevations of the floodplain, where they form intermediate plant communities, with a large share of species typical of dry grassland habitats ( Figure 3). In the Bug valley, lowland hay meadows occur in BG1 at upper terraces, which are not transformed into intensively used meadows or arable fields.

Habitat 6120: Xeric and Calcareous Grasslands
Xeric sand calcareous grasslands are semi-natural plant communities that have originated as a result of extensive grazing of grasslands on permeable, relatively fertile, and sandy soils. Thus, they occur in different landscape types [81] but very often in Poland are spotted in river valleys on extensive flat elevations outside the flood extends. The plant cover is loose, and besides grass species (e.g.,: Koeleria glauca, Festuca polesica, Corynephorus canescens), it contains numerous herbaceous plants (e.g.,: Silene otites, Silene tatarica, Astragalus arenarius, Kochia laniflora), as well as moss-lichen layer [5,82]. Plant communities of dry grasslands in Poland show great diversity and dynamics resulting from landuse, soil and geographical conditions that affect the species composition, and sward density. Currently the most serious threats are related to land abandonment, lack of regular grazing, and overgrowing by shrubs. In the analyzed section of the BN valley, habitat 6120 occurs on small areas (up to 100 m 2 ), and is associated with distinct sandy elevations. Dry grasslands in the Bug valley are definitely more widespread and more typically developed ( Figure 3). An extensive campaign of botanical surveys was conducted by teams of specialized botanists, according to a uniform procedure worked out within the HabitARS project (Habitats Airborne Remote Sensing) [83][84][85]. Botanical research did not last more than seven days from the date of airborne flight campaigns. The procedure for obtaining botanical reference material involved several stages. At the stage of the preliminary reconnaissance, the area was analyzed because of the diversity of natural habitats (based on available data) and to the spatial distribution of habitat patches. Once recognized, GPS points were recorded using a GPS receiver of 0.5 m accuracy. Reference polygons of 3-m radius were established at locations typical and representative for the particular habitat in the

Airborne Data Acquisition and Botanical Field Measurements
Airborne flight campaigns and botanical field data were acquired simultaneously three times in the growing season (Spring/Summer/Autumn). The campaign dates were chosen taking into account the phenological period and agricultural practices. Thus, the optimal terms for obtaining aerial photographs and conducting field survey were planned at the end of May, the middle of July and the For more technical details about the RS platform, the reader is referred to [83]. Field spectro-radiometric measurements were also taken in the field with ASD FieldSpec 4, at selected bright and dark locations to be used in the later stage for atmospheric correction (Section 3.2).
An extensive campaign of botanical surveys was conducted by teams of specialized botanists, according to a uniform procedure worked out within the HabitARS project (Habitats Airborne Remote Sensing) [83][84][85]. Botanical research did not last more than seven days from the date of airborne flight campaigns. The procedure for obtaining botanical reference material involved several stages. At the stage of the preliminary reconnaissance, the area was analyzed because of the diversity of natural habitats (based on available data) and to the spatial distribution of habitat patches. Once recognized, GPS points were recorded using a GPS receiver of 0.5 m accuracy. Reference polygons of 3-m radius were established at locations typical and representative for the particular habitat in the research area, homogenous due to vegetation type and structure of plant community, without physical disturbances caused by animals (e.g., wild boars), agricultural machinery etc. Reference polygons did not encompass trees and shrubs higher than 1.5 m, and were also located at a distance from high elements (e.g., trees, escarpments) in order to avoid shadows. For each reference polygon, a standard description of vegetation and applied agricultural practices was performed. The botanical surveys were also supplemented by photographic documentation and in a later stage were organized into a GIS-database. The same procedure was repeated for as many patches of habitats 6120, 6440, and 6510 as possible, in order to have a good spatial distribution of samples ( Figure 1). Other GPS points were also collected in vegetation patches not considered to be a habitat or representing other land-cover classes and were assigned to the background class (code 9999).

Reference Botanical Data Quality Assessment
The quality of field reference dataset is a crucial aspect that must be considered carefully because it strongly influences the quality of classification results [86,87], especially in highly heterogeneous plant communities. Moreover, because of the non-simultaneity of some field measurements with respect to airborne acquisitions and to the fact that these natural areas are affected by mowing, an appropriate selection and quality assessment procedure was implemented before using the reference dataset for automatic classification.
The first step of this procedure was the visual recognition with the simultaneously acquired RGB-aerial photographs, in order to eliminate defective/erroneous reference polygons. Polygons were removed from the reference set when on the given RGB mosaic they were: mowed, in shadow, inundated, mechanically damaged (by animals or farming), or had other problems. In the second step, a cleaning of the dataset was performed using the t-distributed stochastic neighbor-embedding (t-SNE) algorithm [88]. In a recent study, it has been demonstrated that the t-SNE algorithm is able to increase the final classification accuracy of 6% in Kappa coefficient [85]. This algorithm was created to efficiently visualize high-dimensional data in order to understand the underlying spectral relations between groups of points. Each multidimensional polygon can be plotted in a two-dimensional space, as showed in Figure 4. Points which are close to each other share similar remotely sensed (spectral) features, while different samples are represented by distant points. In this way, a multi-dimensional dataset can be interpreted in an easier way, revealing its global structure, such as the presence of clusters and their spectral characteristics.
In this study, t-SNE algorithm was used to evaluate if reference polygons were correctly assigned to the belonging habitat class [85]. A series of iterative visualizations, with perplexity hyper-parameter varying between 5 and 130, were performed. PCA-based initialization was used to provide stable global layouts. Figure 4 presents the t-SNE plots using a perplexity value of 30 and after removal of detected incorrect points, for all areas and acquisitions analyzed. It depicts the real spectral characteristics of the investigated habitats. Habitats 6120 and 6440 are clearly clustered and it seems that they hold quite different spectral characteristics among both BN and BG2 areas. Instead, habitat 6510 seems slightly more mixed with other classes, being in-between habitats 6120 and 6440 in BN-related plots. However, if we compare habitat 6510 only to habitat 6120 in BG1-related plots, it seems they have more clear differences. Over all plots, the background class (code 9999) is the one having the highest heterogenous spectral characteristics.
The final numbers of reference polygons, obtained after the quality assessment and used for classification, are shown in Table 1 for each habitat and study area analyzed. The higher number of background polygons (class 9999) is justified by the large acquisition areas (such as BG1) and therefore by the much higher spatial distribution of heterogeneous non-habitat classes.

RS Data Pre-Processing
In Figure 5 a scheme of the pre-processing steps is illustrated, together with the experiments design (Section 3.5) and the corresponding addressed objectives in the paper. The HS images were radiometrically and geometrically corrected with the PARGE software [89] and atmospherically corrected using ATCOR4 model [90], using as a verification the ASD field spectral measurements collected at several locations. The result of each acquisition was a mosaicked orthophoto at 1-m resolution, with 430 spectral reflectance bands (SR) in the range 450-2500 nm. Several spectral indices (SI) were calculated using the ENVI 5.3 software. The reader is referred to the ENVI's user guide [91] for a detailed listing of all 65 SI computed for this study. The minimum noise fraction transformation (MNF), a well-known technique for hyperspectral dimensionality reduction and denoising, was also computed for each image [69]. Based on eigenvalues and after testing with different number of bands, we decided to retain the first 30 MNF components, also as suggested in other studies [83]. From the ALS acquisition, 93 LiDAR-based indices were computed after point-cloud pre-processing. They consisted of vegetation structure layers, extracted from the OPALS software [92], and of topographic indices extracted from the SAGA software [93]. All produced rasters were re-projected with final spatial resolution of 1 m. The different RS-based products generated in this phase were then used as an input in the different classification experiments using the RFE-RF system, as depicted in Figure 5 and explained in detail in Section 3.5.

Recursive Feature Elimination-Random Forest (RFE-RF) Classification System
The RFE is a feature selection method that aims at estimating which features are most helpful to discriminate the classes of interest. It is able to eliminate any features that are not useful in this task, in order to obtain the input feature-set having the lowest possible number of layers, at the same time without reducing the final classification accuracy. The algorithm relies on variable importance assessment, which is calculated internally by RF classifiers and requires performing multiple rounds of classification. Each round consists of learning a new RF classification model, assessing its accuracy based on cross-validation, analyzing the feature of importance metrics for every feature used, and modifying the feature-set that would be used for the successive round of the procedure. The first classification round is performed using all the available features. Then, the weakest performing ones are detected using variable of importance metric, estimated by the model during learning. One (or more) of the weakest features are then eliminated from the feature-set, and the next round of the procedure is performed. By doing so, RFE attempts as well to eliminate dependencies and collinearity that may exist in the input features.
One of the objectives of this paper is to develop an RFE-RF system, capable of combining the potentiality of RF classification and RFE feature selection, at the same time automatizing the analysis of HS and LiDAR data for our habitats mapping purpose. A vegetation classification studio (VCS) software was therefore developed [94], based on RF classification and RFE dimensionality reduction technique. The VCS software reads simple text-language commands defined by users, which allow the automatization of the whole RF classification procedure, such as splitting reference data into training/validation sets, model training and validation, feature selection with RFE, accuracy assessment and plotting classification maps results, where needed. In this way, multiple classification cycles can be easily defined and automated, bringing more confidence to the results, since based on numerous experiments and not just a few. Moreover, the performing efficiency in terms of processing time needed to generate classification maps, after the preparation of the input and reference data (described in Sections 3.2 and 3.3), is much higher than other non-automated systems.
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 35 radiometrically and geometrically corrected with the PARGE software [89] and atmospherically corrected using ATCOR4 model [90], using as a verification the ASD field spectral measurements collected at several locations. The result of each acquisition was a mosaicked orthophoto at 1-m resolution, with 430 spectral reflectance bands (SR) in the range 450-2500 nm. Several spectral indices (SI) were calculated using the ENVI 5.3 software. The reader is referred to the ENVI's user guide [91] for a detailed listing of all 65 SI computed for this study. The minimum noise fraction transformation (MNF), a well-known technique for hyperspectral dimensionality reduction and denoising, was also computed for each image [69]. Based on eigenvalues and after testing with different number of bands, we decided to retain the first 30 MNF components, also as suggested in other studies [83]. From the ALS acquisition, 93 LiDAR-based indices were computed after point-cloud pre-processing. They consisted of vegetation structure layers, extracted from the OPALS software [92], and of topographic indices extracted from the SAGA software [93]. All produced rasters were re-projected with final spatial resolution of 1 m. The different RS-based products generated in this phase were then used as an input in the different classification experiments using the RFE-RF system, as depicted in Figure 5 and explained in detail in Section 3.5.

Recursive Feature Elimination-Random Forest (RFE-RF) Classification System
The RFE is a feature selection method that aims at estimating which features are most helpful to discriminate the classes of interest. It is able to eliminate any features that are not useful in this task, in order to obtain the input feature-set having the lowest possible number of layers, at the same time without reducing the final classification accuracy. The algorithm relies on variable importance assessment, which is calculated internally by RF classifiers and requires performing multiple rounds of classification. Each round consists of learning a new RF classification model, assessing its accuracy When the RFE option is used in the VCS software, an RFE report is produced in the results ( Figure 6). This is employed in the next step for the definition of the optimal feature selection to be chosen as the best-to-use feature configuration out of the multiple run sequences, generally being the one having the best compromise in terms of higher Kappa and lower number of features. Figure 6a illustrates the different classification results sorted by classification run sequences (on the x-axis). If the same results are sorted by the Kappa values (Figure 6b), it is possible to locate the point where a major increase in the number of features (green line, fc) will not cause a major change in the Kappa results. In this example, run 100 produced a Kappa of 0.631 using 23 features, while a similar Kappa value of 0.632 was obtained using 430 features. Therefore, run 100 can be retained as the optimal feature selection because by using a much lower number of features a very similar Kappa accuracy can be achieved.
Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 35 dimensionality reduction technique producing the highest classification results was then retained and used for the next experiment (Exp.2).

Exp.2: LiDAR Products and Spectral Indices
After the dimensionality reduction technique had been selected, the following objective was to investigate which other input features could increase the classification accuracy achieved in Exp.1a/1b by using the HS dataset alone. For this purpose, the 93 LiDAR indices and the 65 SI computed in Section 3.2 were all added as input features in the Exp.2 (see also Figure 5). The same ten manually generated training/validation shapefiles used in Exp.1a/1b were also adopted for this experiment ( Table 2). The RFE technique (Section 3.4) was also employed during RF classification, using the same approach as described in Exp.1a.
The aim this time was to analyze the lists of optimal feature selections among the ten different runs and then finding a certain number of commonly selected features which might represent the best-to-use HS and ALS features to map the selected meadows and dry grasslands habitats.

Experiments Design
The different experiments described below and resumed in Table 2 were designed with the aim of addressing the objectives set out for this paper.

Exp. 1: Spectral Features of Importance and Dimensionality Reduction
One of the objectives of this paper is to understand which spectral features are more important for automatic classification of selected meadows and dry grasslands habitats and which dimensionality reduction technique is most effective. For this reason, the RFE technique was compared to the well-established data transformation technique MNF. The 430 SR bands were compared to the 30 MNF components, by running two different experiments ( Figure 5). In the first experiment (Exp. 1a), only the 430 SR bands were analyzed. The RFE technique was used together with RF classification (RFE-RF), in order to identify which spectral bands are mostly used for the classification of the investigated habitats (see also Table 2). In the second experiment (Exp. 1b), only the 30 MNF components were used for RF classification and without the RFE option. In both cases (Exp. 1a/1b), selection of samples was done manually on a 50/50 basis: ten different random sample selections were manually realized and ten different shapefiles produced for the classification (Table 2). In this way, instead of using just one training/validation selection, ten different RF classification runs were realized and a better statistical Remote Sens. 2020, 12, 1842 13 of 33 distribution of both Kappa accuracy and optimal spectral features selection will be possible. The aim was to analyze the lists of optimal features selections among the ten different runs and then finding a certain number of commonly selected HS channels, which might represent the spectral features of importance for the habitats 6120, 6440, and 6510. The dimensionality reduction technique producing the highest classification results was then retained and used for the next experiment (Exp. 2). After the dimensionality reduction technique had been selected, the following objective was to investigate which other input features could increase the classification accuracy achieved in Exp. 1a/1b by using the HS dataset alone. For this purpose, the 93 LiDAR indices and the 65 SI computed in Section 3.2 were all added as input features in the Exp. 2 (see also Figure 5). The same ten manually generated training/validation shapefiles used in Exp. 1a/1b were also adopted for this experiment ( Table 2). The RFE technique (Section 3.4) was also employed during RF classification, using the same approach as described in Exp. 1a.
The aim this time was to analyze the lists of optimal feature selections among the ten different runs and then finding a certain number of commonly selected features which might represent the best-to-use HS and ALS features to map the selected meadows and dry grasslands habitats.

Exp. 3: Feature Selection Validation Attempt
The objective of this last experiment was two-fold: verifying that the optimal features were properly selected in Exp. 2 and attempting a generalization of the selected optimal features among the different study areas analyzed. The RF classification was therefore run using only the selection of best-to-use features resulting from Exp. 2, not on the entire list of features generated by HS and ALS pre-processing ( Figure 5). Obviously, in this case the RFE technique was not used during RF classification (Table 2). Moreover, the splitting between training/validation polygons was done automatically by the software on a 50/50 random basis and repeated for 50 different runs, so to have 50 different RF classification results ( Table 2). In this way, a much higher statistical distribution of Kappa results, based on the optimal features selection, can be investigated and discussed as compared to Exp. 2.
In order to decide which optimal feature selection would be better, single features were used in Exp. 3 if selected at least 50% of times, based on two different approaches: • Exp. 3a: the selection of optimal features was done separately for each area, meaning that a different best-to-use features selection was used for each case study, as a result of Exp. 2; • Exp. 3b: the selection of optimal features was the same for all study areas, meaning that a unique best-to-use features selection was extrapolated from all study areas at the same time, as a result of Exp. 2.
By performing these two separate experiments we can test: (1) How the RF classification performs when using a very limited number of features; (2) how a different selection of best-to-use features could affect the classification results; and (3) attempt to generalize a best-to-use features selection for the classification of the investigated habitats over the three selected study areas.

Results of Exp. 1
The best-to-use spectral channels selected by the RFE technique (based on the approach described in Section 3.4) when classifying the 430 SR bands (Exp. 1a) are plotted in Figure 7, for each study area and for each acquisition campaign (different colors). The results show that the mostly used spectral channels belong to the following spectral ranges: The ƿ-value resulting from the "Wilcoxon signed-rank test" [95] are also plotted. A ƿ-value less than 0.01 (typically ≤ 0.01) indicates that there is a statistically significant difference between the two results, while on the contrary for a ƿ-value higher than 0.01, there is not a statistically significant difference [96].   For each of the ten RF classification runs (scenario 1-10, as explained in Table 2), the Kappa accuracies obtained by the optimal features selection after applying the RFE (Exp. 1a) are plotted in Figure 8 against the Kappa accuracies obtained by classifying the 30 MNF components alone (Exp. 1b). The p-value resulting from the "Wilcoxon signed-rank test" [95] are also plotted. A p-value less than 0.01 (typically ≤ 0.01) indicates that there is a statistically significant difference between the two results, while on the contrary for a p-value higher than 0.01, there is not a statistically significant difference [96]. If we compare all results, we can observe that in all cases the MNF results outperform the RFE results, with statistically significant difference. For example, for the BN spring acquisition, the MNF results (Exp. 1b) have a mean Kappa of 0.72, while the RFE results (Exp. 1a) have a mean Kappa of 0.656. As another example, for the BG2 summer acquisition, the MNF results have a Kappa of 0.59, while the RFE results have a Kappa of 0.521. For this reason, the 30 MNF components were retained for further processing and Exp. 2 was run by adding all others input features (93 LiDAR-products and 65 SI) to the 30 MNF components, as described in Section 3.5.2 and Figure 4.

Results of Exp.2
The Kappa accuracies of Exp.2 are summarized and compared to the Kappa accuracies of the Exp.1b in Figure 9. The ƿ-values resulting from the "Wilcoxon signed-rank test" are also plotted. As

Results of Exp. 2
The Kappa accuracies of Exp. 2 are summarized and compared to the Kappa accuracies of the Exp. 1b in Figure 9. The p-values resulting from the "Wilcoxon signed-rank test" are also plotted. As we can see, the accuracy is increased when adding LiDAR and SI features (Exp. 2) in a statistically significant way for all cases analyzed: the p-values are always lower than the significance level of 0.01. Moreover, the final classification accuracy is now higher than a tolerance level of K = 0.65 (red dotted line in Figure 9) for all Exp. 2 results. For each run of the ten runs in Exp. 2, an optimal selection of features was extracted (as explained in Figure 6, Section 3.4) and a frequency distribution of most selected features was computed. This can be done either considering all the study areas together, or area by area, as illustrated in Figure 10.  The area-related frequency distributions (Figure 10b) reveal some variability in the mostly used features. For example, for the BN study area, also the MNF components [02,04,06] were selected for all scenarios (100% of times). Other features selected more than 75% of times are: • LiDAR products: SAGA_TWI, SAGA Duration of Insolation (SAGA_DurI), OPALS mean

Results of Exp.3
The results of Figure 10 were used as a supportive tool to decide which features are to be used for the next Exp.3a (area-by-are feature selection) and Exp.3b (selection based on all study areas together). Therefore, a list of features was produced and is presented in Table 3, with corresponding explanation of each indicator. For Exp.3b a fixed number of 24 features are to be used, while for Exp.3a between 24 and 28 features are to be used, depending on the study area. The main difference consists in the different SI being selected from one study area to another. In fact, only the NDNI and The area-related frequency distributions (Figure 10b) reveal some variability in the mostly used features. For example, for the BN study area, also the MNF components [02,04,06] were selected for all scenarios (100% of times). Other features selected more than 75% of times are: •

Results of Exp. 3
The results of Figure 10 were used as a supportive tool to decide which features are to be used for the next Exp. 3a (area-by-are feature selection) and Exp. 3b (selection based on all study areas together). Therefore, a list of features was produced and is presented in Table 3, with corresponding explanation of each indicator. For Exp. 3b a fixed number of 24 features are to be used, while for Exp. 3a between 24 and 28 features are to be used, depending on the study area. The main difference consists in the different SI being selected from one study area to another. In fact, only the NDNI and WVWI indices were always selected for the three study areas. Similarly, from the LiDAR products, MRRTF, MRVBF, TPI, MCA, DSM_Sigma0 are all being used for all study areas. From the MNF components, only some of them were used for all study areas, strengthening the spectral heterogeneity of the different habitats' classes within the study areas analyzed.
In Figure 11, the Kappa values of Exp. 3a/3b for all areas and acquisitions are plotted and compared with respect to Exp. 2. The p-values resulting from this comparison are also plotted. The Kappa values are very similar for most classification results and in most study areas analyzed. Apart from some exceptions in BG1 area, the p-values are mostly higher than the significance level of 0.01, meaning that in general terms there are no significant differences between Exp. 2, 3a, or 3b. For the summer acquisition of the BG1 area, Exp. 2 and Exp. 3a have no significant difference, while Exp. 3b resulted in a statistically significant lower accuracy (K = 0.697 instead of K = 0.712). The opposite occurs for the autumn acquisition of the BG1 area, where Exp. 2 and Exp. 3b have no significant difference, while Exp. 3a produced a significantly lower accuracy (K = 0.661 instead of K = 0.675). On the other hand, for the spring acquisition of BG1 area, Exp. 3a and Exp. 3b generated a slightly higher mean Kappa, with a statistically significant difference in respect to Exp. 2. This is the only case in which Exp.  Considering that only in one case (BG1-Summer) the Exp. 3b produced significantly lower accuracies, while in all other cases no significantly reduction in classification was remarked, the final classification maps were produced with the Exp. 3b selected features (Table 3) and are plotted in Figure 12. Table 3. List of features selected more than 50% of times from Exp. 2 and used in Exp. 3a/3b, with corresponding explanations and references. Colors in the right columns mean that the feature was used in the corresponding experiment (3a or 3b). The different colors are associated to different feature groups, using same coloring as in Figure 10.

The Importance of Different Hyperspectral Channels and Dimensionality Reduction Techniques for Mapping Meadows and Dry Grasslands Habitats in the Selected River Valleys
The Exp.1 was designed with the aim of identifying the most relevant hyperspectral channels, required to classify the selected meadows and dry grasslands habitats (code 6440, 650, and 6120) and to choose the most efficient dimensionality reduction technique between RFE and MNF.
Because the spectral heterogeneity of each habitat class was pretty high due to strong variability in habitats' species composition, we decided to realize ten different training/validation selections and therefore run ten different RF classifications. In this way, results can be considered more reliable and less influenced by an individual sample selection. When applying the RFE technique on the 430 SR bands, a similar pattern of spectral bands was selected among the ten different classification scenarios realized for each study area and acquisition period. The obtained results emphasized the importance of HS information for mapping the investigated habitats, especially of the SWIR channel. Both habitats are in fact characterized by very different soil moisture conditions, therefore it is probable that the spectral information contained in the SWIR channel is a key source of information to enable distinguishing such types of grasslands habitats.

The Importance of Different Hyperspectral Channels and Dimensionality Reduction Techniques for Mapping Meadows and Dry Grasslands Habitats in the Selected River Valleys
The Exp. 1 was designed with the aim of identifying the most relevant hyperspectral channels, required to classify the selected meadows and dry grasslands habitats (code 6440, 650, and 6120) and to choose the most efficient dimensionality reduction technique between RFE and MNF.
Because the spectral heterogeneity of each habitat class was pretty high due to strong variability in habitats' species composition, we decided to realize ten different training/validation selections and therefore run ten different RF classifications. In this way, results can be considered more reliable and less influenced by an individual sample selection. When applying the RFE technique on the 430 SR bands, a similar pattern of spectral bands was selected among the ten different classification scenarios realized for each study area and acquisition period. The obtained results emphasized the importance of HS information for mapping the investigated habitats, especially of the SWIR channel. Both habitats are in fact characterized by very different soil moisture conditions, therefore it is probable that the spectral information contained in the SWIR channel is a key source of information to enable distinguishing such types of grasslands habitats.
When comparing the Kappa accuracies, we found that MNF transformation is a better and more effective technique for removing the inherent noise in the hyperspectral data cube without losing the relevant spectral information required to obtain a higher level of classification accuracy. On the other hand, the RFE technique, although reducing the number of spectral bands during classification, loses part of the inherent spectral information in a way that the classification performances are reduced in a statistically significant manner: Kappa values were rather low and around 0.50-0.55 in most cases analyzed. The major drawback of MNF transformation is that the transformed bands have no physical meaning and cannot provide insights on the spectral channels of importance for the classification of the selected habitats.
Despite of the fact that RFE is not able to maintain classification accuracies as high as MNF method, the results of Exp. 1 have proved its usefulness as a tool to highlight and identify which hyperspectral channels are mostly used for the investigated habitats classification. Something that is not possible to perceive with the MNF technique alone. This is an important information to shed lights on the spectral configuration requirements for monitoring these Natura 2000 habitats, for example in view of future hyperspectral data acquisition mission planning. Nevertheless, the rather low levels of Kappa accuracies obtained even when using the MNF components alone (mostly K<0.65), suggest the need of other sources of data in order to get a higher and more satisfiable level of classification accuracy for mapping the selected habitats.

LiDAR Products and Spectral Indices Selection to Enhance Habitats Classification Performances
The Exp. 2 was designed with the aim of understanding which LiDAR and SI products are to be produced for the investigated classification problem and if these extra input features might enhance the rather low classification accuracies obtained by using the HS dataset alone (mostly with K < 0.65, Exp. 1b). In line with other studies focusing on mapping similar habitat types [50], the integration of LiDAR-derived products and SI obtained from different spectral bands has proven to be a good solution for getting higher classification accuracy. In fact, in most study areas, the Kappa values have scored a statistically significant increase, reaching a quite satisfiable result (K > 0.65). This is a great achievement, considering the very challenging task of being able to automatically distinguish different types of selected Natura 2000 habitats, characterized by a high within-class heterogeneity. These natural grasslands are in fact composed by a high variability of species, described by different spectral characteristics within the same habitat class. At the same time, they can be very similar to each other and present similar spectral characteristics with respect to other non-habitats vegetation patches (background class 9999) that are found in the studied river valleys. The combination of all these aspects, makes their automatic identification in a natural landscape a very challenging task.

LiDAR-Based Features of Importance
The frequency distributions of the selected features of importance revealed the great influence of the LiDAR-based products, especially computed from the SAGA software. In fact, the three morphological indicators (TPI, MRRTF, and MCA) proved to be fundamental inputs to reach a good classification level, because they are selected 100% of times. Among them, the multiresolution index of the ridge top flatness (MRRTF) is a topographic index designed to identify high-flat areas at different scales. It complements the multiresolution index of valley bottom flatness (MRVBF) that instead is designed to identify areas of deposited material in flat valley bottoms [99].
The dependence of habitats types occurrence to different topographical characteristics of the landscape, justifies the selection of these indicators as the most important selected features. In fact, the lowland hay meadows (habitat 6510) found in Biebrza and Bug river valleys, are characterized by plant communities which often arise outside the floodplain, in higher terraces where floods or inundation do not occur [80], therefore with probably higher MRRTF values. On the other hand, alluvial meadows (habitat 6440) are usually found in shallow depressions of the floodplain [5], that might be depicted by low MRVBF values. Finally, dry grasslands (habitat 6120), are associated with sandy deposits, forming extensive flat elevations both on the floodplain and on higher terrace levels, described by spatially different values of MRRTF and MRVBF. Similarly, the topographic position index (TPI), has been reported to be useful to classify the landscape into slope position (i.e., ridge top, valley bottom, etc.,) and landform category (i.e., gentle valleys, plains, open slopes, etc.) [100]; moreover, also species distribution models showed significant relationships to TPI in [116].
The modified catchment area (MCA), together with the topographic wetness index (TWI), have been developed for predicting the spatial distribution of soil moisture contents in hilly terrains [101] and have been reported to be important indicators for describing floods [117], which are a very recurrent and important phenomenon in the Biebrza and Bug river valleys investigated in this work (Section 2.1). In fact, habitats 6510, 6440, and 6120 are all strictly connected to periodic floods. Plant communities of lowland hay meadows (habitat 6510) and dry grasslands (habitat 6120) are beyond the reach of the flood and inundations [80,81], while alluvial meadows (habitat 6440) communities are able to adapt to changing moisture contents and are found in shallow depressions of the floodplain, characterized by temporary flooding [5].
It is remarkable seeing that the most important features being selected by the RFE-RF system are exactly those topographic and wetness indicators (TPI, MRRTF, MRVBF, MCA, and TWI) that better portray the characteristics of the Natura 2000 habitats found in our investigated study areas, along both the Biebrza and Bug river valleys.

Spectral Indices of Importance
The spectral indices (SI) computed under the ENVI software, using different HS channels, have also proven to be among the most important features to be used for classification. However, their importance and impact are less evident than the topographic-based features described above. Figure 10 clearly showed that the selection of each index is more scattered and linked to specific areas and therefore probably to local habitats conditions. In fact, by their nature, SI are meant to picture very specific vegetation conditions, that might be difficult to generalize for large and highly heterogenous plant communities, as they might be more suited to describe local-specific phenomena or only some specific habitats. Different categories of indices were selected by RFE-RF during Exp. 2 (Table 3). In this section we will discuss the most relevant and common ones.
The canopy nitrogen indices measure the nitrogen concentration in vegetation, which is an important component of chlorophyll and which is present in high concentrations in quickly growing vegetation. Among this group, the normalized difference nitrogen index (NDNI) was designed with the aim of measuring the relative estimate of nitrogen in vegetation canopies, the sensitiveness of leaves to nitrogen concentration as well as the overall foliage biomass of the canopy [109]. It is based on the 1510 nm and 1680 nm spectral bands, both in the SWIR part of the electromagnetic spectrum. The NDNI is one of the most important SI, used more than 75% times to classify our habitats among all study areas analyzed. This result is in line with the results obtained when analyzing the mostly used spectral channels of the HS dataset (Section 5.1) and highlights the importance of having the SWIR part of the electromagnetic spectrum measured.
Leaf pigments indices measure the stress-related pigments, which are present in higher concentrations in weakened vegetation. Among this group, the carotenoid reflectance index 1 (CRI1) is a measure of stressed vegetation, based on the level of carotenoid concentration, which regulate the light absorption processes in plants [102]. The CRI1, computed with 510 nm and 550 nm spectral bands, is another important feature selected by the RFE-RF system most of the times for the investigated study areas. Another important selected SI is the clay mineral ratio (CM), a band ratio working in two different regions of the SWIR spectral range, aimed at highlighting the presence of clay. Lowland hay meadows (habitat 6510) have been reported to originate in mineral soils, with high share of silt or clays [5,79], while dry grasslands (habitat 6120) have been reported to develop on sandy soils [81].
The dependence of habitats types occurrence to different soil types justifies the selection of CRI1 as an important selected feature. Moreover, the WorldView water index (WVWI), designed to highlight areas of standing water [115] and working in the NIR and blue spectral ranges, is also selected among the most important features for the selected study areas, reflecting the different moisture content of the investigated habitats [5].
A number of Narrowband Greenness indices is also selected in the different acquisitions and study areas analyzed (Table 3). They are a combination of reflectance measurements in the green, red and NIR spectral ranges aimed at measuring the overall amount and quality of photosynthetic material in vegetation, which is essential for understanding the state of vegetation [107,108]. Among these, three types of indices dedicated to the estimation of chlorophyll abundance are selected (MCARI, MCARI2, TCARI) and two types dedicated to the estimation of green leaf-area-index (TVI and MTVI).
The light use efficiency indices are highly linked to the carbon uptake, somewhat related to fractional absorption of photosynthetically active radiation (fAPAR). Among this group, the photochemical reflectance index (PRI), mainly selected for BG1 area, is used in studies of vegetation health and agricultural crops productivity and stress [111,112]. Besides, the structure insensitive pigment index (SIPI), developed to assess canopy stress [113], is selected several times for BG2 area.
The dry or senescent carbon indices provide an estimate of the amount of carbon in dry states of lignin and cellulose. Among this group, the cellulose absorption index (CAI) highlights dried plant materials, because of absorptions in the 2000-2200 nm range to cellulose [103]. Similarly, the normalized difference lignin index (NDLI), using 1680 nm and 1754 nm spectral bands, provides an estimation of the relative amounts of lignin contained in vegetation canopies [109]. Both NDLI and CAI are selected only for BG2 study area. This index might operate as an indicator of dryness, a characteristic of habitats 6120 [81]; this might explain why it was selected. Finally, two geology indices, iron oxide ratio (IO), and WorldView new iron index (WVII) are selected for BG1 and BG2 areas respectively.
The scattered selection of different spectral indices has well pictured the spectral heterogeneity nature of the investigated habitats. In fact, habitats 6120, 6440, and 6510 are both characterized by a high variability in terms of species composition as well as physiognomy [75,76], which both contributes to the selection of different SI depending on the area and habitats analyzed. However, among the selected SI, the most important ones are computed using the SWIR channels and subsequently the NIR and visible parts of the spectrum.

Best-to-Use HS+ALS Products for Mapping Meadows and Dry Grasslands Habitats in the Selected River Valleys
The last set of experiments (Exp. 3a and Exp. 3b) was designed with the aim of verifying the appropriateness of the optimal features selection. In the first case (Exp. 3a), the optimal features were selected area by area, while in the latter case (Exp. 3b) a common feature selection was performed considering the three study areas together. From this feature selection exercise, it was revealed that spectral indices are more dependent on the individual study areas (different SI selected for different study areas), while the LiDAR products are more uniform over the three study areas. This could be explained by the fact that the topographical microrelief variability is more homogeneous among the analyzed river valleys, while SI, depicting different vegetation characteristics and conditions (e.g.,: vegetation stress, chlorophyll/nitrogen concentrations, moisture content, etc.), have by nature a higher local spatial variability.
The classification results showed very similar outcomes in most cases analyzed (both Exp. 3a and 3b), demonstrating the effectiveness of the RFE-RF system for classification and feature selection. With a much lower number of features, 24-28 instead of 188, it was possible to reach similar performances in terms of classification accuracies, at the same time reducing processing time and data handling. Moreover, reliability of results was also increased by using 50 classification runs, based on random training/validation selections, instead of only 10.
Because of the fact that almost no significant reduction of classification accuracy was recorded when generalizing the feature selection to the three areas (Exp. 3b), the final classification results were produced using the 24 features identified by this common feature selection, in particular ( Therefore, we can claim that the optimal feature selection generalization attempt was a reliable choice. Similar classification accuracies can be produced when using the same selected features for the three study areas, without the need to calculate several other products area by area. By generalizing the feature selection, however, several other area-dependent SI are not used for classification (Table 3). This might cause, in some exceptional cases, slight reductions in Kappa values, as reported for example in BG1-summer acquisition (from K = 0.712 to K = 0.697).

Considerations on the Computational Efficiency of the RFE-RF System
The RF-RFE system proved to be beneficial from the computational performance point of view. The only upfront cost is that over hundred models had to be learned in sequence to obtain the optimized set of features. However, this needs to be done only once per dataset, then the successive prediction phase can benefit from it. The prediction phase is the most demanding one, because it is applied on the entire pixels (and features) in the whole study areas. Therefore, every eliminated feature results in faster prediction time, since less data have to be read from disk, transferred over the network, and processed by the CPU. Typically, by using only 1/3 of the original number of features, we can expect roughly even 3x decrease of our prediction time. In our case, classification times were in the range of minutes for training and below half an hour for predicting the whole study area (on a modern, budget machine), so performance considerations were not very critical. But especially for larger and more complex datasets we consider implementing the RF-RFE procedure a win-win situation, providing considerable savings in computational effort without losing in terms of final classification accuracy.

Conclusions
The main objective of this study was to develop an automated system based on recursive feature elimination and random forest classification (RFE-RF), capable of complementing the potentiality of RF classification and RFE feature selection, for the mapping of selected Natura 2000 meadows and dry grasslands (habitats 6120, 6440, and 6510) along Polish river valleys by the fusion of airborne Hyperspectral and LiDAR data. For this purpose, multiple flight acquisitions were performed along Biebrza and Bug river valleys (Poland) and extensive field-botanical activities were conducted to built-up a reliable database to be used as a reference for the classification. Several experiments were performed to test the different potentialities of the RFE-RF system and to tackle several research questions. The conclusions we can draw from the results of this work are as follows: • The MNF dimensionality reduction method outperformed RFE: part of the inherent spectral information was lost during RFE feature selection and therefore the classification accuracy significantly reduced. On the other hand, with the RFE technique, it was possible to highlight the important HS channels that are necessary for mapping the investigated habitats: VIS, NIR, and SWIR are all required spectral channels; • By selecting and using only the original input bands, without any kind of data-transformation, the RFE-RF system proved to be a very efficient and useful setup for automated hyperspectral and LiDAR data processing, highlighting the added value of the fusion between these complementary and diverse datasets. It is therefore possible to use a common selection of 24 features (instead of 188) to distinguish the investigated habitats of this study and still obtain very similar classification accuracies among all considered study areas; • LiDAR-based products, depicting the variable topographic micro-reliefs of the investigated river valleys, proved to be the most selected features producing also a significant enhancement in the classification accuracy. In particular: topographic position index (TPI), multiresolution index of the ridge top flatness (MRRTF), multiresolution index of valley bottom flatness (MRVBF), modified catchment area (MCA), topographic wetness index (TWI), DSM_Sigma0 and DTM_Signma0 proved to be necessary adjuncts for mapping Natura 2000 habitats 6120, 6440 and 6510; • The meaningfulness of the selected products is strongly linked to the habitats' characteristics. It was remarkable noticing that since habitat 6510 is mostly found in higher terraces, while habitat 6440 is usually found in low depressions of the floodplain, the MRRTF and MRVBF products were retained as the most relevant classification features. In fact, MRRTF and MRVBF have been conceived for the identification of high-flat areas and flat valley bottoms respectively [99]. Likewise, since habitats 6120, 6510, and 6440 are all connected to periodic floods in different ways, the TWI and MCA products, reported in literature to be good indicators of floods [101], have also been selected as best-to-use features; • The common feature selection also showed the importance of using HS data with high spectral resolution covering a broad part of the spectral range, so to compute specific spectral indices (SI) which significantly contributed to enhancing the final classification accuracies. The high heterogeneity of the habitats is well pictured by the selection of different SI in different study areas. However, it was proved that using only CRI1, CM, NDNI, and WVWI, together with the other topographical products and MNF components, it was possible to obtain satisfiable classification accuracies over all investigated areas. All of them are strictly linked to the highly variable characteristics and conditions of the diverse habitats analyzed: vegetation stress and health (CRI1), presence of different soil types and minerals (CM), nitrogen concentrations (NDNI), and moisture content (WVWI). The selection of these specific indicators highlight also the importance of the SWIR, NIR and visible channels for mapping Natura 2000 habitats 6440, 6120, 6510 and confirms the selected HS channels in the first part of our analysis; • A great time-effort needs to be envisaged to collect a high number of field-based samples (between 1000 and 1500 polygons) in order to achieve similar classification results. This is a very important step in order to embark on a classification problem of this kind.
To conclude, the present study showed the feasibility of the RFE-RF system to map Natura 2000 meadows (habitats 6440 and 6510) and dry grasslands (habitat 6120) from the fusion of airborne HS and LiDAR data, with a great level of accuracy (max K = 0.774, min K = 0.652). For future acquisition planning aimed at monitoring these Natura 2000 habitats, it is strongly recommended to select proper optical sensors with high spectral resolution (covering the SWIR, NIR and visible spectral ranges), to include LiDAR data when acquiring HS data and to implement a feature selection system similar to the RF-RFE developed here, so to process and classify efficiently such a high dimensional dataset.