remote sensing Categorizing Grassland Vegetation with Full-Waveform Airborne Laser Scanning: A Feasibility Study for Detecting Natura 2000 Habitat Types

: There is increasing demand for reliable, high-resolution vegetation maps covering large areas. Airborne laser scanning data is available for large areas with high resolution and supports automatic processing, therefore, it is well suited for habitat mapping. Lowland hay meadows are widespread habitat types in European grasslands, and also have one of the highest species richness. The objective of this study was to test the applicability of airborne laser scanning for vegetation mapping of different grasslands, including the Natura 2000 habitat type lowland hay meadows. Full waveform leaf-on and leaf-off point clouds were collected from a Natura 2000 site in Sopron, Hungary, covering several grasslands. The LIDAR data were processed to a set of rasters representing point attributes including reflectance, echo width, vegetation height, canopy openness, and surface roughness measures, and these were fused to a multi-band pseudo-image. Random forest machine learning was used for classifying this dataset. Habitat type, dominant plant species and other features of interest were noted in a set of 140 field plots. Two sets of categories were used: five classes focusing on meadow identification and the location of lowland hay meadows, and 10 classes, including eight different grassland vegetation categories. For five classes, an overall accuracy of 75% was reached, for 10 classes, this was 68%. The method delivers unprecedented fine resolution vegetation maps for management and ecological research. We conclude that high-resolution full-waveform LIDAR data can be used to detect grassland vegetation classes relevant for Natura 2000.


The Natura 2000 Network
Natura 2000 is a network of protected nature areas established by all European Union (EU) member states in response to the Fauna-Flora-Habitats directive of the European Council in 1992 [1] as an implementation instrument of the 1979 Bern Convention on the Conservation of European Wildlife and Natural Habitats.The member states have committed to monitor these areas in order to follow their conservation status, and to report this every six years to the Commission.Remote Sensing is expected to play and increasing role in standardizing and streamlining this monitoring procedure [2].

