The study was conducted in the Oostvaardersplassen nature reserve, a eutrophic wetland in the polder Zuidelijk Flevoland located in the centre of the Netherlands (Figure 1). It is both a Ramsar site (No. 427) and a Natura 2000 site (site code: NL9802054) for breeding and wintering wetland birds. The wetland consists of a marsh (36 km²) and a “dry” border zone with grasslands (20 km²). To create a diverse landscape in the border zone, cattle and horses were introduced in 1983/1984. In 1992, red deer (Cervus elaphus) were introduced to help decrease the ongoing establishment of shrubs and trees at that time. The area is fenced and there are no large predators such as wolves or bears. As the number of ungulates was not controlled, their numbers grew exponentially and in 2010 maximum numbers were reached, leading to very high densities of >200 animals per km² in the border zone. In 2018, the management changed because the landscape was no longer diverse and instead dominated by shortly grazed grasslands, which resulted in a decrease of bird diversity. From 2018 onwards, the ungulate populations were controlled to much lower densities of <75 animals per km². Although all the ungulates can potentially enter the marsh from the border zone, only red deer make use of it. They enter the reedbeds from the grassland areas and look for shrubs and trees to forage on. Our workflow development focuses on the reed-dominated marsh of the study area (Figure 1) which has been subject to red deer trampling and grazing.

