Recursive Feature Elimination and Random Forest Classification of Meadows and Dry 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 field surveys. In this study we assess the potentiality of Recursive Feature Elimination (RFE) in combination with Random Forest (RF) classification 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 different Spectral Indexes (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 classification 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 classification 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 ecosystem services for our societies [4]. These habitats are important due to the preservation of natural values occurring in the agricultural landscape of lowland river floodplains, including fauna and flora species strictly connected with periodic floods and with extensively used agricultural to be computed from ALS and HS data respectively, for the classification of Natura 2000 alluvial meadows 6440, lowland hay meadows 6510 and xeric and calcareous grasslands 6120 in the selected 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) occurring in such valley systems. They are showed in Fig. 2.1, and are located respectively: • in the northeast of the country in the Biebrza and Narew rivers valleys, which we will refer to as "BN"; • in the central-eastern part of Poland, in the middle section of the Bug river valley, which we will refer to as "BG1"; • in the central part of Poland, in the lower section of the Bug river valley, which we will refer to as "BG2". The BN study area is of approximately 16 km 2 , covering the estuary of the Biebrza river and its floodplain, within the boundaries of a protected Natura 2000 area. On the floodplain there is a varied fluvial microrelief, visible on the Digital Terrain Models (DTM) computed from the ALS flights ( Fig.  2.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) [76]. 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).
The study areas in the Bug river valley represent two sections of the river (BG1, BG2) with a total area of about 100 km 2 , typical for the middle and lower course of the river, located in the Natura 2000 area called Ostoja Nadbużańska. Both sections of the valley are characterized by a highly diverse fluvial microrelief ( Fig. 2.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 [77]. 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 Fig.2.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. They arise in places where flooding or inundation occurs periodically [5]. Thus, they have the character of communities adapted to changing moisture conditions, usually found on fertile, alluvial soils in river valleys. They show high variability in species composition [78,79]. In plant communities, in addition to the characteristic and distinctive plant species, there may occur species that require higher moisture content, as well as those that are associated with dry habitats. The characteristic species of Cnidium meadows includes: Cnidium dubium, Allium angulosum, Scutellaria hastifolia, Viola stagnina, and Gratiola officinalis.
In river valleys, alluvial meadows are usually found in shallow depressions of the floodplain, where water does not stagnate for a long time, and thus where there is no accumulation of organic material, while alluvial mineral material is deposited. In such morphological systems, Cnidium dubium and Viola stagnina often occur in sward. There are also numerous sedges in the species composition. Alluvial meadows also develop on the slopes of small sandy or sandy-silty hills formed on the floodplain-former river bars, as well as on slightly elevated flat fragments of the plain built of silty material. Allium angulosum is a common species in such morphological systems. Species of lowland hay meadows and even xeric sandy grasslands are common among co-dominant species.
In Poland, the habitat occurs mainly in the valleys of large and medium rivers, ie: Vistula, Bug, Biebrza, Narew, Warta, Odra and others. Its exact spatial distribution in the country is not fully known [80]. The habitat usually occurs in small patches and due to the specifics of the species composition and agricultural use in some periods is difficult to identify it straight in the field. The main current threats include no periodic flooding, too intensive agricultural use or abandonment of landuse [80,81]. Habitat 6440 occurs both in the selected Bug, Biebrza and Narew valleys. A characteristic feature of plant communities of the analyzed section of the Biebrza and Narew valley (BN) is the widespread occurrence of alluvial meadows associated with flat-bottomed, extensive and shallow depressions of the floodplain, in which alluvial sediments are stratified with organic material. Among the characteristic species in the sward is Cnidium dubium and Viola stagnina (Fig. 2.3). They are associated with both local longitudinal depressions that accompany old river bars, but also occur on slopes of elevations and on extensive flat surfaces. The community is characterized by high variability in terms of species composition as well as physiognomy, and also shows a clear geographical variation. In the analyzed sections of the Bug valley (BG1 and BG2), habitat 6440 occurs less frequently.

Habitat 6510: Lowland hay meadows.
Lowland hay meadows are semi-natural communities where hay is most often obtained twice. They originated on fertile mineral soils (with a high share of silt or clays), in various morphological conditions and due to long-term not intensive agricultural use as mesic hay meadows [5,82]. They are composed of plant communities which often arise in river valleys usually outside the floodplain, on upper terraces where floods and inundation do not occur [83]. In addition, they also occur outside 7 of 36 river valleys, mainly in the mountainous areas and highlands. Their location in Poland is well known in Natura 2000 areas, but poorly recognized outside protected areas. Lowland mesic meadows are one of the richest species meadow habitats, very valuable particularly due to invertebrates, but also at high risk of transformation into arable land [82]. The community is characterized by high variability in terms of species composition (Galium mollugo, Campanula patula, Knautia arvensis, Rumex thyrsiflorus) as well as physiognomy [82]; to a small extent the community shows geographical variability.  Xeric sand calcareous grasslands (also called dry 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 occur in different landscape types [84]. Plant cover is loose and in addition to grass species contains numerous herbaceous plants, as well as moss-lichen layer typical for well-developed sandy dry grassland communities [5,85]. Plant communities of dry grasslands in Poland show great diversity and dynamics resulting from soil and geographical conditions that affect the species composition and sward density.
Habitat 6120 found in Poland in river valleys is mostly associated with sandy deposits forming extensive flat elevations both on the floodplain and on the upper terrace levels. They are beyond the reach of the flood and inundations, as occupy the highest parts of the landforms. The dry grasslands are characterized by a variable density of vegetation cover, numerous patches of exposed sandy deposits, a large share of mosses and lichens. The density of vegetation cover and its species composition reveal variability depending on the landuse, phase of plant community development, soil type and geographical location.
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 on the floodplain. Dry grasslands in the Bug valley are definitely more widespread and more typically developed.

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 beginning of September. The RS flight campaigns were conducted during the best possible weather conditions, using a RS platform with the following sensors installed: • Hyperspectral scanners (Hyspex: VNIR 0.4-0.9 µm & SWIR 0.9-2.5 µm): 470 spectral bands with 1 m spatial resolution; • Airborne Laser Scanner (Riegl Lite Mapper LMS-Q680i): point cloud data acquired with 7 points/m 2 ; • Medium format RGB camera (50Mpix) with 0.1 m spatial resolution. For more technical details about the RS platform, the reader is referred to [86]. Field spectroradiometric 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) [86][87][88]. The surveys were synchronized as much as possible with the acquisition of aerial photographs and therefore also carried out three times during the growing seasons. 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 due to 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 meters 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 in 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 GISdatabase. 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 ( Fig.2.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 is strongly influencing the quality of classification results [89,90], especially in floristically rich and heterogeneous non-forest plant communities. Moreover, due to 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 [91]. 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 [88]. This algorithm was created to efficiently visualize high-dimensional data in order to understand underlying spectral relations between groups of points. Each multidimensional polygon can be plotted in a two-dimensional space, as showed in Fig. 3.1. Points which are close to each other share similar remotely sensed (spectral) features, whilst 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 the correctly assignment of polygons into homogenous groups of the investigated habitats [88]. Incorrect points, wrongly assigned to another habitat class, were detected and removed where necessary.
The final numbers of reference polygons, obtained after the quality assessment and used for classification, are shown in table 3.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. Fig. 3.1 presents the t-SNE plotted after quality assessment, for all areas and acquisitions analyzed, depicting 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 BG1related 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.

RS data pre-processing
The HS images were radiometrically and geometrically corrected with the PARGE software [92] and atmospherically corrected using ATCOR4 model [93], using as a verification the ASD field spectral measurements collected at several locations. The result of each acquisition was a mosaicked orthophoto at 1m resolution, with 430 spectral reflectance bands (SR) in the range 450-2500nm ( Fig.  3.2). Several Spectral Indexes (SI) were calculated using the ENVI 5.3 software. The reader is referred to the ENVI's user guide [94] 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 [72]. From the ALS acquisition, 93 LiDAR-based indexes were computed after point-cloud pre-processing: vegetation structure layers, extracted from the OPALS software [95] and topographic indexes extracted from the SAGA software [96]. All produced rasters were re-projected with final spatial resolution of 1 m. The different RSbased products generated in this phase were then used as an input in the different classification experiments using the RFE-RF system, as depicted in Fig. 3.2 and explained in details in Section 3.5.

Feature selection and the Recursive Feature Elimination-Random Forest (RFE-RF) classification system
Dimensionality reduction techniques allow to handle the multi-dimensionality nature of hyperspectral data by removing the redundancy of noisy data without losing the original information content and in some cases increasing final classification accuracy [74]. Many studies can be found in literature for applying such techniques to HS data and to the simultaneous exploitation of HS and LiDAR, also integrated with the use of Machine Learning (ML), in particular Random Forest (RF) and/or Support Vector Machines (SVM) [62,63]. The Minimum Noise Fraction (MNF) [75] is a wellestablished linear transformation technique that consists of two separate PCA rotations that transform the noisy hyperspectral data cube into output channel images with steadily increasing noise levels. In general, the first few bands contain already more than 85% of the inherent spectral information with very limited noise, however the exact number of components to be selected is variable and unknown. In this work, 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 [86]. One of the major drawbacks of MNF is the lack of ''explainability'' of the resulting transformed data from a physical point of view, since the original spectral information is not pre-served [59] and therefore it could be difficult to handle and interpret.
On the other hand, 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 comprises 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.
The Recursive Feature Elimination combined with Random Forest (RFE-RF) classification has been reported as a promising technique to effectively handling the fusion of diverse data sources, such as HS and LiDAR data, at the same time generating unbiased and stable classification results in different application fields [64][65][66][67][68]. One of the objectives of this paper was therefore 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 ALS data for our habitats mapping purpose. A Vegetation Classification Studio (VCS) software was therefore developed [97], based on RF classification and RFE dimensionality reduction technique. The VCS software reads simple textlanguage 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, several procedures performing multiple cycles of classification can be easily defined and automated, bringing more confidence to the results, because based on numerous experiments, not just on a few as before. Moreover, the performing efficiency in terms of processing time needed to generate classification maps, after the preparation of the input data and reference data (described in Section 3.2 and 3.3), is much higher than other non-automated systems.
When the RFE option is used in the VCS software, an RFE report is produced in the results ( Fig.  3.3), which is employed for the definition of the optimal feature selection to be chosen as the best-touse 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. Fig. 3.3 (a) illustrates the different classification results sorted by classification run sequences (on the x-axis). The blue line (fc) shows the steady decrease in number of features (from 430 to 0 in this case) used for classification, while the green and red lines represent the produced Kappa and fuzzy Kappa respectively for each run sequence. If the same results are sorted by the Kappa values ( Fig. 3.3 (b)), 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. The RFE report provides also the features IDs used for each classification run in order to recover and produce the final list of best-to-use features. By selecting and using only the original input bands, without any kind of data-transformation, the method is a useful technique to identify which spectral channels or additional input features (i.e.: LiDAR products and Spectral Indexes) are relevant and necessary to distinguish the investigated classification problem of this paper.

Experiments design
The different experiments described below and resumed in table 3.2 were designed with the aim of addressing the objectives set out for this paper.
3.5.1 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, described in Section 3.4. The 430 SR bands were compared to the 30 MNF components, by running two different experiments (see also Fig. 3.2). In the first experiment (Exp.1a, table 3.2), 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 classification of the investigated habitats. In the second experiment (Exp.1b, table 3.2), 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 but randomly on a 50/50 basis: ten different random sample selections were manually realized and ten different shapefiles produced for the classification (table 3.2). In this way, instead of using just one training/validation selection, ten different RF classification runs were realized and a better statistical 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 find 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 Indexes and the 65 SI computed in Section 3.2 ( Fig. 3.2), were all added as input features in the Exp.2. The same ten manually generated training/validation shapefiles used in Exp.1a/1b were also adopted for this experiment (table 3.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 find a certain number of commonly selected features which might represent the bestto-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 twofold: 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 not performed on all the features generated by HS and ALS pre-processing (Fig. 3.2), but only using the selection of best-touse features resulting from Exp.2. Obviously, in this case the RFE technique was not used during RF classification. 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 3.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, Exp.3 was divided in two classification experiments: • 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
The feature selected by the RFE technique when classifying the 430 SR bands (Exp. 1a) are plotted in Fig. 4.1, for each study area and for each acquisition campaign (different colors). The best-to-use spectral channels were selected for each of the ten classification runs based on the approach described in Section 3.4 and results of this selection are plotted in Fig. 4.1. The results show that the mostly used spectral channels belong to the following spectral ranges: • 400-800nm of the visible spectral range (mainly red and blue) • 1050-1100nm of the Near-Infrared For each of the ten RF classification runs (scenario 1-10, as explained in table 3.2), the Kappa accuracies obtained by the optimal features selection after applying the RFE (Exp.1a) are plotted in Fig. 4.2 against the Kappa accuracies obtained by classifying the 30 MNF components alone (Exp.1b). The ƿ-value resulting from the "independent samples t-test" are also plotted. A ƿ-value less than 0.05 (typically ≤ 0.05) indicates that there is a statistically significant difference between the two results, while on the contrary for a ƿ-value higher than 0.05, there is not a statistically significant difference [98].
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 (Section 3.5.2 and Fig. 3.1). The Kappa accuracies of Exp.2 are summarized and compared to the Kappa accuracies of the Exp.1b in Fig. 4.3. The ƿ-values resulting from the "independent samples t-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 ƿ-value is always lower than the significance level of 0.05. Moreover, the final classification accuracy is now higher than a tolerance level of K=0.65 (red dotted line in the figure) for all Exp. 2 results. On the contrary, if only the 30 MNF components are to be used, in several cases the Kappa is below 0.65. As explained in Section 3.5.2, also in this experiment ten classification runs with different manual splitting of training/validation samples were realized (table 3.2). For each run an optimal selection of features was extracted (as explained in Fig. 3.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 Fig. 4.4. The results show that there are features being selected 100% of times for all study areas (Fig.4.4a In order to perform the Exp.3a/3b, aimed at verifying the classification performances when using only the optimal features selected as a result of Exp.2, we decided to screen the features selected at least more than 50% of times in two different ways: area-by area (for Exp.3a) and considering all study areas together (for Exp.3b), based on Fig. 4.4. Therefore, a list of features selected and used for Exp.3a/3b was produced and is presented in table 4.1, 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 indexes were always selected for both 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 Fig. 4.5 the Kappa values of Exp.3a/3b for all areas and acquisitions are plotted and compared with respect to Exp.2. The ƿ-value resulting from this comparison are also plotted. The Kappa values are very similar for most classification results and in most study areas analyzed. Apart some exceptions in BG1 area, the ƿ-value is mostly higher than the significance level of 0.05, meaning that in general there are no significant differences between Exp.2, 3a or 3b. For the summer acquisition of the BG1 area, Exp. 2    Considering that only in one case (BG1-Summer), the Exp.3b produced significantly lower accuracies, while in all other cases no reduction in classification was remarked, the final classification maps were produced with the Exp.3b selected features (table 4.1) and are plotted in Fig. 4.6.

The importance of different HS spectral 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 high 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 most used spectral channels are: visible spectral range (400-800 nm), NIR (1050-1100 nm) and four different ranges in the SWIR part of the electromagnetic spectrum (1250-1400 nm, 1650-1800 nm, 1950-2050 nm and 2250-2400 nm). These figures show the importance of HS information, especially in the SWIR channel, for mapping the investigated habitats. Both of them in fact are 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 distinguish such types of grasslands habitats.
When comparing the Kappa accuracies, we found that MNF results outperformed the RFE technique in all cases analyzed, in a statistically significant way. In other words, the 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 are rather low and around 0.50-0.55 in most cases analyzed (except spring and summer acquisitions of BN study area). The major drawback of MNF transformation on the other hand, 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 Indexes 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). Therefore, a total of 188 features were computed and used for this RFE-RF classification exercise: 30 MNF components, 93 LiDAR-products and 65 SI were classified together, using ten selections of training/validation sets, the same employed for the previous test.
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, as compared to using the HS data alone. In fact, all Kappa accuracies have scored a remarkable and statistically significant increase. Although in some cases the classification results are still a bit low, for example the summer acquisition of BG2 area (K=0.656), in most study areas, the Kappa values reached a quite satisfiable result (K>0.65), as for example the spring acquisition of BN area (K=0.773), the summer acquisition of BG1 (K=0.712) or the spring acquisition of BG2 (K=0.7). 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 and at the same time very subtle spectral mutual differences both between each other's and to other non-habitats vegetation patches (background class 9999) that are found in the studied river valleys. Worth noticing is that in both study areas investigated, the autumn acquisition is the one generating the lowest Kappa, while the highest have been alternatively the spring or summer acquisitions.

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 selected for all classified datasets. Moreover, also the DSM standard deviation of the unit weight (OPALS_DSM_sigma0) was selected for all datasets by the RFE-RF.
The Topographic Position Index (TPI) is calculated as the difference between real elevation of each cell and the averaged elevation in a selected window of neighborhood cells. This index 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.) [102]; moreover, species distribution models showed significant relationships to TPI in [118]. 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 [101]. Unlike MRVBF, the MRRTF index does not have a clear link to landform processes but it has been found to be a useful adjunct to MRVBF in landform classification [101].
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 [83], 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.
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 [103] and have been reported to be important indicators for describing floods [119], 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 [83,84], 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 the Biebrza and Bug river valleys.

Spectral Indexes of importance.
The Spectral Indexes (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 then the topographic-based features described above; the selection of each index is more scattered and linked to specific areas and probably also 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 indexes were selected by RFE-RF during Exp.2 (table 4.1). In this section we will discuss the most relevant and common ones.
The Canopy Nitrogen indexes 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, by using the 1510nm and 1680nm spectral bands both in the SWIR part of the electromagnetic spectrum, and the sensitiveness of leaves to nitrogen concentration, as well as the overall foliage biomass of the canopy [111]. 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 indexes 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 [104]. The CRI1, computed with 510nm and 550nm 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,82], whilst dry grasslands (habitat 6120) have been reported to develop on sandy soils [84]. 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 [117] 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 indexes is also selected in the different acquisitions and study areas analyzed (table 4.1). 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 [109,110]. Among these, three types of indexes 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 Indexes 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 [113,114], and the Structure Insensitive Pigment Index (SIPI), developed to assess canopy stress [115], is selected several times for BG2 area.
The Dry or Senescent Carbon indexes 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, due to absorptions in the 2000-2200nm range to cellulose [105]. Similarly, the Normalized Difference Lignin Index (NDLI), using 1680nm and 1754nm spectral bands, provides an estimation of the relative amounts of lignin contained in vegetation canopies [111]. 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 [84]: this might explain why it was selected. Finally, two Geology Indexes, Iron Oxide Ratio (IO) and WorldView New Iron Index (WVII) are selected for BG1 and BG2 areas respectively.
The scattered selection of different Spectral Indexes 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 [78,79], 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.

5.3
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 by repeating the classification procedure using only the best-to-use features without the RFE technique and increasing the classification runs from 10 to 50. 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 using those layers that were selected for more than 50% of times among the three study areas. From this feature selection exercise, it was revealed that Spectral Indexes 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, whilst SI, depicting different vegetation characteristics and conditions (e.g.: vegetation stress, chlorophyll/nitrogen concentrations, moisture content, etc..), have by nature a higher spatial local 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 different training/validation selections, instead of only 10.
Due to 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 (table 4.1): • LiDAR (SAGA) products: DurI, MRRTF, MRVBF, TPI, MCA, TWI; • LiDAR (OPALS) products: DSM_Sigma0, DTM_Sigma0; • SI: CR1, CM, NDNI, WVWI; • MNF: 1-11. Therefore, we can claim that the optimal feature selection generalization attempt was a reliable choice, because 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, hence reducing computing time and efficiency. By generalizing the feature selection, however, several other area-dependent SI are not used for classification (table 4.1). This might cause, in some exceptional cases, some slighter reductions in Kappa values, as reported for example in BG1-summer acquisition (from K=0.712 to K=0.697).

Conclusion
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 required 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 producing a 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. 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 [101]. 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 [103], have also been selected as best-to-use features; • The common feature selection showed also the importance of using HS data with high spectral resolution covering a broad part of the spectral range so to compute specific Spectral Indexes (SI), which also 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 and 6510 and confirmed the HS channels selection resulting in the first part of our work; • The remarkable results obtained in this study wouldn't have been possible without a good quality field-based sampling activity, which was conducted extensively by a dedicated team of botanist. A high number of samples (between 1000 and 1500 polygons) were collected in the field for each area, which required a great time-effort.
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 regions of the electromagnetic spectrum and to include LiDAR data when acquiring HS data, due to the important role in achieving significantly higher results .