Conservation of Grasslands
Grasslands are vital elements of the European landscape [3,4], they have traditionally been managed by grazing and mowing for centuries [5][6][7].Grasslands also deliver essential ecosystem services such as erosion protection, groundwater recharge, carbon sequestration, and recreation [8].From the nature conservation point of view, grasslands are crucial in maintaining landscape-scale habitat and species diversity [3,9].Species richness is outstandingly high in many types of grasslands compared to other habitats, especially at small grain sizes [9][10][11].In addition to the high overall species diversity, grasslands harbour several rare plant and animal species [4,12,13]; 18.1% of the endemic species of Europe can be found in grasslands [14].Their current total extent is 1.80 million km 2 in Europe [15] with 608,410 km 2 in the EU [16].However, due to increasing urbanization and intensification of agriculture resulting in cease of traditional land use and also due to infrastructure development leading to fragmentation, the area of grasslands in Europe is constantly decreasing [3,[17][18][19].

Lowland Hay Meadows
Lowland hay meadows are one of the most widespread semi-natural grasslands of Europe.They are included in the Annex I of the Habitats Directive of the European Union [1] as habitat of community interest: habitat type 6510, lowland hay meadows (Alopecurus pratensis, Sanguisorba officinalis).In the European Union, these habitats cover 4402 km 2 (fifth largest extent of Natura 2000 grassland habitats) with 222 km 2 in Hungary [20].
Lowland hay meadows are species-rich mesophilous grasslands formed on lightly or moderately fertile soils [21].They are typical for the higher parts of floodplains, which are only rarely inundated, and also in lowland areas with a relatively high groundwater level.They are one of the most diverse habitats of Europe and harbor several protected and endangered species [5,13,22].The species pool of lowland hay meadows is rather variable, depending on the geographical location, neighboring habitats and management.The tall grass layer (up to 120-150 cm above ground according to our field observations) is dominated by Arrhenatherum elatius; Poa pratensis and Leucanthemum vulgare are typical for the middle layer (50-60 cm above ground) and the lower layer (0-25 cm) consists mainly of the ground leaves of rosette plants such as Plantago lanceolata and Ranunculus species.
Generally, lowland hay meadows were created by forest cutting and were developed and maintained by mowing, typically twice a year: early summer and late summer/early autumn, which contributes to a species-rich vegetation community [21,23].The major threats for these habitats are changes in land use, such as afforestation or conversion to arable lands.Changes in the intensity of management (both intensification and neglect) lead to the decline of species richness [24,25].Abandonment of traditional mowing regimes is most typical in Central and Eastern Europe [22,26] due to considerable decrease in livestock numbers and accordingly a decreased demand for hay production [27][28][29].Discontinuing management rapidly alters habitat structure and ecosystem functions, resulting in a decline of total species richness [22,28,30].This is mainly caused by organic litter accumulation [31] and encroachment of herbaceous and woody competitor species [18,32,33].

Grassland Habitat Assessment
Mapping and assessment of grasslands for conservation, management and research is an especially tedious task due to fine-scale dynamics, high sensitivity of vegetation to season, and often the lack of clearly defined categories [34,35].Typical methods for grassland monitoring and assessment are based on phytosociological relevés, repeated according to random, regular or nested sampling layouts [36].While such surveys provide a good understanding of vegetation composition and patterns within the surveyed samples, they cover small areas and can be difficult to extrapolate to wider areas [37].For conservation purposes, ideally, both vegetation type and management regime should be mapped [38][39][40].There is a clear demand for robust and repeatable vegetation mapping tools in grasslands that can aid and complement field methods [41,42].This demand is further increased by the requirement of the Habitat Directive [1] that member states shall monitor the status of their protected areas and report to the Commission every six years: repeatable, reliable, efficient, and area-covering monitoring technologies are needed for this.

Remote Sensing in Habitat Monitoring
Remote sensing has often been recommended as a way to obtain broad-scale ecological data, and typically allows high repeatability even over large areas [43].The high spatial heterogeneity and seasonality of grasslands continues to challenge remote sensing technology, with very few tools routinely applied for mapping or assessment [2,44].At the current level of technology, most Natura 2000 assessments are done by field mapping, but, due to the need for harmonization between member states [2,45], an increasing involvement of remote sensing is foreseen [36,39,46].The most frequent (but probably least published) method for remote-sensing assisted grassland mapping is using an aerial photo for aiding navigation and as a basis for drawing the final map [41].Georeferenced aerial images have been used in a GIS environment for outlining grasslands [47], and with the onset of high resolution multi-temporal satellite images, habitat mapping schemes, based on such data, have grown to include various types of grassland as well [48,49].

Remote Sensing of Grassland Vegetation
Remote sensing studies of grasslands are relatively rare compared to other habitats [36].Recent surveys of grassland or scrubland vegetation were typically carried out with multispectral or hyperspectral sensors [50][51][52].Water content of the leaves was identified from spectral signatures [53,54], and based on the resampling of field spectra, it is anticipated that satellite images would contribute to mapping floristic gradients [55].On hyperspectral images, up to three classes of grassland vegetation and other, non-grassland classes were differentiated [56].Object-oriented image analysis of multispectral aerial or satellite images [57][58][59], and classification of high temporal resolution satellite images have also been used for similar purposes [60].Heathlands were successfully analyzed for conservation-related mapping by hyperspectral and multispectral images, and these studies also included grassland vegetation classes [50,58,61,62].Using multi-temporal RapidEye and separately TerraSAR-X data, Schuster et al. [60] and Neumann [49] mapped grassland, reaching very high accuracies (Kappa > 0.8), and Franke et al. [40] successfully categorized grassland use intensity (four classes, overall accuracy > 80%).While these results are certainly ground breaking in their own field and cover a certain aspect of monitoring (detecting mowing regime for the radar-based studies, species composition principle components for the hyperspectral study), they are less compatible with the vegetation classification schemes normally used in conservation mapping, and their spatial resolution is also limited to several meters.

LIDAR for Ecological Mapping
LIDAR (Light Detection And Ranging, also referred to as Airborne Laser Scanning, ALS) has acknowledged potential for mapping canopy structure but it is rarely considered an alternative to imaging methods for vegetation classification [63].Instead of a more intuitive "picture", LIDAR samples terrain elevation and vegetation height and creates a dense set of points [64].However, with the onset of full waveform recording and later radiometric calibration, it was proven that vegetation classification relying only on LIDAR data is possible [65].LIDAR was first used for such purposes in forests, where vegetation height and canopy structure are more diverse [66], but applications in shrublands [67,68] and wetlands were also successful [69][70][71][72].LIDAR is increasingly used for mapping conservation relevant variables, including coverage of different species or associations for biodiversity assessment [73,74], but also human activities or natural environmental variables relevant for habitat quality [46].Meanwhile, one of the major advantages of LIDAR compared to airborne multi-or hyperspectral data is the versatile use of such data, which has led to more widespread surveys and easier data access than the case for passive multi-or hyperspectral airborne surveys.In many locations, LIDAR data is collected for purposes including topographic mapping, forestry, or in preparation for construction.Exploiting this data for automated vegetation classification of grasslands would allow cheap access to high-resolution habitat maps that could be within the reach of conservation management.

LIDAR in Grasslands
To our best knowledge, LIDAR has not been tested so far for direct vegetation classification in grasslands.In a prairie setting, LIDAR was used to quantify aerodynamic resistance based on vegetation height [75].LIDAR was also used for predictive modeling of potential grassland vegetation using water availability based on LIDAR-derived terrain height models [76].However, this approach did not involve mapping the vegetation itself based on its signature in the remote sensing data, but modeling its most probable pattern based on the governing environmental variable.All in all, LIDAR as a standalone tool was so far not considered suitable for classification of grasslands [36].

Objectives
The main objective of this study was to check the feasibility of using LIDAR as a tool for mapping different types of grassland vegetation.We also intended to establish a methodology for automatic grassland vegetation classification based on LIDAR, and to evaluate the accuracies of different surveying and processing approaches.

Study Site
Our study areas are located in the Soproni-hegység Natura 2000 site, in western Hungary (47°41ʹN, 16°34ʹE, Figure 1).The Natura 2000 site has an area of 52 km 2 in a hilly region about 150-500 m above sea level.The local climate is sub-alpine with a mean annual precipitation of 750-900 mm, and a mean temperature of 8.5 °C [77].The site was designated to preserve 14 different Annex I habitat types: mainly beech and oak forest types, dry, mesophilous and wet grasslands, fens, fringes, and small heaths [78].Meadows are present throughout the area with patch sizes ranging from 0.1 ha to 27 ha.These grasslands are mainly situated in the valleys, which are the floodplains of small rivers and in some cases on other flat, mesophilous areas.Most of the studied grasslands have an anthropogenic origin, except for the extremely wet patches in the largest meadows [77].These clear-cuts were formerly mainly orchards or vegetable gardens, and grassland development supported by traditional land use by regular mowing resulted in species-rich habitats.The most widespread grasslands in our study site are lowland hay meadows, which build a diverse mosaic structure with other grassland types such as semi natural dry grasslands with Bromus and Festuca species, wetter Molinia meadows, sedge stands, disturbed and abandoned grasslands, weed patches, and shrubs.Most of these meadows are mown twice a year (late spring and late summer or early autumn) regardless of vegetation type.

Field Data Collection and Survey
Field campaigns were carried out during May and September 2012 and April and May 2013, timed to match the local vegetation climax as determined by weather conditions of each year.Most of the grassland patches were visited several times in different years to confirm that both vegetation and management corresponds to the assigned class.Field reference data were collected in 140 homogenous vegetation patches with a minimum size of 50 m 2 , giving a total ground truth area of 16,162 m 2 for 7 km 2 of meadows within a total flight area of 90 km 2 .Sampling polygons were typically circles of 4-meter radius, or more or less rectangular patches of similar minimum area.If homogeneous areas of this size could not be found, polygons of different (more elongated) shape were also used, especially for fringes and field roads.The study area was divided into a set of 9 different studied meadow sites based on the LIDAR-derived tree mask (Section 3.6), and we aimed to capture at least one polygon of each class for each of these areas as far as possible.We specifically searched for field examples of the less frequent classes (e.g., lawn) to reach an optimum distribution of the calibration and validation data in space across the study area.We selected the methodology of the fieldwork based on the objective as a feasibility study, balancing the need for similarity to Natura 2000 assessments with efficient and focused data collection.Coordinates of the polygons were marked by a differential GPS.While full phytosociological relevés were not carried out, we registered a list of characteristic vascular plant species, and noted water availability (mesophilous, wet, dry, seasonally wet or dry patches) and the current management (mown, grazed, or abandoned) for each polygon.Polygons representing field roads and non-grassland classes such as shrubs and non-vegetated surfaces were manually digitized from the raw leaf-on NDSM (Normalized Digital Surface Model) and Reflectance LIDAR product rasters (Section 3.5) validated by differential GNSS (Global Navigation Satellite System) in the field afterwards.The classes were also defined to test whether LIDAR is suitable for grassland vegetation mapping.We did not adhere to a pre-existing classification system, such as EUNIS [79] or the Hungarian Á-NÉR system [80] since these were not developed for remote sensing, but created categories with the intention of keeping them relevant for habitat management and grassland ecology, and covering the full range of land cover encountered in the studied sites ("mutually exclusive and totally exhaustive" [81]).However, for grassland categories, we included three classes where we followed Natura 2000 habitat characterization manuals [82][83][84] as far as possible.For outlining these polygons, we strictly relied on the diagnostic species and specific structures outlined in the manual (and in the respective category definition, Section 3.3.1).Since not all grasslands fall into Natura 2000 categories, other classes had to be specified to cover areas of the grassland mosaic that do not fulfill the criteria of any Annex I habitat type.In addition to lowland hay meadows we defined seven grassland classes altogether, some of these including (but not exclusive to) other Natura 2000 habitat types.We also introduced categories for non-vegetated surfaces, and shrubs and single trees.The final set consisted of 10 classes that were hierarchically merged to five larger categories as an alternative scenario (Table 1).Half of the reference polygons were used for calibration and half for validation, based on random selection and checked in a GIS to ensure maximum independence of the calibration and validation sets."Shrub": Shrubs were defined as woody plants below 1.5 meter height (threshold of the meadow mask, Section 3.6) They are important for separation from meadow vegetation, but also as an element of habitat diversity.Shrub encroachment is often part of the degradation process of abandoned meadows.In our case, these were mostly Crataegus, Rosa or Salix species, or fruit trees.Young trees planted for forest renewal purposes were also included in this category.

Description of Classes
"Fringe": Our class fringe pools three different types of vegetation: Patches dominated by non-native neophytes, such as Fallopia japonica (Figure 2a), Impatiens glandulifera or Solidago canadensis; native fringe communities represented by nitrophilous, fresh patches of Urtica dioica (Figure 2b) and hydrophilous forb fringe communities on forest edges including Natura 2000 habitat type 6430 or Phragmites australis.Therefore, this category contains the habitat 6430 but is not exclusive to it.
Fringe vegetation is rich in or dominated by tall growing forb species, and typically occupies an intermediate setting between (mown) grassland habitats and forests or single shrubs."Abandoned": We mapped grassland patches as abandoned if there was no sign of any management during the field campaigns and species composition confirmed this.Calamagrostis epigejos dominated patches were also assigned to this category even if some of them were regularly mown, since its presence reflects deteriorating quality of the meadow.Usually, species of grasses and herbs typical for mown meadows still occur in these areas, but ruderal species, such as Tanacetum vulgare, Artemisia vulgaris, Bromus inermis are also present, and a high proportion of Dactylis glomerata is typical.Tree and shrub encroachment is also characteristic (Figure 3a,b)."Meadow like": Grasslands composed of grass and/or herb mixtures but not belonging to one of the other classes were categorized "meadow like".This is a relatively heterogeneous class: meadows where the species composition reflected degradation and did not allow identification of the original grassland association, without a meadow typical mowing regime (such as road embankments), or covered by weeds but regularly mown, were included in this category, some of them very similar to a lowland hay meadow but lacking some characteristic species and not fulfilling the Natura 2000 criteria.Examples for meadow like areas in our case are intensely used or undergrazed pastures, vegetated field tracks (Figure 4a) and an open, unmanaged areas (Figure 4b).Meadow like patches are characterized by a high cover of tall grass and fringe species, while typical meadow herb species are missing or just occur with very low coverage, or ruderal indicator species can be present.Meadow like implies occasional mowing, nevertheless, the transition to the species composition of abandoned meadows is fluent.[84,85].Since this definition refers to a specific combination of species composition, management and abiotic parameters, the term lowland hay meadows was applied strictly and only to areas that truly corresponded to all criteria (within the limits of identification) in the field (Figure 5a).A detailed species list was collected for all patches categorized as lowland hay meadows, and most such areas were visited several times in order to confirm that they are correctly managed.Nevertheless, field identification was often problematic, since at the studied spatial scale of a few meters, lowland hay meadows form smooth transitions with dry meadows and Molinia meadows.Differences in definition of the habitat category according to national Natura 2000 evaluation schemes were also a problem, but wherever contradictions were found, the Hungarian national system was applied [83,85]."Molinia": All areas where Molinia species were among the dominant grasses belong to this class, which therefore roughly corresponds to the Natura 2000 category 6410 (Planar to montane Molinia meadows [82][83][84]).In our study site, Molinia patches occurred in the transition zones between wet grassland types and lowland hay or dry meadows which are wet in spring but dry in summer, sometimes without developing a well-defined association.
Molinia can form small patches or mingle with other grasses but shows a characteristic tussock shape (Figure 7).Unlike most of the other grass species, Molinia develops quite late in the year and reaches its maximum height with the flower stands only by the end of summer.[82][83][84]) and Viscario-Festucetum rubrae [86] patches belong to this category.These xerotherm communities grow under warmer and drier conditions than lowland hay meadows, often on southern exposed slopes.Dry meadows usually grow lower in height than lowland hay meadows and the layer of tall grasses is not that dense, therefore they are usually rich in short growing, flowering herbs (Figure 5b).Both types form fluent transitions to the lowland hay meadows in our study site.
"Lawn": Lawns are frequently mown grasslands (more than three times a year).They are generally species poor, sown grassland areas.Lawns are much shorter in height than all other meadow types.Typical lawn areas were playgrounds, sports grounds and gardens.