We utilized the 3D point clouds collected during the fourth Dutch national ALS flight campaign, known as Actueel Hoogtebestand Nederland (AHN4), conducted between 2020 and 2022 (https://www.ahn.nl). The raw AHN4 point clouds, released in LAZ format (https://rapidlasso.com/laszip/), include x, y, z coordinates, intensity, the number of returns, and the GPS time of data acquisition. The AHN4 point clouds for the study area were captured in April 2020 and have an average point density of approximately 25 points per square meter. The original point clouds were downloaded from the AHN official public portal (https://hub.arcgis.com/maps/esrinl-content::ahn4-download-kaartbladen-1/explore). This included a total of seven original tiles (20DZ1, 20DZ2, 26AN2, 26BN1, 26BN2, 26AZ2, and 26BZ1). Each of these tiles has a spatial coverage of 5  km × 6.25  km . The total number of points of these tiles is 6,243,044,410 with a total volume of 42.3 GB. These original point clouds were used as the input into the presented workflow.

The downloaded AHN4 point clouds (i.e., seven tiles) of the study area were used as input for workflow step 1 (pre-processing). The pre-processing of the point clouds consists of three steps: (i) re-tiling; (ii) clipping, and (iii) outlier removal. The point cloud for input is large in storage space (i.e., the seven original tiles amount to ∼42.3 GB tile volume in LAZ format), which makes manipulating the data difficult. We therefore implemented a spatial extent-based re-tiling approach, which segments the original input point clouds ( 5  km × 6.25  km tiles ) into 50   m × 50   m tiles. The generated tiles were then clipped using the polygon shapefiles of the marsh area. For the clipping, we implemented the winding number algorithm which parses points within polygons (Kumar and Mallikarjun, 2018).

Statistical outlier removal (SOR) was implemented to remove noise, e.g., points from water surface reflection and sparse higher vegetation points (Vosselman and Maas, 2010), using the method proposed in Rusu et al. (2008). Suppose P is a set of 3D points, and for each query point p q u e r y ∈ P , d ¯ is the mean distance of a query point to its k nearest neighbors. For all points in P , the mean distance and standard deviation of the distances of their k nearest neighbors are then determined. Only those points are kept which have distances that are close to the mean distance of the closest neighbours, using Equation (1). P k = p q ∈ P   μ k − α σ k ≤ d ¯ ≤ μ k + α σ k (1) Here, α is a density threshold coefficient, and μ k and σ k are the mean and standard deviation of the distance from a query point to its k closest neighbors. P k is the point set that is kept, i.e., after removing the outliers.

In summary, the raw ALS point clouds are first re-tiled into smaller spatial units to manage the large file sizes. Then, the tiles are clipped to the marsh boundary using polygon shapefiles. Finally, outliers are removed using the SOR method, which filters points based on their local neighbourhood distance distribution.

In workflow step 2 (near-terrain filtering), the reedbed points are removed and only the terrain and near-terrain points (i.e., ground and low vegetation points) within the study area are kept. A virtual 3D grid layer is used in this filtering approach which applies a horizontal boundary of the 3D point cloud with a user defined grid size S i z e g (Figure 3).

To construct a 3D grid layer from a given 3D point cloud, the first step is to determine the 3D bounding box that encloses the entire point cloud. This bounding box is then discretized into a structured array of 3D grid cells by dividing it along the x , y , z axes. The number of divisions (bins) in each direction is determined based on a user-defined grid size. The origin of this grid, which is referred to as the “lower left point” of the bounding box, is denoted as P min , which represents the minimum coordinate values in each axis direction. Each point in the point cloud is assigned a 3D grid index by comparing its coordinate P i , where i ∈ x , y , z with the corresponding coordinate of P min i using Equation 2. n i = P i − P min i S i z e g i   i ∈ x , y , z (2) Here, n i denotes the grid index of a point along axis i , and S i z e g i is the user-defined grid size in the respective direction.

Grids that contain no points are marked as empty and excluded from further processing. Figure 4 illustrates this process of subdividing the bounding box into grid cells. Initially, the grid size is chosen to be relatively coarse to facilitate efficient organization of the raw point cloud. In later stages of the workflow, particularly during iterative filtering, the grid size is progressively reduced to extract finer terrain detail.

After the creation of the 3D grids, only the bottom layer of the grids will be used for the height interpolation and iterative subdivision. The elevation of each grid in the bottom layer is interpolated by the points inside the grid using Equation (3). E g = ∑ E p L 2 − D g ∑ L 2 − D g (3) Here, E g is the interpolated elevation of a grid at the geometric center, E p is the elevation of point p , L is the size of the grid and D g is the distance from point p to the geometric center of the grid.

The initial terrain surface is defined as the geometric centres of a set of coarse grid cells covering the area of interest. This coarse grid is constructed by dividing the bounding box of the point clouds at a coarse resolution level. The geometric centre of each cell is calculated as the midpoint of its horizontal extent X , Y , formatting a regularly spaced point set. In each iteration, the bottom layer is subdivided into finer cells, and the centres of these finer cells generate an updated virtual terrain surface for filtering (near-)terrain points (Figure 4).

For low vegetation like herbs and sedges, the terrain points and vegetation points might be difficult to separate, especially in places where trails are located. To prevent segmenting vegetation points as terrain points, a local terrain slope is introduced. Rather than using the geometric definition in Sithole and George (2004), the interpolated elevation of grids in each subdivision is used to compute a compensated local slope using Equation 4. S p = H p − H g x p − x g 2 + y p − y g 2 (4)

Here, S p is the local slope value, H p is the elevation value of point p , H g is the interpolated elevation value of the grid in which point p is in, x p , y p are the horizontal coordinates of point p , x g and y g are the geometric center of the grid. The procedure for filtering terrain and near-terrain points in an iterative approach is described in Algorithm 1.

Algorithm 1

Iterative terrain points filtering from ALS point clouds.

Input: ALS points {P}, initial and minimum grid size S i z e i n i t , S i z e min , height threshold T h

In summary, a coarse 3D grid is created over the bounding box of the point clouds and elevation values at each grid cells are interpolated using local point distributions. Through iterative grid subdivision and terrain filtering, vegetation points are separated from terrain and near-terrain points based on height thresholds.

In workflow step 3 (DTM reconstruction), the obtained terrain and near-terrain points are used to generate DTMs. We apply the Inverse Distance Weight (IDW) interpolation to estimate the elevation of each vertex of the DTM (Figure 5).

As shown in Figure 5a, the points in red are the obtained terrain and near-terrain points using the method described above (Section 3.2). The blue grids are the DTM with a user-defined resolution. For a specific vertex P of the DTM (orange point), its elevation is estimated by all red points within the neighborhood radius r using Equation 5. H p = ∑ H i d i ∑ 1 d i (5) Here, H p is the estimated elevation of vertex P on the DTM, H i is the elevation of i -th point within radius r , and d i is its distance to vertex P . The generated DTM will then be used for workflow step 4 (trail extraction).

In summary, the filtered terrain and near-terrain points are used to construct DTMs using IDW interpolation. Elevations at DTM vertices are computed using nearby points within a defined radius.

In workflow step 4 (trail extraction), the Laplacian smoothing method is first employed to smooth the DTM mesh (Häufel, Böge, and Bulatov, 2020; Nealen et al., 2006). As shown in Figure 5b, point P t is one vertex on the DTM (orange point), and points P k , P l , P m and P n are its neighbouring vertices (green points). The smoothing procedure is defined by Equation 6. P i t + 1 = P i t + λ ∆ P i t (6) Here, P i t + 1 is the smoothed point of vertex P i t + 1 after t iteration, λ is a coefficient and Δ P i t is defined by Equation 7. ∆ P i t = 1 N i ∑ j ∈ N i P j − P i (7) Here, N i is the number of neighboring points used for smoothing P_i is the query point and P_j is one of the neighbouring points. As illustrated in Figure 5c, the Laplacian smoothing is applied iteratively on the smoothed DTM mesh. D T M t is a profile of the DTM in the previous iteration and D T M t + 1 is the DTM of the next iteration. Note that only the elevation values z are used for smoothing (Figure 5c). Thus, the horizontal distances between the DTM vertices remain unchanged.

For classifying trail grid cells, the smoothed DTM is then used as input. As illustrated in Figure 5c, D T M t is a profile of a trail at iteration t and D T M t + 1 is the smoothed profile at the consecutive iteration using Laplacian smoothing. The vertices in the shallow area of the trail were smoothed such that their elevation (z) is higher compared to that of the previous iterations. For example, the elevation of point P i t + 1 is higher than that of point P i t . Thus, subtracting the elevation of the two DTMs results in a map of elevation difference where areas of trails show the largest elevation changes. Equation 8 defines the elevation changes of DTMs from two consecutive Laplacian smoothing iterations. Δ H = D T M t − D T M t + 1 (8) Here, Δ H is the calculated height difference. The grid cells of trail areas can then be determined by Equation 9. P t r a i l s = p k ∈ P     Δ H ≤ μ Δ H − κ σ Δ H (9) Here, P t r a i l s are the points in the trail area, p k are the trail points, P is the entire point sets, μ Δ H and σ Δ H are the mean and standard deviation of Δ H , and κ is the coefficient. The trail extraction step thus results in each grid cell being classified as “Trail” or “Non-Trail” (Figure 2h).

The obtained “Trail” points are then cleaned using the SOR method described in Section 3.1. After that, a spatial clustering algorithm is applied to identify patches (i.e., clusters) of connected points based on the results from the Laplacian smoothing of the DTM (Cohen and Deschamps, 2001). This algorithm performs spatial clustering of a given set of 3D points P = P 1 , P 2 , … , P N based on a connectivity criterion defined by a radius r . The goal is to partition P into clusters C 1 , C 2 , … , C k such that for any two points p i and p j , they belong to the same cluster if there exists a path of points p i , p n 1 , p n 2 , … , p j where each consecutive pair of points p a and p b satisfies: d p a , p b ≤ r , here d * is the Euclidean distance computation function. The algorithm begins by constructing a KD-tree for efficient nearest neighbor searching. For each point p i , a radius search retrieves the set of neighors N i , such that N i = p j ∈ P     p i − p j 2 2 ≤ r 2 . Then a union-find data structure also known as Disjoint Set Union (DSU) is used to efficiently merge points into connected components, where each point p i initially represents its own cluster. When p i finds a neighbor p j with j > i , then a union operation is performed using Equation 10. u n i o n i , j → f i n d i = f i n d j (10) After processing all points, clusters are identified by extracting the connected components where all points sharing the same root representative in the DSU form a cluster. This results in a partition of P into k clusters C 1 , C 2 , … , C k using Equation 11. C i = p j   f i n d j = r i } , ∀   j ∈ P (11) The clusters are then used to extract linear features and to refine the trail extraction results. To achieve this, a 3D structure tensor S ∈ R 3 × 3 is calculated at the neighborhood and 3D sparse tensor voting is applied (Jörgens and Moreno, 2015) to enhance the linear feature extraction and to remove small, non-linear and isolated clusters. The dimension of a neighborhood is determined as follows: Suppose p i = x i , y i , z i T are the coordinates of a point p i , and T denotes the transpose of the matrix. Then the centroid of all the neighborhood points is determined by Equation 12. p ¯ = 1 n ∑ i = 1 n p i (12) Here, n is the number of points in the neighbourhood. The 3D structure tensor M of those points is then defined by Equation 13. M = 1 n Q T Q (13) Here, Q = p 1 − p ¯ , p 2 − p ¯ , … , p n − p ¯ T . M is a symmetric matrix and can be decomposed as M = R I R T , R is a rotation matrix and I is a diagonal positive definite matrix in this case. The elements of I are the eigenvalues of matrix M . The three eigenvalues are positive, can be denoted by λ 1 , λ 2 , λ 3 , and are sorted such that λ 1 > λ 2 > λ 3 . The corresponding eigenvectors are v 1 , v 2 , v 3 , respectively. (Weinmann et al., 2015). Thus, Equation 13 can be rewritten to Equation 14 as follows: M = ∑ i = 1 3 λ i   v i v i T (14) The linearity, planarity and sphericity of the structure tensor are defined as follows (Equation 15): L = λ 1 − λ 2 λ 1 + λ 2 + λ 3 v 1 v 1 T P = λ 2 − λ 3 λ 1 + λ 2 + λ 3 v 1 v 1 T + v 2 v 2 T S = λ 3 λ 1 + λ 2 + λ 3 v 1 v 1 T + v 2 v 2 T + v 3 v 3 T (15) The sparse 3D structure tensor defined in Equation 14 can then be reformulated as Equation 16: M = L + P + S (16) Based on the above decomposition, sparse 3D structure tensor voting is applied on the obtained point clusters (Medioni et al., 2000). In this workflow, the linear voting field is computed as given in Figure 6.

