{'fraction_non_alphanumeric': 0.04050830005776137, 'fraction_numerical': 0.02574148019789548, 'mean_word_length': 5.018742442563482, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 18, 'https://': 12, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Classifying wetland‐related land cover types and habitats using fine‐scale lidar metrics derived from country‐wide Airborne Laser Scanning

Mapping 3D vegetation structure in wetlands is important for conservation and monitoring. Openly accessible country‐wide Airborne Laser Scanning (ALS) data—using light detection and ranging (lidar) technology—are increasingly becoming available and allow us to quantify 3D vegetation structures at fine resolution and across broad spatial extents. Here, we develop a new, open‐source workflow for classifying wetland‐related land cover types and habitats using fine‐scale lidar metrics derived from country‐wide ALS data. We developed a case study in the Netherlands with a workflow consisting of four routines: (1) pre‐processing of ALS data, (2) calculation of lidar metrics (i.e. 31 features representing cover, 3D shape, vertical variability, horizontal variability and height of vegetation as well as microtopography), (3) assessing feature importance of lidar metrics for classifying wetland habitats, and (4) applying a Random Forest algorithm for mapping and prediction. We used an expert‐based vegetation map for annotation and generated 100, 500 and 1000 annotation points for each class. Using a three‐level hierarchical approach, we differentiated at level 1 planar surfaces (e.g. roads and agricultural fields) from wetland vegetation with 82% mean overall accuracy, using predominantly height and horizontal variability metrics. At level 2, we classified wetland vegetation into four land cover types (forest, grassland, reedbeds, shrubs) with 71% mean overall accuracy, using lidar metrics related to vegetation height and horizontal and vertical variability. At level 3, we differentiated two types of land reed as well as water reed with 78% mean overall accuracy, using predominantly vertical variability metrics. Our results demonstrate that lidar metrics (related to vegetation height, cover, vertical and horizontal variability) derived from country‐wide ALS data can differentiate land cover types and habitats within wetlands at high resolution. Given appropriate annotation data, our workflow can be up‐scaled to a country‐wide extent to allow the comprehensive mapping and monitoring of wetlands at national scales.


Introduction
Wetlands are key habitats for animals and plants, are among the most vital ecosystems on Earth, and provide multiple ecosystem services to humans (Mitsch et al., 2015). Wetlands consist of a large variety of different land cover types and habitats, including meadows and wet grasslands formed by herbaceous vegetation (e.g. grasses, sedges and bulrush), patches of trees and shrubs (e.g. European alder, willows), and reedbeds (i.e. mono-dominant stands of the common reed Phragmites australis) (Ellenberg, 1996;Weller, 1999). Since plants in wetlands vary in life form and physical stature, they form different structural habitat components which are important for many wetland animals, including migratory and breeding birds (Weller, 1999;Strayer and Findlay, 2010). Wetlands in general, and reedbeds in particular, have been strongly modified by changes in water chemistry, climate, land-use, urban development and water level regulations (Ostendorp, 1989;Davidson, 2014). Consequently, international agreements such as the Ramsar Convention on Wetlands (https://www.ramsar.org/) urge to monitor and map wetlands for effective planning and natural resource management, and to foster national action and international cooperation for the conservation and wise use of wetlands and their resources.
Remotely sensed data provide an efficient method for mapping and monitoring biodiversity and habitat change across broad spatial extents and with fine resolution (Pettorelli et al., 2016;Guo et al., 2017;Navarro et al., 2017). Manual interpretation of stereographic aerial photographs has traditionally been used to map wetland extent (Howland, 1980). More recently, satellite-based multispectral (Shadaydeh et al., 2017) and radar (Muro et al., 2016) data, as well as hyperspectral imagery (Hirano et al., 2003), has allowed the (semi-)automatic mapping and monitoring of wetlands. Light detection and ranging (lidar) technology is increasingly used because it captures the 3D physical structure of terrestrial vegetation and habitats (Lefsky et al., 2002;Vierling et al., 2008) and provides fine-scale digital elevation models for characterizing microtopography (Brubaker et al., 2013). A common type of lidar measurement is Airborne Laser Scanning (ALS) where the sensor is mounted, for instance, on an airplane or helicopter. The technique uses the time difference between the pulse emitted from the sensor and the return signal from the ground from which x,y,z coordinates of points characterizing objects can then be calculated. Besides the coordinates of the points, the intensity of the returned pulse is recorded which depends on the backscattering properties and distribution of the objects. The resulting 3D point cloud typically varies in point density depending on the type of the scanner, acquisition parameters and whether single or multiple return or full-waveform data are recorded by the instrument. The resulting high-resolution information on the vertical and horizontal distribution of vegetation enables to quantify 3D habitat structures across large spatial extents (Coops et al., 2016;Moeslund et al., 2019). Because country-wide ALS data are now increasingly becoming available, there is potential to characterize finescale habitat structures not only in forests (where much of lidar-based research has been conducted) but also in lowstature habitats such as wetlands.
Lidar can be used to derive ecologically relevant information on 3D vegetation structure and requires the processing of the 3D point cloud into suitable metrics, for example, by statistically aggregating the 3D point cloud information within raster cells (Davies and Asner, 2014;Bakx et al., 2019). These lidar metrics can quantify vegetation height (e.g. maximum of height values within a cell), vertical vegetation variability (e.g. the standard deviation of height values within a cell) or horizontal vegetation variability (e.g. the standard deviation of height values between neighbouring cells). Moreover lidar metrics have also been operationally used in forestry to derive stand inventory parameters such as biomass or basal area (Hyypp€ a et al., 2008;Cao et al., 2014). Within wetlands, commonly used lidar metrics are those derived from Digital Surface Models (DSMs), Canopy Height Models (CHMs) or Digital Terrain Models (DTMs) (Lindsay et al., 2004;Goodale et al., 2007;Gilmore et al., 2008;Onojeghuo et al., 2010;Onojeghuo and Blackburn, 2011;Lang et al., 2013;Chasmer et al., 2016;Rapinel et al., 2018). Height metrics such as the mean and maximum percentiles of height (Alexander et al., 2015;Ene et al., 2018), canopy density measures (Ding et al., 2011;Kope c et al., 2016;Corti Meneses et al., 2017;Ene et al., 2018), the standard deviation of points within cells (Alexander et al., 2015;Kope c et al., 2016;Ene et al., 2018), or the variance of the DSM have also been used in wetland studies (Zlinszky et al., 2012;Alexander et al., 2015;Chasmer et al., 2016). Beyond these examples, ALS is not yet widely used for characterizing wetland-related land cover types and for differentiating habitats within wetlands, probably because robust and open-source workflows and ground-based annotation datasets are largely lacking.
Here, we develop a new open-source workflow that uses fine-scale lidar metrics derived from country-wide ALS data for classifying land cover types and habitats. We test this with a case study in a Dutch wetland to make a first step towards upscaling the use of country-wide ALS data. To the best of our knowledge, this is the first study that uses country-wide ALS data for classifying land cover types and habitats within wetlands in the Netherlands. We use a wall-to-wall digitized expert-based map for annotation because ground-based annotation data were not available. The openly accessible ALS dataset comes with multiple returns and a point density of 6-10 points/m 2 . We derived 31 lidar metrics to capture the 3D structure of the vegetation by quantifying aspects related to vegetation cover, 3D shape, vertical variability, horizontal variability and height, as well as fine-scale topographic variability. We then used a Random Forest algorithm to evaluate feature importance and to classify land cover types and habitats at different hierarchical levels. We expect that from the available set of lidar metrics only a few metrics are needed to classify wetland-related land cover types and habitats, and that metrics related to vegetation height and horizontal as well as vertical vegetation variability play an important role for classification.

Materials and Methods
Airborne Laser Scanning data The ALS data (AHN2) was captured during the second Dutch national ALS flight campaign (Actueel Hoogtebestand Nederland, https://www.pdok.nl/) which took place from 2007 to 2012. Our study area was captured in 2008 and 2009 in the leaf-off season (Northern hemisphere winter). The AHN2 dataset covers the whole Netherlands with a point density of 6-10 points/m 2 from multiple return data acquisition and with an overall vertical accuracy of 20 cm (https://www.ahn.nl/). The dataset is openly accessible (https://ahn.arcgisonline.nl/ahnviewer/) and has been pre-processed by "Rijkswaterstaat" (the executive agency of the Dutch Ministry of Infrastructure and Water Management). The available dataset comes with a ground classification that separates ground points from non-ground points. The dataset consists of the x,y,z coordinates of all points which are stored in LAZ format files (https://rapidlasso.com/laszip/). Other characteristics such as the number of returns, the intensity values and other flight-related parameters (e.g. scan angle, pulse reputation frequency and flight height) are not provided by the publisher of the dataset.

Study area and annotation dataset
We selected the study area Lauwersmeer in the northern part of the Netherlands in the province of Groningen because it is a relatively large (5754 ha) wetland area with diverse land cover types and habitats (Fig. 1). This area is recognized as a Ramsar site (site no. 1247) and has the status of a Nature 2000 site (site code NL9802012) because it is an important breeding site for coastal, water and reedbed birds. Since 2008, field data on vegetation and birds have been collected for monitoring the potential consequences of soil subsidence caused by gas extraction (Kleefstra et al., 2016). Based on these surveys, a digital vegetation map has been produced by an ecological research and consultancy firm (Altenburg & Wymenga, http://www.altwym.nl/). The mapping of vegetation structure in Lauwersmeer started in 2005 on the basis of vegetation height and coverage of the shrub and forest layers as well as the coverage of Phragmites australis within reedbeds. Updates were made in 2008/2009 and 2012. The final vegetation map was made in 2015 and consists of twenty-one digitized vegetation types (Fig. 1). These capture important habitats for animals (especially wetland birds). The vegetation map was derived from an expert-based interpretation of 3D stereoscopic aerial photographs using both textural and actual height information in combination with field visits (Kleefstra et al., 2016). We used the 1:5000 expertbased vegetation map as a source for our annotation data (training and testing points). The final vegetation map was made 6-7 years after the acquisition of the lidar data. The vegetation monitoring between 2007 and 2012 indicates that no substantial changes in vegetation structure have occurred within the annually visited transects (Bijkerk et al., 2013).
We used this expert-based vegetation map (Fig. 1) to create annotation data for the Random Forest algorithm. The algorithm was applied at three different levels which were hierarchically nested (Table 1). At level 1, the wetland class was separated into planar surfaces and wetland vegetation; at level 2, four vegetated wetland-related land cover types (forest, grassland, reedbed and shrubs) were separated from each other; and at level 3, the reedbed category was separated into three different types of reedbed habitats (water reed, structurally poor land reed, and structurally rich land reed). To apply the machine learning algorithm, we aggregated and simplified the original vegetation units for the classification (Table 1). We generally excluded polygons that describe vegetation units as "open" from the training step (see Fig. 1) because these consist of mixtures of vegetation (grasslands, shrubs, trees, reedbed patches) within the user-defined polygons. Furthermore, at level 1 we included low-stature vegetation (pioneer vegetation and grassland <0.3 m; Fig. 1) into 'planar surfaces'. At level 2, we included vegetation polygons described as high grassland (>0.3 m), dense forest (only broad-leaved trees), dense reedbeds and dense shrub (Table 1). At level 3, we used the reedbed habitats dense water reed (flooded) and dense land reed (non-flooded). Land reed was further separated into structurally poor land reed vs. structurally rich land reed (Table 1). Structurally rich land reed consists of 25-50% cover of the common reed Phragmites australis, along with herbs as co-dominant vegetation (associations Phragmitetea and Convolvulo-Filupenduletea) (Kleefstra et al., 2016). In contrast, structurally poor land reed and water reed consist of >50% of common reed. The reclassified habitat types from the expert-based map (Table 1) served as input for generating the annotation dataset.
We performed a sensitivity analysis to test how the number of training points influences the overall accuracy of the Random Forest algorithm. For this, we randomly placed either 100, 500 or 1000 annotation points in each vegetation class (as defined in Table 1). We tested between 100 and 1000 annotation points because 100 points could be feasibly collected in the field using a GPS device (e.g. for applying the workflow in another area) whereas 1000 points provided a very dense sampling of each vegetation class.

Workflow
To process, analyse and classify the ALS data, we developed an open-source workflow ( Fig. 2) with four main routines: (1) pre-processing of the ALS data, (2) calculation of lidar metrics, (3) selection of lidar metrics and (4) classification using the Random Forest algorithm. We used the R package lidR (https://github.com/Jean-Roma  in/lidR) for handling and processing the ALS data and for calculating the lidar metrics, and additional R libraries (raster, randomForest, caret) for analysing and classifying the raster layers. The scripts for running the developed workflow are available via GitHub (https://github.com/ eEcoLiDAR/PhDPaper1_Classifying_wetland_habitats).

Pre-processing of ALS data
The first routine included three pre-processing steps: (1) tiling, (2) water masking and (3) normalization of the height above the ground (Fig. 2). First, the ALS data of the whole study area was tiled into 2.5 km 9 2.5 km tiles to allow parallel processing within the lidR environment. Second, a water mask was derived based on the official Digital Terrain Model product (based on AHN2 at 0.5 m resolution obtained from https://ahn.arcgisonline.nl/ahn viewer/). In this product, not-available (NA) values indicate the presence of water and we used this information to exclude water bodies from our dataset. Third, the height (z value) of the point cloud was normalized using the officially published ground classification product. We selected the nearest 20 ground points in the neighbourhood of each non-ground point and interpolated a local terrain model using the inverse distance weighting interpolation method with default settings (i.e. power equal to 2). The height (normalized z value) of each non-ground point was then subtracted from the height of the underlying terrain surface to estimate the absolute (normalized) height of each non-ground point in the point cloud.

Calculation of lidar metrics
In the second routine of the workflow (Fig. 2), we calculated thirty-one lidar metrics (Table 2). These lidar metrics were calculated at four different spatial resolutions (1 m, 2.5 m, 5 m and 10 m grid cells). The 31 lidar metrics (Table 2) belong to six feature classes (vegetation cover, 3D shape, vertical variability, horizontal variability, height and topography) which reflect different aspects of the 3D vegetation structure and microtopography. All metrics were calculated using the normalized point cloud. Furthermore, buildings were removed from the datasets before the analysis by using a cadastre dataset from the Netherlands as a mask (https://zakelijk.kadaster.nl/-/top10nl). Finally, we included ground points for lidar metrics at level 1 (planar surfaces vs. vegetation) and for the topography-related metrics at all levels. The code for the implementation of each lidar metric is available from GitHub (https://github.com/ eEcoLiDAR/PhDPaper1_Classifying_wetland_habitats).

Lidar metrics selection
In the third routine of the workflow (Fig. 2), the calculated lidar metrics were analysed with the Random Forest algorithm (Breiman, 2001). This was done in three steps.
First, a sensitivity analysis was performed to identify the influence of (1) the spatial resolution of the lidar metrics (1, 2.5, 5 and 10 m), and (2) the number of annotation points on the overall accuracy of the classification. We used the randomly placed annotation points (100, 500 and 1000 points per vegetation class), split them into training (75%) and testing (25%) data, and then used the Random Forest algorithm (run 100 times and splitting the annotation data each time randomly) to classify land cover types and habitats at each of the three levels (level 1-3). The number of trees was set to 100 and the number of variables randomly sampled at each split was set to the square root of the input variables across all RF runs. We applied the RF using lidar metrics calculated either at 1 m, 2.5 m, 5 m or 10 m resolution. We then calculated the mean AE SD of Overall Accuracies  (Congalton, 1991) of the Random Forest algorithm for each spatial resolution, each number of testing points and each classification level. Subsequent analyses were only done with the most appropriate spatial resolution and number of annotation points (see Results).
In the second step, we performed a collinearity analysis of the 31 lidar metrics using the Variance Inflation Factor (VIF) (Naimi, 2015). While the Random Forest algorithm itself is not sensitive to multicollinearity among predictor variables, the subsequent ranking of the lidar metrics is, given that highly correlated metrics have a similar feature importance value which can affect a backward feature selection procedure. We set the VIF threshold to 5 (James et al., 2013) to avoid multi-collinearity, and all lidar metrics with VIF >5 were excluded from further analyses (see Appendix S1 Table S1).
In the third step, we applied the Recursive Feature Elimination (RFE) procedure (Guyon et al., 2002) to identify the most important lidar metrics. The RFE procedure uses the Random Forest algorithm with a backward selection of lidar metrics based on ranking the features (i.e. lidar metrics) by their importance. The feature importance was assessed using the mean decrease in impurity (MDI) (Breiman, 2001). The MDI is defined as the total decrease in node impurity, averaged over all decision trees within the random forest algorithm (https://cran.r-project.org/web/packages/randomForest/random Forest.pdf). We parametrized the RFE using cross-validation by dividing the annotation points into training (75%) and testing (25%) data (100 times). We then extracted the mean and confidence interval (CI) of the accuracy values across the 100 runs. The RFE then identifies the minimum set of lidar metrics, using a 5% loss in accuracy around the highest accuracy value. This analysis was done separately for each of the three hierarchical levels.

Classification
In the fourth routine of the workflow (Fig. 2), we used the selected set of lidar metrics at each hierarchical level to classify different land cover types and habitats (Table 1). We first build the Random Forest algorithm with 100 decision trees using randomly split points for training (75%) and validation (25%) over 100 runs. We then used the confusion matrix to assess the accuracy of the classification. We note that the accuracy of the classification is quantified relative to an expert-based vegetation map which also contains unreported classification errors. For the error assessment of the classified vegetation map, we calculated the Producer's Accuracy (precision) and User's Accuracy (recall) for each vegetation class, as well as the Overall Accuracy and Cohen's Kappa values (Congalton, 1991). We finally used the Random Forest algorithm to make a prediction and map of the vegetation classes across the wetland. For the most important lidar metrics (as identified by the RFE) we used partial dependence plots (Friedman, 2001) to visualize the relationship (i.e. response curves) between the lidar metrics and the predicted probability from the Random Forest algorithm. Partial dependence plots illustrate the probabilities (ranging from 0 to 1) that a particular vegetation class is present along the gradient of the selected lidar metric. The partial dependence plots can therefore support the interpretation of the ecological relevance of each lidar metric.

Lidar metrics selection
The sensitivity analysis (Table 3, Appendix S1 Fig. S1) confirmed that the highest accuracies for the classification were achieved when the spatial resolution of the lidar metrics (5 m) matched the spatial scale of the annotation data (1:5000 map). An increase in the number of annotation points for the classification at 5 m resolution resulted in a 2% increase in Overall Accuracy at level 3 (79% for 100, 80% for 500 and 81% for 1000 annotation points), and a 1% increase at level 1 (86% for 100 and 500 and 87% for 1000 annotation points). Hence, a smaller number of annotation points still provided good classification results. Moreover the use of 1000 points resulted in a dense clustering of points, especially for vegetation classes with a small areal coverage, thus introducing a high level of pseudo-replication. We thus opted for 500 annotation points and a 5 m spatial resolution for the remaining analyses.
The results of the collinearity analysis suggested that 17 out of 31 metrics were highly correlated with other metrics (VIF > 5). The remaining fourteen metrics (Table S1) had VIF < 2.3. The maximum Pearson correlation in the remaining lidar metrics was r = 0.6 between linearity (S_lin) and density of vegetation points below 2 m (C_b2). All other lidar metrics had lower correlations (r < 0.6).
The RFE procedure identified a minimum set of lidar metrics (out of the 14) for the classification at each of the three hierarchical levels (Fig. 3a-c). At level 1 (Fig. 3a), the first 2 lidar metrics already achieved 82 AE 2% Overall Accuracy (mean AE CI, n = 100 runs). These two most important lidar metrics represented vegetation height measured as 25th percentile of normalized Z (H_25p) as well as horizontal variability measured as variation of the DSM (HV_var) (both with MDI values >5; Fig. 3d). At level 2, three lidar metrics were sufficient to separate different wetland-related land cover types, reaching 71 AE 2% Overall Accuracy (Fig. 3b). The three metrics represented measures of horizontal variability (HV_var), vertical variability (VV_var) and vegetation height (H_25p) (MDI values >9; Fig. 3e). At level 3, the  Wilson et al. (2007) Terrain roughness* T_rough The maximum difference of the digital terrain model based on 4 neighbouring pixels The feature class is grouping lidar metrics according to the type of 3D vegetation structure (or topography) that they represent. Metric abbreviations include the feature class abbreviation (C = cover, S = 3D shape, VV = vertical variability, HV = horizontal variability, H = height, T = topography) and an acronym for each lidar metric name. The asterix (*) behind the name of the lidar metric indicates which metrics were included after the collinearity analysis (see Methods and Appendix S1 Table S1). The eigenvalues were calculated based on the structure tensor which describes the principal directions of the points within the raster cell by determining the covariance matrix of the x, y and z coordinates. The specific implementation and calculation method of each lidar metric is available from GitHub (https://github.com/eEcoLidar/PhDPaper1_Classifying_we tland_habitats). DSM, digital surface model; z, normalized height value.
separation of different types of reedbeds required four lidar metrics (Fig. 3c) and achieved 78 AE 2% Overall Accuracy. The four selected lidar metrics represented vertical variability (VV_var), two metrics related to vegetation cover (C_can and C_b2) and one related to vegetation height (H_25p) (MDI values >5; Fig. 3f).

Classification
The fine-scale land cover and habitats predicted from the Random Forest algorithm with the minimum set of lidar metrics was mapped across the wetland (Fig. 4). A visual comparison of the lidar-based prediction (Fig. 3) and the expert-based vegetation map (Fig. 1) suggested that the general distribution of forest, shrub, and reedbed habitats was similar. One apparent misclassification was that structurally poor land reed was wrongly classified as structurally rich land reed on the lidar-based map in the southern part when compared to the expert-based map (E6°11'', N53°21'') (compare Fig. 1 and Fig. 4). The Overall Accuracy of the classifications (Table 4) ranged from 71-82%, and Cohen's Kappa ranged from 0.61 to 0.67. At level 1, the User's and Producer's Accuracy were both high (>81%; Table 4). At level 2, forest was classified with a Producer's and User's Accuracy of >83%, but grassland, reedbed and shrub achieved lower Producer's and User's Accuracy (41-75%; Table 4). At level 3, the Random Forest algorithm achieved the separation of structurally rich land reed from the other two reedbed classes with ≥92% Producer's and User's Accuracies (Table 4). The structurally poor land reed and the water reed were classified with lower success (53-78% Producer's and User's Accuracy; Table 4).
The partial dependence plots of the most important lidar metrics revealed how land cover types and habitats within the three hierarchical classification levels were separated along gradients of vegetation structure (Fig. 5). At level 1, vegetation and planar surfaces clearly differed in their 25th percentile of normalized height (H_25p), with planar surfaces having a high probability below a threshold of 0.5 m (Fig. 5a). At level 2, the variance of the DSM (HV_var) separated different wetland-related land cover types, with reedbed and grassland showing a high probability when variability in vegetation height across neighbouring cells was low (<2.5 m; Fig. 5b). Instead, shrub and forest showed a high probability when HV_var was >2.5 m (Fig. 5b). Grassland and reedbeds also showed a high probability when the variance of height (VV_var) was <2.5 m whereas shrubs and forest had a high probability >2.5 m (Fig. 5c). At level 3, water reed could be distinguished from both land reed classes when the variance of normalized height (VV_var) was >7.5 m (Fig. 5d). Water reed was also well separated from structurally rich land reed when the 25 th percentile of height (H_25p) was <2 m (Fig. 5e) or when the density of vegetation below 2 m was >10% (Fig. 5f).

Discussion
We tested 31 lidar metrics for classifying and mapping land cover types and habitats in a Dutch wetland using country-wide ALS data. The RFE procedure together with the hierarchical classification approach revealed that five lidar metrics were of key importance to separate fine-scale habitats within this wetland. These lidar metrics reflected horizontal variability (HV_var), vertical variability (VV_var), height (H_25p) and cover (C_can, C_b2) of the vegetation. The results confirmed our expectation that the total set of lidar metrics can be substantially reduced for classifying wetland-related land cover types and habitats, and that lidar metrics related to vegetation height and horizontal as well as vertical Table 3. Results of the sensitivity analysis based on a Random Forest algorithm using varying amounts of annotation points (100, 500 and 1000) at different spatial resolutions (10 m, 5 m, 2.5 m, and 1 m). The Overall Accuracies (mean AE SD across 100 runs) are given for three hierarchical levels (see Table 1, also for the vegetation classes). At each resolution, the 100, 500 and 1000 annotation points were placed randomly in each vegetation class. At 1 m resolution, only 100 and 500 annotation points were tested. A visualization of the data from this table is presented in Appendix S1 Figure S1.   . Fine-scale land cover and habitat map of the Lauwersmeer area predicted from one example tree of a Random Forest algorithm using lidar metrics derived from country-wide Airborne Laser Scanning data. Classification was performed separately at three hierarchical levels (levels 1-3) and predictions were then combined into one map, with colour coding showing the predicted land cover types and habitats at the finest hierarchical level (white boxes in the legend indicate that the respective habitat type is classified at a finer level). See Table 1 for the vegetation classes, and Figure 1 for the expert-based vegetation map that was the basis for creating training and testing data.

Lidar metrics selection
The spatial resolution of lidar metrics may affect the classification accuracy. Many lidar-based wetland studies have used spatial resolutions below 5 m (Onojeghuo et al., 2010;Zlinszky et al., 2012;Luo et al., 2015;Kope c et al., 2016;Chasmer et al., 2016;Corti Meneses et al., 2017). However, the expert-based vegetation map which we used for annotation in our study does not capture the fine-scale heterogeneity of vegetation structure within polygons. This was confirmed by our sensitivity analysis which showed that classifications at 1 m or 2.5 m resolution showed less accuracy compared to classifications at 5 m or 10 m resolution. Other annotation data (e.g. from vegetation plots) may allow better classification accuracies at resolutions below 5 m (e.g. Alexander et al., 2015;Corti Meneses et al., 2017). Our results therefore suggest that the best use of the available expert-based vegetation map for classification in our study is at a moderate (5-10 m) spatial resolution. Studies using expert-based maps for annotation could additionally explore fuzzy classification methods which may allow to handle classification uncertainties along vegetation class boundaries (Zlinszky and Kania, 2016). In our study, the availability of a wall-to-wall expert-based vegetation map allowed us to test the influence of the number of annotation points on the overall accuracy of the classification. Such expert-based maps allow to generate hundreds or thousands of annotation points (Millard and Richardson, 2013), while studies using field-based annotation data usually rely on less because they have to be collected in situ (Alexander et al., 2015;Kope c et al., 2016). Our sensitivity analysis implies that using 1000 instead of 100 annotation points only slightly increases the overall accuracy of the classification (from 86% to 87% at level 1, from 70% to 73% at level 2 and from 79% to 81% at level 3). This suggests that a smaller amount of annotation points (e.g. 100) could be sufficient to classify vegetation classes within wetlands in the Netherlands. This would make the collection of in situ measurements for annotation in other wetlands feasible. It would only require to record the land cover type or habitat with a GPS at each annotation point which could be rapidly done within a few days at a given study site. The collinearity analysis in combination with the RFE procedure allowed us to reduce the number of lidar metrics in the classification at each hierarchical level. A similar approach to metric reduction has been used in other lidar-based vegetation classifications (Millard and Richardson, 2013;Alexander et al., 2016;Koma et al., 2016). In our study, only five lidar metrics related to vegetation structure where needed to separate different land cover types and habitats within this wetland. Similar to the wetland study of Zlinszky et al. (2012), we found that vegetation height (25th percentile of height) and horizontal variability (variance of the DSM) were the most important metrics to separate planar surfaces from vegetation. However, 3D shape metrics did not play a particularly important role in separating planar surfaces from vegetation. This is in contrast to studies in urban settings (Weinmann et al., 2017) and agricultural landscapes (Lucas et al., 2019) where neighbourhood-based features based on eigenvalue characteristics of the point cloud seem to be crucial for identifying planar surfaces. For separating wetland-related land cover types such as forest, grassland, reedbed and shrubs, the variance of height Table 4. Classification accuracies based on the confusion matrices for each of the three hierarchical levels (mean AE SD across 100 runs) using 500 annotation points and lidar metrics calculated at 5 m resolution. became additionally important. This mostly reflects the key structural differences between tall and low-stature vegetation (Alexander et al., 2015;Ene et al., 2018). Finally, metrics related to vegetation cover (e.g. canopy cover and density of vegetation below 2 m) together with metrics capturing vertical variability (variance of height) and vegetation height (25th percentile of height) were crucial to separate water reed from land reed. This supports evidence from other studies that have used point density, height and variability of the DSM to quantify the structure and extent of reedbeds along lakeshores (Millard and Richardson, 2013;Corti Meneses et al., 2017). The reduced number of lidar metrics facilitated the ecological interpretation of the Random Forest algorithm by focusing on the most important lidar metrics and their specific response curves. The partial dependence plots revealed particular thresholds along gradients of vegetation structure that separate one type of vegetation class from another. For instance, the separation of water reed from land reed showed thresholds for land reed along gradients of variance of height (VV_var > 7 m), 25th percentile of height (H_25p between 0 and 2.8 m) and density of vegetation below 2 m (C_2b > 18%). Clear thresholds also existed for separating planar surfaces from vegetation (H_25p < 0.4 m) and for using horizontal variability for separating grasslands and reedbeds (HV_ var < 2.5 m) from shrubs, forests patches and individual trees (HV_var > 2.5 m). These identified thresholds together with repeated ALS surveys could guide wetland management and monitoring, for example, in terms of vegetation structural changes that are important for breeding birds. Furthermore, the most relevant information on vegetation structure derived from ALS (i.e. lidar metrics related to vegetation height and vegetation density) corresponds to the criteria used for deriving the expert-based vegetation map.  Table 2 for details of the lidar metrics and Figure 3D-F for the feature importance of each metric in the classification.

Classification
For wetland vegetation classifications, some studies have reported >85% Overall Accuracies (Goodale et al., 2007;Gilmore et al., 2008;Zlinszky et al., 2012;Kope c et al., 2016;Shadaydeh et al., 2017;Rapinel et al., 2018). This is slightly higher than in our study and could be explained by the inclusion of spectral information in the classification (e.g. Onojeghuo and Blackburn, 2011;Zlinszky et al., 2012;Alexander et al., 2015). In comparison to other sensors (e.g. optical), lidar captures the structural aspects of vegetation which is its main advantage. The ALS data were acquired in the leaf-off season which could result in lower accuracies because the variability of vegetation structure might not be fully captured (Onojeghuo et al., 2010;Wasser et al., 2013). To the best of our knowledge, the effect of seasonality on capturing vegetation structure from ALS data (i.e. leaf-off vs. leaf-on season) has not yet been investigated within wetlands, and needs to be addressed in the future. A further classification error could be introduced by the use of the expert-based vegetation map for annotation because it is derived from an interpretation of 3D stereoscopic aerial photographs together with field visits. This means that we validate our automatic (ALS-based) classification against a vegetation map derived from human cognitive processes. This may contain unreported classification errors. Furthermore, there was also a time difference between the mapping of the annotation data (2015 and before) and the collection of the lidar data (2008-2009) which could further lower the classification accuracies. Nevertheless, the achieved accuracies (71-78%) for separating level 2 and level 3 vegetation classes (i.e. different wetland-related land cover types and reedbed habitats) correspond to other wetland studies in which accuracies range from 57% to 88% (Goodale et al., 2007;Zlinszky et al., 2012;Millard and Richardson, 2013;Alexander et al., 2015;Corti Meneses et al., 2017). The combination of lidar data with hyperspectral information (Onojeghuo and Blackburn, 2011) or with Synthetic-Aperture Radar (SAR) (Millard and Richardson, 2013;Montgomery et al., 2019) could further improve the classification of land cover and habitat types.

Conclusion
With the increasing availability of open-access, countrywide ALS datasets, our developed workflow offers a promising step towards monitoring vegetation structure in wetlands at fine spatial resolution and over broad spatial extents. Such monitoring of habitat structures is essential for developing an observation system that informs decision makers on global biodiversity change (Navarro et al., 2017). Most of the currently existing studies using ALS data to classify wetland-related land cover types and habitats have so far focused on relatively local sites (<100 km 2 ), with the exception of two studies covering wetland areas of >1000 km 2 (Zlinszky et al. 2012;Chasmer et al. 2016). A main bottleneck for upscaling local ALS-based habitat studies is the lack of consistent in-situ annotation data within and across countries and regions. The use of expert-based vegetation maps (as in our study) is limited for upscaling because mapping schemes and their vegetation units often differ among sites. Nevertheless, additional annotation data could be derived from other sources such as vegetation inventories (e.g. georeferenced vegetation plot data obtained from standardised field sampling protocols), or from sampling information on land cover and habitat types with a GPS across multiple wetland sites. Furthermore, in many cases the collection of field data is not necessarily done simultaneously with the lidar flight campaigns, which introduces some additional uncertainty in the cross-validation procedure. A stronger coordination between communities of ecologists and remote sensing scientists is therefore desirable for monitoring biodiversity change and for the conservation and management of wetlands and other habitats.