Set of Five Classes
The set of only five classes was composed by hierarchically fusing the previously defined categories (Table 1).It is more focused on meadow identification for management and more specific to lowland hay meadows.
"Not vegetation": Not vegetation corresponds to the not vegetation class of the 10 category scheme."Shrub": Shrub is exactly the same as the shrub class of the 10 category scheme."Not mown": Not mown vegetation is composed from the categories "Fringe", "Abandoned", and "Meadow like" and thus includes all patches dominated by mainly unmanaged non-woody vegetation.
"Lowland hay meadow": Lowland hay meadow is the same as for the 10 category survey, thereby focusing the five categorization scheme on the detection of this key habitat."Mown meadow": Mown meadow vegetation was fused from the following classes: "Wet high", "Molinia", "Dry meadow", and "Lawn".This category covers all regularly managed grasslands except lowland hay meadows.

Airborne Survey
LIDAR data were acquired in two separate flights, both covering the same study area, under early spring (March 2012) and summer (July 2011) conditions, which we respectively call "leaf-off" and "leaf-on" datasets in the following.The data were collected in a broader project with the main purpose of forest mapping, which explains why the dates are less than optimal for grassland mapping.A Riegl LMS-Q680 sensor was used [87], flown at an altitude of ca.500 m above ground, with a pulse repetition rate of 400 kHz and 50% strip overlap.The emitted laser pulse had a wavelength of 1550 nm, which means the reflectivity is strongly influenced by water content of the plant cell structure [54].The sensor had a footprint diameter of 0.25 m and full waveform recording was applied [65].The resulting point density was 12.8 echoes/m 2 .3.5 LIDAR Data Product Rasters Data products (point cloud "derivatives") were all generated in a 0.5 m × 0.5 m raster size, using the Opals software package [88,89].0.5 m × 0.5 m is a widely accepted plot size for grassland studies [83] and has been used successfully before for grassland remote sensing [57].The final set of LIDAR-derived raster data was the following (Table 2): Mean surface reflectance: The LIDAR sensor emits a short pulse of light, and the amplitude of the reflected signal is influenced by the reflectivity of the target surface.In case of vegetation, calibrated LIDAR reflectance has already been proven to contribute significantly to classification [70], but it is also known that reflectance is only partly correlated with species composition and is influenced by plant health and chemical structure together with soil surface parameters [52].Based on homogeneous reference surfaces with known reflectance (asphalt roads), a mission-specific calibration constant was calculated and LIDAR amplitude was converted to the real surface reflectance coefficient using the OpalsRadioCal software module [90,91].The mean reflectance value of all points within the raster cells was output as a result.
Echo width: The outgoing laser pulses are characterized by a uniform width in time, while the width of the backscattered echo is modulated by properties of the target surface within the footprint [92].If the area illuminated by the pulse footprint is flat, smooth, and homogeneous in height, the recorded signal closely resembles the emitted pulse.If there are differences in elevation within the illuminated footprint, this results in the energy of the pulse being spread over a longer time in the echo.Therefore, areas with heterogeneous vegetation height or a rough canopy or terrain surface will have higher echo widths than smooth and flat terrain [93].Echo width was stored as an attribute for each recorded LIDAR point during the airborne survey.We calculated the mean echo width for each 0.5 m cell.
NDSM height: In typical vegetation mapping studies, the most important LIDAR-derived parameter is the vegetation canopy height [57,94,95].This is represented by the Normalized Differential Surface Model (NDSM), which is generated by creating both a Digital Terrain Model (DTM) from the LIDAR points reflected from the actual terrain surface, and a Digital Surface Model (DSM) from the points corresponding to the top of the canopy, and then subtracting the DTM height from the DSM height.In our case, we used hierarchical robust filtering implemented in the SCOP++ software package [96] to generate a DTM (cell size 0.5 m) from the leaf-off data.In order to avoid eventual noise points, the second highest point within the 0.5 m cells was rasterized to create a DSM separately for leaf-off and leaf-on point clouds and the DTM was subtracted from each.
LIDAR-derived canopy height models generally slightly underestimate vegetation height [97], since the echoes do not always sample exactly the highest point of the canopy.Therefore, NDSM height was not expected to deliver the absolute vegetation height, but to represent it consistently enough for accurate classification.
Sigma Z: Surface roughness is relevant for vegetation classification for two different reasons.On one hand, it characterizes the ability of the sensor to penetrate the canopy, with strong penetration resulting in point returns both from the canopy surface and the terrain below.On the other hand, it quantifies the local-scale variation in canopy height, referring to a neighborhood of points (d = 3 m) instead of the points in a single raster cell (as is the case for echo width).Both of these variables (penetration and height variation) characterize different plant growth strategies [70].During moving planes interpolation we used for calculating this variable, an inclined plane was fitted to the eight nearest LIDAR points within a 1.5 m distance from the central point.The standard deviation of the vertical distances of each point to the fitted plane is called Sigma Z and describes the roughness of the surface sampled by the ALS points.
Openness: Surface texture is one of the characteristic parameters of different vegetation classes, and especially clumped or tussock grasses may be well identified, based on their shape in 3D.Openness describes the roughness of a surface (represented by a raster) through a kernel filter where a cone shape is fitted to each pixel, recording the maximum opening angle where the cone touches one of the neighboring pixels [98,99].Openness was calculated both by fitting the cone from above and from below (positive and negative openness).The minimum openness is described as the angle where the cone touches the closest neighboring pixel while the maximum openness refers to the angle corresponding to the furthest.
Datasets were generated separately for leaf-on and leaf-off flight data.The cell-by-cell difference between leaf-off and leaf-on measurements was produced in all cases.
Part of the information in the LIDAR data is the pattern of point data within a neighborhood.Raster products that sample a neighborhood are a way to extract this information, and they also have a denoising property.Therefore we checked various kernel filtering algorithms and evaluated how they influence classification accuracy.The best results were produced by Bilateral Filtering [100].This algorithm preserves prominent edges while integrating neighboring pixels, using the non-linear combination of image values in a local neighborhood.It is based on measuring local image variation and creating a weighted average within a neighborhood of variable size and shape depending on the local variation of the pixel values.Therefore, instead of simply averaging pixels that are similar in terms of horizontal spatial coordinates, it works on pixels that are similar both in spatial position and in pixel value.We applied bilateral filtering with two different similarity settings to all the LIDAR-derived rasters described above before using them for classification.Therefore, for each of the leaf-on and leaf-off data products, two filtered versions and the unfiltered original were all used as separate bands for classification.The full list of LIDAR data products used for classification is available in the supplementary material.