The direction of the voting tensor O to position P is obtained by rotating the voting tensor direction with angle 2 α , as denoted by Figure 6a. Here, θ is the angle between vector O X → and O P → . Equation 17 is used to define the weighting function which reduces the strength of the vote with the arc length O P ˘ and the curvature k . W s v = l 2 2 σ 2 − k 2   α ≤ π 4 0   o t h e r w i s e (17) Here, l = v → 2 α sin ⁡ 2 α , k = 2 ⁡ sin ⁡ 2 α v → , and σ is a scale parameter which is defined by the radius that is used for computing the local structure tensor within its neighborhood. The voting tensor propagates its saliency to the neighboring locations and aggregates the saliency that associated with the 3D structure tensor.

In summary, a Laplacian smoothing procedure is iteratively applied to the generated DTMs to capture subtle depressions that are indicative of “trails.” Then, trail locations are identified by computing elevation differences between two consecutively smoothed DTMs. Grid cells with high residuals are classified as “trail.” A spatial clustering algorithm is then applied to group connected trail cells, and sparse 3D structure tensor voting is applied to refine the output by enhancing linear features and removing noisy, non-linear clusters.

To validate the accuracy of the workflow for extracting ungulate trails, we randomly placed a total of 100 plots of 30   m × 30   m size in our study area. All 100 sample plots were not located in water areas nor at the water-reedbed boundary to avoid incomplete plots. Half of the plots were randomly placed in those reedbed areas of the marsh in which only red deer are grazing (red colour in Figure 7), and the other half where both red deer and graylag geese (Anser) are grazing (yellow colour in Figure 7). Areas in which only geese are grazing were not considered (green colour in Figure 7). Besides red deer, graylag geese are also grazing reedbeds, but they mainly feed on young reed shoots that occur near open water (Bakker et al., 2018). To test whether goose grazing has an influence on the trail extraction results in areas where both red deer and geese co-occur, we performed the method validation on areas of the reedbed that have either been subject only to red deer activity or to both red deer and graylag geese grazing. These areas are mapped yearly as part of the breeding bird monitoring in the marsh (Beemster et al., 2020).