Pre-Processing and Explorative Data Analysis
We calculated areas with canopy height above 1.5 m based on the NDSM, and after a morphological closing operation [101] (kernel size 10 pixels) they were excluded from further analysis as non-meadow areas.In order to gain a deeper understanding of the data properties and to perform a basic a priori input variable selection, the independence of the candidate LIDAR derivatives from each other was checked through covariance matrices, and the distributions of each input variable within the reference polygons corresponding to each vegetation class were calculated.Boxplots and violin plots [102] were used for comparative analysis of the distribution of LIDAR derivative values characteristic for each vegetation category.Scaled map visualizations of the input rasters were also checked together with the reference data polygons for initial qualitative assessment of feasibility.

Multi-Band Image Generation and Reference Data Processing
The different data products were processed by a script implemented in Python language.The first step was to create a multi-band pseudo-image from the data, where the bands were the different LIDAR derivatives.This was then overlain with the vectors of the field reference polygons, and the pixels within the polygons clipped and assigned to the corresponding vegetation class.Next, the distribution of each vegetation class in the feature space (defined by the pixel values in the LIDAR derivative bands) was investigated and used for training a classifier.

Training a Classifier Using Machine Learning Algorithms
The set of rules for classification (called the "model") was derived using machine learning algorithms implemented in the Scikit-learn library [103].Decision trees were used for initial analysis and refining different algorithm hyper-parameter settings.An automatic decision tree classification algorithm was employed [104] that created a model predicting the class of each input pixel based on the corresponding local values of input LIDAR data derivatives (pseudo-image "bands") and optimized the various thresholds for separation between classes ("trained the model") using the ground truth pixels assigned for calibration.The trained model then classified the pixels of the ground truth polygons assigned to validation, and accuracy was assessed by comparing for each pixel the result of LIDAR-based classification and the ground truth.The validation results were summed in confusion matrices and the band importance percentages output in separate report files for each run.These percentages describe how much each of the input rasters influenced the final outcome of the classification.The insight we gained from the consecutive evaluation was used iteratively to refine the calculation settings of input LIDAR derivatives (e.g., interpolation and point filtering parameters).
After testing various machine-learning algorithms, random forest classification was found to perform best.Therefore this method was used for training the final classification models and at a later stage to create the output maps [105].In this method, random bootstrap subsets of the training pixels are selected and used to construct a large number of decision trees (the "random forest").Nodes of each decision tree are split using values from a random sample of data bands (input rasters) and the split resulting in the largest information gain is selected.Random selection of input pixel sets increases classification accuracy, decreases sensitivity to noise and minimizes over-fitting [105].The ensemble of trained decision trees is then applied to each pixel of the validation dataset.Output of all individual decision trees takes part in the voting process, which selects the final class by the majority vote.The method also outputs the probability of each class for the classified pixels.

Image Classification and Rendering
Finally, the classification model with the best validated accuracies was applied to the multi-band pseudo-image covering the whole study area, and a classified vegetation map was produced with 0.5 m × 0.5 m resolution and 10 classes.
Two different visualization methods were used: basic "hard category" visualization assigns the most probable class to each pixel, and delivers a thematic raster map in the classical sense, where each of the pixels belongs to one (and only one) of the vegetation classes represented by the pixel value in the raster.While this is the default method for vegetation mapping by remote sensing, in some cases it fails to provide an appropriate representation of the species composition pattern [106].
A possible alternative is a fuzzy visualization.Random forest classifiers estimate for every pixel the probability that it belongs to each of the classification categories.Since the categories were designed to exclude each other, in the typical case a single class has a high probability for the pixel and all the rest have low or zero probabilities.However, given the smooth transitions between grassland classes and the similarities in structure and species composition, in our case it frequently happened that two or more classes had similar and high probabilities for a particular pixel.The fuzzy visualization we, therefore, produced was based on the probability output assigned by the random forest classifiers to the vegetation classes for each pixel.The colors corresponding to the vegetation classes were blended by creating a combination of their RGB (red, green, blue) values weighted according to the class probability (as already suggested by Foody [106]).

New Method for Grassland Vegetation Mapping and Accuracies of Different Baseline Datasets
The machine learning approach together with the generation of a large number of LIDAR product rasters allowed good exploitation of the information contained in the point cloud.The concurrent use of waveform attributes and calibrated radiometric information allowed a step beyond the usual LIDAR products of vegetation height and texture, and were the key to classification into a large number of classes.The relative contribution of each LIDAR component to the classification accuracy is reported in full detail in the supplementary material both for the 5 and 10-class scheme.
The difference between leaf-off and leaf-on variables added the dimension of seasonality to the information in the image, which is especially important for tall grassland vegetation.Combining leaf-off and leaf-on data improved overall accuracy by 10 pp (percentage points) compared to using only data from one season.Noise filtering proved to be crucial since the rasters were generated from a point cloud and not from direct imaging sensors: the bilateral filtering algorithm proved to be suitable for this rather unusual data type: overall accuracy increased by 8 pp compared to the accuracy obtained from non-filtered rasters.The necessary amount of reference data was comparable to the amount needed for classifying passive optical airborne data to a similar level of detail.