The accuracy of trail extraction was evaluated by comparing the results from the presented workflow with the results from manually created ground truths. Manual labelling of the point clouds was done by using the segment tool of the CloudCompare software (https://cloudcompare.org/). The 3D point cloud data were first loaded and then segmented using the DB panel of CloudCompare. Then, the segmentation mode was activated (using the Scissors icon in the toolbar), and a closed polygon was drawn around the desired points by clicking sequentially. Once the selection was complete, the class labelling button was used to classify the selected “trail” points as class “1” and the rest of the points as “non-trail” points (classified as “0”). Areas of “non-trail” depict higher elevations compared to areas of “trail.” After manually labelling the discrete points, a DTM of 10 cm resolution was generated for each plot using the same method as described in Section 3.2 (DTM reconstruction). The class of each grid point was assigned to that of the closest point from the manually created ground truth. A total of 9,000,000 grid points were finally labelled as ground truth across the 100 sample plots. For each plot, the overall accuracy A c c was used as the evaluation metric, which is defined in Equation 18. A c c = T P + T N T N + T N + F P + F N (18)

Here, T P (true positives) are the correctly extracted trail cells of 10 cm size, F P (false positives) represent non-trail cells that were wrongly extracted as trails, F N (false negatives) are trails cells that are incorrectly segmented as non-trails, and T N (true negatives) are the correctly segmented non-trail cells. The overall accuracy was calculated for each of the 100 sample plots.

The presented workflow uses seventeen parameters (Table 1). The specific settings of the parameters for the application in our study area are given in Table 1. The choice of these specific parameter values was informed by sensitivity analyses or constrained by the specific characteristics of the point cloud data from our study area.

A total of twelve parameters were tested with a sensitivity analysis (Table 2). For each of these twelve parameters, ten experiments with different parameter values were conducted (Table 2). All other parameters were kept constant (as given in Table 1). Five of the parameters ( S i z e , k , α , M a x i m u m   g r i d   s i z e , and R e s o l u t i o n ) were not subject to a sensitivity analysis for the following reasons. S i z e , M a x i m u m   g r i d   s i z e , k , and R e s o l u t i o n are the parameters which influence the computation time rather than the output accuracy of the presented workflow. Their values were chosen to make the workflow computationally efficient. The parameter α was kept constant because the used point clouds were all collected with the same flight campaign, i.e., the overall quality of the data in our study area is similar.

For the near-terrain filtering (workflow step 2), the relevant parameters are M a x i m u m   g r i d   s i z e , M i n i m u m   g r i d   s i z e , and H e i g h t   t h r e s h o l d (Table 1). To evaluate the accuracy of this workflow step when changing the parameters (Table 2), we manually labelled all vegetation and near-terrain points in the 50 sample plots which were exclusively grazed by red deer (see example in Supplementary Figure A1). The manually labelled point clouds where then compared to the automatically segmented point clouds from the workflow using different parameters settings for the M i n i m u m   g r i d   s i z e (Supplementary Figure A2) and for the H e i g h t   t h r e s h o l d (Supplementary Figure A3).

For the DTM reconstruction and trail extraction (workflow steps 3 and 4), the accuracy was evaluated based on the manually created trail and non-trail ground truth of the 50 sample plots that are placed in the area that is exclusively grazed by red deer.

Filtering near-terrain and vegetation points is a fundamental step in the proposed workflow to ensure accurate and efficient extraction of ungulate trails. The filtering process begins by taking cleaned raw ALS point clouds as input. By structuring this input into grids of multiple scales, the virtual 3D layer is iteratively refined until near-terrain points are accurately identified. Several existing filtering methods are available for segmenting terrain and vegetation points from raw ALS point clouds (Wan et al., 2018). For example, a mathematical morphology-based filtering method, which uses a moving window to perform sequential opening and closing operations, was initially proposed to separate ground and non-ground points (Lindenberger, 1993). However, this method is sensitive to fine-scale variation in terrain and vegetation points and can thus influence applications such as the extraction of ungulate trails in marsh vegetation. Applying this method to filter point clouds of marsh vegetation would result in the misclassification of vegetation points as near-terrain points, thereby introducing additional errors in the subsequent DTM reconstruction and trail extraction processes. Other algorithms have attempted to improve the overall accuracy of segmenting ground and non-ground points by incorporating additional parameters, such as slope and height threshold adaptation functions (Vosselman, 2000; Sithole, 2001; Hu et al., 2014). An example is an easy-to-use algorithm with three parameters which was proposed to simulate the terrain as a piece of cloth (Zhang et al., 2016). However, this method takes predefined fixed resolutions and is thus not very flexible. Moreover, the marsh vegetation that we studied has a relatively flat terrain, and introducing more parameters not only increases the computational effort but also provides limited improvement in overall accuracy. Some filtering algorithms first identify a proportion of ground points and then progressively expand until all points are processed. For instance, an irregular triangulated network (TIN) can be used to simulate the terrain surface, and points within a predefined distance threshold are then progressively added to the terrain points (Axelsson, 2000; Sohn and Downman, 2002). However, constructing a TIN is less efficient as it costs substantial computational time and memory. Our method using iterative 3D grid filtering is thus more efficient than other filtering methods, especially regarding two aspects. First, it only uses a proportion of the points that are within their corresponding 3D grids rather than the entire point cloud. Second, it does not require to construct a TIN and hence reduces computation time. The method we implemented for filtering of near-terrain points is therefore the most appropriate solution for our application.

We constructed a DTM with the obtained terrain points of the marsh habitat using an IDW interpolation method with two parameters, i.e., the resolution of the desired DTM and the radius by which the height of the DTM grid is interpolated. We predefined those two parameters based on the average point density of the used point clouds, which is 20–30 points per square meters for AHN4. Since the elevational changes in the study area are relatively small and the DTM has a high resolution (10 cm), we did not have to take additional parameters such as the local slope and roughness into account. However, if a point cloud has a lower point density compared to the desired resolution of the DTM, holes and gaps could be generated in the DTM. Other algorithms should then be considered, e.g., those that allow hole filling and more complexed elevation interpolation. An example is the thin plate spline (TPS) interpolation which combines a hierarchical multiresolution approach to construct the DTM (Mongus and Žalik, 2008). However, the accuracy of the resulting DTM can be compromised, e.g., the algorithm needs further weight adaptation if there are strong topographical changes in a study area. Another choice could be the algorithm proposed by Maguya et al. (2013) for complex terrain. This adaptive algorithm adjusts to terrain complexity by partitioning the input data and applying linear, quadratic, or cubic spline interpolation based on the characteristics of the terrain. This adaptive approach can be particularly effective in steep and forested terrains. However, the complexity of the adaptive algorithm results in higher computational demands, which can slow down the processing for large datasets. Users who apply our workflow in other study areas should therefore be aware that the steepness and complexity of the terrain can influence the DTM reconstruction.

The trails created by red deer in our study area have shallow surrounding elevations and resemble the geometric shape of valleys and canals. Common methods for extracting such features from 2D grids involve the watershed delineation, which is widely used in hydrological modeling and management (Lai et al., 2016). For instance, Tarboton (1997) introduced a watershed delineation method for determining water flow directions and upslope areas in grid-based DTMs. This method calculates the flow direction as a single angle, representing the steepest downward slope from each DTM grid cell. By assigning flow between two downslope cells, this approach improves the traditional D8 (eight directions) method by reducing grid bias and minimizing dispersion. However, the method struggles to extract valleys in flat areas and water pits. By connecting local sinks using the accumulated water inflow and direction, Al-Muqdadi and Merkel (2011) addressed this challenge of watershed delineation in flat and arid terrains, where traditional methods often struggle due to low relief and poorly defined drainage patterns. However, the accuracy of their watershed delineation method heavily depends on the resolution and quality of the DTM data and involves a high computational cost. Moreover, applying their watershed delineation to the 10 cm resolution of the generated DTM in our study area would lead to a huge computational cost. Zhang et al. (2020) proposed a watershed merging algorithm designed to resolve the issue of fragmented watersheds caused by traditional delineation methods. This method merges smaller, disconnected watersheds into cohesive units based on the channel network, resulting in more realistic hydrological models. The merging process, however, is computationally expensive, especially when applied to large datasets with numerous small watersheds.

Our proposed trail extraction method uses the height residuals of two differentiated DTMs with a Laplacian smoothing technique, which has lower computational complexity and only two parameters to fine-tune the results. This potentially improves the applicability of the method and reduces computational costs but occasionally could also lead to lower accuracies. Several factors may contribute to reduced trail extraction accuracy in certain areas. First, trails with variable width or low depth may not produce sufficient elevation residuals to surpass the detection threshold, especially when the terrain is extremely flat. Second, dense or homogeneous vegetation can mask subtle terrain depressions, leading to non-classified trail grid cells. Third, variations in point cloud density or noise, such as water surface reflections, may affect interpolation quality and influence the classification of grid cells. Also, this will most likely happen if the coefficient κ is not optimal for the plots, e.g., when a large proportion of the plot area is covered by trails. We further refined the trail extraction results by employing sparse 3D structure tensor voting. The 3D stick tensor voting that we applied is a powerful technique for enhancing linear features in 3D point clouds (Laidlaw and Vilanova, 2012). It propagates local orientation information to preserve geometric continuity and emphasizes structures such as ridges, edges, and skeleton-like formations. To our knowledge, we applied this technique for the first time to extract linear ecological features. Most current applications are for medical image processing and man-made object extraction (Liu et al., 2012; Moreno et al., 2012; Soni et al., 2021). Our application could be extended to the extraction and refinement of other ecological structures such as tree lines, hedges, and stone walls. The application of sparse 3D stick tensor voting to refine the trail extraction results not only removed outliers of the obtained trail points but also allowed combining points that were located between two trail segments.

The results from our trail extraction workflow offer insights into the habitat use of a large mammalian herbivore, in this case red deer, in vast, inaccessible and remote habitats, such as reedbeds. The results suggest that 11.6% of the marsh area is covered by red deer trails. Based on the output of the workflow, the trail density of red deer and the resulting reedbed fragmentation can be mapped using a very high (10 cm) resolution trail map. This allows ecologists and land managers to study the effects of ungulates on other species (e.g., breeding birds) and to adjust the management of large herbivores in case of negative effects. We found the highest density of trails at the edge with grasslands from which red deer are entering the reedbeds. Our study also suggests that it can be important to know which other herbivores are present in the area and affect the vegetation. In our study area, greylag geese also affect the reed vegetation as they forage on the reed leaves during moult. While doing so, they create paths starting from the open water into the reedbeds and eating the reed vegetation away alongside these paths, creating large patches of up to several ha with very low homogeneous reed vegetation that borders the water. However, the accuracy of our workflow was only slightly reduced in areas where both geese and red deer are grazing, suggesting that the effect of geese on red deer trail extraction was minor in our study area. Nevertheless, applications in other study areas should be aware that the effects of other herbivores should be taken into account.

For studying the effects of ungulate trails on the development of other species such as breeding birds, trails should be monitored over time. In the study area, it is known from annual breeding bird monitoring that some breeding bird species, e.g., the Savi’s warbler (Locustella luscinioides), appear to be in decline in areas where red deer trails are more abundant. Furthermore, visual inspection of annual aerial photographs indicate that the density of paths increased as the red deer population grew from 50 individuals in 1992 to over 3,500 individuals in 2017. After 2018, as the red deer population decreased to 900 in 2024, a visual decrease in the number of trails can be seen on the aerial photographs. However, a detailed quantitative and spatial assessment of reedbed fragmentation caused by red deer, and its effects on breeding birds, is still lacking in the study area. For effective management, it is important to understand the extent to which different herbivore densities have negative or positive impacts on nature conservation goals, such as the conservation of certain bird species. This knowledge is essential so that timely action can be taken when biodiversity begins to decline as a result of intensive grazing by ungulates. Remote sensing applications with ALS data, like the workflow we have developed here, can thus provide useful tools for habitat condition monitoring (Kissling et al., 2024).

While the workflow presented in this study has a high accuracy in the study area and shows promising results to inform habitat management, there are a few limitations that need to be considered.

1. Transferability. This workflow has so far been developed and applied only in one study area, the Oostvaardersplassen nature reserve in the Netherlands, and with one ALS dataset (AHN4). Its transferability should be tested in other wetland areas with reedbeds, for different species of ungulates, in a variety of habitats (e.g., outside wetlands), and for ALS datasets from different time periods and with different characteristics. Certain parameter settings, such as the M i n i m u m   g r i d   s i z e and the H e i g h t   t h r e s h o l d , are likely to vary depending on the used LiDAR point cloud data and the type of vegetation. K e r n e l   s i z e might also be related to the specific behaviour of the ungulate and its impact on vegetation structure. R t e n s o r should be related to the width of the trails to ensure the linearity of a trail point in its neighbourhood can be characterized. The parameters were optimized based on the manually created ground truth using the point clouds of the study area. Those parameters probably must be re-optimized when applying the presented workflow to other study areas or datasets.

Addressing these limitations is essential for enhancing the robustness, transferability, and applicability of the presented workflow. Our trail extraction method could also be adapted to extract other linear features such as stone walls, hedges, and tree lines. This would make it a versatile tool for ecological applications and habitat management across diverse environments, using 3D point clouds from airborne laser scanning as an input.