Map of 10 Classes
Overall accuracy 68.0% and Cohen's Kappa is 0.64, resembling again a "good agreement" according to Altmann [107] (Table 3).The most accurate classes are (not surprisingly) not vegetation, shrubs and trees, artificial lawn and wet high, all with producer's and user's accuracy near 80%.These categories are well defined by their typical NDSM heights.The Molinia and Dry meadow classes also have producer's and user's accuracies between 60% and 80%, the category abandoned has producer's and user's accuracies close to 65%.Lowland hay meadows and meadow like areas have accuracies between 40% and 50%.The vegetation map of these 10 categories reveals the fine-scale mosaic pattern of different grassland habitats (Figure 8).Both natural patterns created by micro-relief and soil composition and artificial features such as field roads and boundaries between management regimes stand out clearly in the hard-boundary visualization.Strip boundary effects sometimes cause artefacts in the classification, but these can be detected if information on the flight pattern is available.
The fuzzy visualization introduces a further level of detail by showing smooth transitions, and allows even more detailed identification of gradients, boundaries and objects (Figure 9).

Shrub
Fringe Abandoned

Map of 5 Classes
Overall accuracy from the 5-category confusion matrix is 75.5%, with a Cohen's Kappa value of 0.66 ("Good agreement" according to Altmann [107]) (Table 4).However, it has to be kept in mind that these values are grouped across all categories, and are therefore influenced by the differences in the number of validation pixels for the 5 categories, since they were fused from the 10-category scheme where more similar amounts of validation pixels were used.Most classes have both producer's and user's accuracies above 70%, the only exception being lowland hay meadows themselves.While these only have producer's and user's accuracies around 45%, the values are well balanced (producer's and user's accuracies are similar) and indicate relatively low overestimation.This map allows quantification of mown meadow vegetation in general, together with abandoned areas, which are both of high interest for conservation and ecosystem service assessment.

Survey Flight and LIDAR Data
Meadows have pronounced seasonality due to the mowing regime.In our case, the timing and data collection parameters of the airborne survey could not be optimized with the intention of mapping grassland vegetation: for some of our studied meadows, the grass was unfortunately mown shortly before the flight.Nevertheless, the high point density and the collection of leaf-on and leaf-off data together with waveform and radiometric attributes resulted in sufficient information for automatic classification of grasslands.While we believe that survey timing is crucial for grasslands and optimal accuracies can probably be reached when the flight overlaps with the optimum in meadow vegetation development, this study represents a real-world case: LIDAR data collected for non-scientific purposes (construction preparations, topographic mapping, urban surveying) is becoming available for more and more areas, and these can be used as input data even for such detailed vegetation studies.

Ground Mapping and Reference Data
The size and distribution of ground truths we surveyed was similar to other remote sensing studies using multispectral data [108].Reference data mapping in grasslands can require several visits over the course of a year for identification of characteristic species as a basis for clear categorization of reference polygons.The trade-off between reference data coverage and detail applies to all remote sensing studies and this was no exception.We found that (not surprisingly) classification accuracy increases steeply with the number of reference pixels for each category, but for the category lowland hay meadows we may have reached a threshold where the main limit was really the similarity to other grassland categories.Full quality assessment of the fuzzy classification by fieldwork was outside our scope and would have required much more detailed reference data collection [109], in our case this was merely a visualization tool.

LIDAR Data Products and Their Importance
During explorative data analysis, it was clear that the initial variables can be grouped into three main categories: those that depend closely on vegetation height (Echo width, NDSM height, sigma Z), those that refer to vegetation height texture but are theoretically independent from height itself (openness variables) and reflectance bands, which are mostly independent from vegetation height patterns.
Reflectance: Reflectance at 1550 nm is closely related to the water content of the leaves and the soil, and to the material of artificial surfaces [54].The reflectance raster is a single-channel infrared image, and since LIDAR is an active sensor, it delivers a better signal-to-noise ratio than most passive sensors.While sun shadow effects are usually problematic in passive imaging, the LIDAR-derived reflectance raster is free from shadow, which also underlines its use for classification.Very fine differences in soil and vegetation properties are well represented in the data: we suspect that even single circular clones of Brachypodium sylvaticum could be detected within the meadow.Unvegetated surfaces typically have a much higher reflectance than green leaves, therefore especially the leaf-on reflectance separates non-vegetation surfaces well from the meadows, allowing field roads with dry bare soil to be clearly traced.Tall and relatively dry vegetation reflects stronger than mown areas or very wet patches.The difference between leaf-on and leaf-off reflectance seems to be especially sensitive to grass and herb biomass and is a strong discriminator between mown and unmown areas.
Echo Width: Echo width resembles vegetation height structure consequently, but at an even finer scale.Differences in the order of magnitude of 0.1 ns in mean echo width are typical between mown and abandoned meadows, or Molinia-dominated patches and lowland hay meadows (both mown).Since the waveform data is delivered with a precision of 0.1 ns these patterns are at the very edge of use for classification.Non-vegetated surfaces have the lowest echo widths, followed by regularly mown lawns, and grasslands that do not correspond to any Natura 2000 habitat type.Leaf-on echo width has a broader range, but leaf-off echo width also shows some characteristic patterns, especially in the unmown places and where some growth is still possible after the second mowing.Given the very small differences that are exploited for classification based on echo width, it is no surprise that strip edge effects are also prominent, and even influence the final categorization of the raster.Echo width varies slightly with the scan angle, increasing by 0.1-0.2ns towards the strip edges.This difference could not be corrected, it was probably one of the major sources of classification artifacts.
NDSM height: As expected, NDSM height is proportional to vegetation height in general, but while it is close to the true heights for shrubs and fringes where the canopy is really closed, the height of grassland vegetation is strongly underestimated.In lowland hay meadows or dry grasslands where the vegetation height at the time of flying was around 50-120 cm, we found that the leaf-on NDSM shows heights around 20-30 cm.In a meadow that was only partly mown at the time of flight, an NDSM height difference of 20 cm represented an approximate difference of 50-100 cm in vegetation height.Abandoned grasslands and wet, high vegetation still have consequently higher NDSM values than mown wet or dry meadows despite the fact that the NDSM strongly underestimates vegetation height.In many cases, small variations in NDSM height down to a few centimeters were shown to convey valuable information.While the leaf-on NDSM was more useful for classifying natural differences in vegetation height, the leaf-off NDSM showed vehicle tracks, mowing strips and historic plowing marks more clearly.However, these small values are in the error range of the LIDAR measurement (2-5 cm), and also the typical uncertainties of co-registration between leaf-off and leaf-on flights (2-7 cm).Nevertheless, since the classifier does not rely on the absolute NDSM values but on their relative differences, even these very fine scale vertical patterns could be used for categorizing the data pixels.
Sigma Z: Sigma Z does not enhance very fine differences in point heights: very flat surfaces such as asphalt roads have quite similar values to mown meadows, both in the range of 2-5 cm.However, sigma Z contributes to discriminating between taller classes such as fringes, abandoned patches and wet high vegetation, since each of these classes have their characteristic canopy height variability and typical laser signal penetration rate.Sigma Z is also less sensitive to artifacts, but it smoothes away some fine detail since it is calculated based on the neighborhood of each cell and not the cell on its own.Openness: From the openness bands, especially minimum openness showed interesting detail as a measure closely resembling the steepest slope within the immediate neighborhood of each cell.Openness is especially strong in detecting artificial boundaries between different management regimes.However, band importance calculations have shown that this parameter delivers limited information beyond the other vegetation height-related variables.
According to the band importance reports produced by the classification software (see supplementary material for an example), the most important LIDAR derivatives were typically leaf-off echo width, the difference between leaf-on and leaf-off reflectance, and leaf-off NDSM height.

Uncertainties, Errors, Accuracies
Machine learning classification of LIDAR-derived parameters delivered high-resolution maps of grassland vegetation.The 0.5 m raster resolution provides an unprecedented fine scale for grassland vegetation maps, hardly achievable by state of the art methods such as airborne hyperspectral imaging (typically lower spatial resolution, lower number of classes), hand digitizing of high-resolution aerial or satellite photos (less fine scale of mapping), hypertemporal satellite imaging (lower spatial resolution or spectral resolution) or field maps (lower area coverage).A primary source of uncertainty is that definition of classes and quantification of habitat health is generally a problem in grasslands [34].We believe that most of the inaccuracy of lowland hay meadow detection can be allotted to the fact that the habitat is not very clearly defined, and the definition itself was developed with field identification and not remote sensing in mind.The high spatial resolution also contributes to uncertainty.While large pixels tend to smooth the inherent variability of the measured objects and, thus, result in relatively similar (average) values of a given variable, using small pixels leads to more prominent differences between pixel values During initial tests, it was shown that higher accuracies can be reached with larger pixel sizes, but we decided to keep the 0.5 m × 0.5 m resolution for the sake of higher information content.
The most widespread source of error we could identify was the effect of flight strip edges on Echo Width, which led to over-estimation of Lowland Hay Meadows, Molinia and Meadow like in these areas.However, these artifacts can be well identified if the strip boundaries are known, and the fuzzy visualization can be used to interpret the correct class in most cases.
The 5-category scheme produced accuracies (overall 75%, Kappa 0.66) that are slightly lower than the repeatability of field vegetation mapping [110], but similar to hyperspectral image classification [111] or LIDAR classification of tall wetland vegetation [70] or forests [112].
The set of 10 classes also delivered accuracies (overall 68%, Kappa 0.64) in the range of many largescale remote sensing methods and also similar to the accuracy of LIDAR DTM-based plant community prediction [76].The accuracies reached in the Molinia and Dry meadow classes are especially interesting for conservation since they are closely related to the Natura 2000 habitat types 6410 and 6210 (the latter is one of the most widespread Annex I grassland habitat types in the EU [113]).The slightly lower reliability of the category Abandoned is mostly due to the smooth transition with fringes, and to the fact that this category covers a wide range from meadow patches abandoned a year ago to agricultural fields unused for a decade and encroaching with shrubs.Fringes, lowland hay meadows and meadow like areas are often difficult to delineate even by experts in the field: fine differences of species composition control their occurrence, and especially the definition of lowland hay meadows is different in various national habitat evaluation instructions [82,85,114].

Feasibility and Potential Applications of LIDAR-Based Grassland Classification
Our results show that the feasibility of LIDAR-based grassland vegetation mapping is quite similar to that of most remote sensing based vegetation mapping studies.Given a certain set of ground truths and a suitable classification algorithm, reasonable accuracies are produced.In our case, high-resolution sensor data with several additional attributes was required, together with a machine-learning algorithm not (yet) available in off-the-shelf software.Grasslands are an especially difficult target for remote sensing, and it is demonstrated that LIDAR has very strong potential for mapping them.The habitat type where this feasibility study was tested is one of the most common grassland types of Europe; therefore we expect our method to be widely applicable in different geographic settings.While this level of accuracy does not directly correspond to the needs of Natura 2000 reporting, such maps can nevertheless be used in support of dedicated field mapping and habitat quality assessment.The non-grassland classes, which are slightly more accurate, can also be used for outlining habitats of interest for focusing field visits.For Natura 2000 reporting obligations, remote sensing is expected to contribute to mapping habitat area and distribution, detection of changes and information on habitat quality [41].It is expected that in case of LIDAR mapping of grasslands, the "hard-boundary" map would be used for this purpose: it delivers clear indications of the area occupied by each vegetation category, including fine-scale patterns and boundaries (Figure 8).The maps can also be used for identifying key habitats, for modeling the preferences of protected wildlife in order to find them in the field [115], or as a basis for precise targeting of conservation management.Since grasslands are hotspots for pollinating insects [116], it is expected that these maps would also be a valuable input for quantifying pollination ecosystem service potential.
For scientific applications in the field of vegetation ecology, or for investigating the result of possible conservation measures, the fuzzy visualization is more applicable (Figure 9).It provides a more detailed view of the probability of lowland hay meadow presence, and generally incorporates a better representation of nature since in reality classes are not usually pure and boundaries are not strict [117].The fuzzy renderings represent the fine texture of grassland habitats, including the transitions, boundaries and succession sequences, and therefore deliver much of the information in the data that would be omitted by hard classification.In context of mapping habitat quality, we expect that the probability-based maps might be suitable as an indication of species composition: the higher the probability of a given habitat class, the more the local species composition corresponds to it.This is confirmed by the observations we made at mixed vegetation reference patches that were not included in the learning and evaluation process.Given its high spatial resolution, it can be correlated in space with environmental variables and quantitatively evaluated by factor analysis.True color aerial image showing a studied meadow with ground truth polygons, and the vegetation map resulting from fuzzy visualization of the LIDAR-derived grassland classification.For areas excluded from classification (forests), the vegetation map is left transparent.Note that the overall vegetation pattern is the same as in Figure 8, but fine-scale differences and transitions are represented in more detail due to the fuzzy visualization.

Recommendations for Future Studies
Since the patterns we were searching for in the data were subtle, some general caveats of airborne survey planning proved to be especially important for our study.Dual leaf-on and leaf-off flights already introduce some variability between surveys, and this is further enhanced if the ground campaign takes place at a third time.In our case, this clearly meant that the timing of the flights was less than optimal; we had more freedom in timing the fieldwork, which successfully surveyed the local vegetation optimum at least for our targeted lowland hay meadows.
We expect that better flight timing and exactly simultaneous field surveys would have resulted in better accuracy.For dedicated grassland flights in the future, it is advised to fly with high strip overlap and a full pattern of lengthwise and cross strips (contrary to the typical flight pattern of many parallel and few cross strips), to ensure accurate geodetic control.Noise removal from the point cloud is a routine processing step, but has to be done with special care in case of a grassland target, and cannot be restricted to the removal of too high or too low points: erroneous amplitude and echo width values also have to be excluded in advance.In addition, data acquisition should be timed to match exactly the mature phase of the meadows but take place before the first mowing, to ensure that the vegetation and therefore the differences between vegetation classes are fully developed [55].This requirement could not be met in our case due to the logistics of the project, but the classification accuracies were nevertheless satisfactory for a feasibility study.
We expect that the continually increasing accuracies and point densities delivered by new generations of LIDAR sensors, together with enhanced information content such as multiple wavelengths [69,118] will result in increasing use of LIDAR as a standalone tool for mapping grassland vegetation.

Conclusions
We developed a new method for classifying grassland habitats based on raster products from dual season LIDAR point clouds and machine learning classification.The method was validated on more than 35000 pixels independent from the training data.For five categories, this initial feasibility study has delivered accuracies directly comparable with conventional remotely sensed methods.For 10 categories (eight of them grassland types), and five categories (three of them grassland) the statistical evaluation also delivered a good agreement with field references (Kappa = 0.64 and 0.68 respectively).This means that LIDAR-based mapping is in good agreement with the field-based method that is currently used for mapping and monitoring.The most important novelty of the study is that it demonstrates that LIDAR can be successfully used for grassland classification.Our results show that full-waveform radiometrically calibrated multi-season LIDAR data holds sufficient information for mapping a large number of categories, and that random forest classification is a suitable tool for extracting this information.The spatial resolution reached in this study is unprecedented for automatic remote sensing classification of grasslands, and is close to the scale of the natural vegetation mosaic found in these habitats.The various categories used are also unique among grassland remote sensing studies in their close (although not direct) relevance for Natura 2000: our results are directly comparable with field-based maps created by botanists.We expect that the use of LIDAR datasets in grassland mapping will be tested in other sites in the future, and especially surveys planned from the ground up for grassland vegetation would produce better accuracies.The availability of high-resolution vegetation maps with good agreement to field references will also significantly advance grassland ecology and conservation.

Figure 1 .
Figure 1.Location of the study area in the northwestern part of Hungary with the LIDAR flight pattern and the boundaries of the meadow study sites.

3. 3 . 1 .
Set of 10 Classes "Not vegetation": In our case, non-vegetation surfaces were asphalt roads, building roofs, open water surfaces and bare soil including some field tracks.In case of the studied meadows, open soil reflects various kinds of disturbance: vehicle traffic, animal trampling, and wild boars foraging are some examples.

Figure 2 .
Figure 2. Field photographs showing examples of the category "Fringe".(a) Tall fringes formed by invasive species Fallopia japonica and Impatiens glandulifera.(b) Lower and broader fringe dominated by Stinging nettle (Urtica dioica).

Figure 3 .
Figure 3. Field photographs showing examples of the category "Abandoned".(a) A recently abandoned meadow already showing changes in structure.(b) In the foreground, the abandoned side of the meadow has high coverage of woody species, while the area in the background is still managed.

Figure 4 .
Figure 4. Field photographs showing examples of the category "meadow like".(a) A field track with a high vegetation cover.(b) An open, unmanaged grassland heavily influenced by timber storage and traffic from logging

Figure 5 .
Figure 5. Field photographs showing examples of Natura-2000 related classes.(a) Typical Lowland hay meadow (6510) with well-structured canopy and high species diversity.(b) Typical "dry meadow" patch.Note that the canopy is lower and less structured than for 6510.

Figure 6 .
Figure 6.Field photographs showing examples for patches of the class "wet high".(a) A wet area dominated by Phalaris arundinacea.(b) Characteristic tall wetland vegetation formed by Carex acuta with Iris pseudacorus.

Figure 7 .
Figure 7. Field photograph showing example of the class "Molinia".Note that Molinia caerulea forms a typical tussock shape.

Figure 8 .
Figure 8. True color aerial image showing a studied meadow with ground truth polygons, and the vegetation map resulting from hard-boundary classification of LIDAR data.For areas excluded from classification (forests), the vegetation map is left transparent.Note strictly defined boundaries between classes.

Figure 9 .
Figure 9.True color aerial image showing a studied meadow with ground truth polygons, and the vegetation map resulting from fuzzy visualization of the LIDAR-derived grassland classification.For areas excluded from classification (forests), the vegetation map is left transparent.Note that the overall vegetation pattern is the same as in Figure8, but fine-scale differences and transitions are represented in more detail due to the fuzzy visualization.

Table 1 .
Overview table of classification categories and their properties

Table 2 .
Overview of input datasets and processing steps

Table 4 .
Confusion matrix of 5 vegetation classes and accuracies