Individual tree detection using UAV-lidar and UAV-SfM data: A tutorial for beginners

2.1 UAV-SfM data

For this tutorial, we considered two study sites where point clouds were obtained using different techniques (see Figure 1). In the first study site, point clouds were built from a UAV carrying a RGB sensor through the stereo matching of multiple overlapping aerial images. This method is referred to as the SfM technique, and it helps solve issues associated with the geometry, camera positions, and orientations [23,24,25,26]. For this study site, trees were not georeferenced in the field. This is shown as a low-cost alternative that is useful for multiple applications and for regions where resources for purchasing sensors or extensive field data collection are not feasible, although limitations related to accuracy should be considered [25,34]. The site is approximately 700 ha located at the E.O. Siecke State Forest, East Texas, and is managed by the Texas A&M Forest Service.

Figure 1

The study sites and respective point clouds from the UAV surveys. (a) United States; (b) Texas; (c) Florida; UAV surveyed area of (d) study site 1 and (e) study site 2; example of theSfM point cloud for the (f) study site 1 in Texas and the (g) study site 2 in Florida.

At this site, we selected an area of 11.95 ha for this study, which was dominated by the oldest slash pine (Pinus elliottii Engelm.) plantations. The aerial imagery was acquired in August 2020 using a DJI Mavic Pro quadcopter. The RGB imagery was collected using the Pix4Dcapture flight planning app (https://pix4d.com/) with the specifications described below (see Table 1). Pix4D Mapper was used for initial image processing, point cloud generation, and orthomosaic creation. For the point cloud densification, an image scale of ¼ and an optimal point density were used. The output point cloud was generated in the LAS format. The raster digital surface model (DSM) was created using inverse distance weighting method.

2.2 UAV-lidar data

The second study site was surveyed by a UAV-lidar system and had trees precisely mapped in the field, representing a best-case scenario regarding point cloud generation and field validation. The site is the Ordway Swisher Forest Dynamics Plot (OSFDP) at the Ordway-Swisher Biological Station in Florida. The station is operated as a long-term research facility by the University of Florida and is part of the Global Earth Observation Network (ForestGEO) (https://forestgeo.si.edu/). The plot has an area of 23.04 ha established and mapped from March 2019 to February 2020. All the trees (dbh > 1 cm) were tagged and had their species, dbh, height, status (living or dead), crown light exposure, and position recorded. Tree mapping was done by first locating each one into a 40 by 40 m quadrat and measuring the azimuth (Suunto K-14 sighting compass) and distance with a metal tape to the georeferenced quadrat center point. The quadrat center point geolocation was acquired in a real-time kinematic (RTK) survey using two survey-grade GPS receivers tape to measure antenna height above the monument, a tripod for the base, a bipod for the rover, a range pole, and a Topcon Data Collector. The dominant species were longleaf pine (P. palustris) and turkey oak (Quercus leavis Walter). The plot measurements are publicly available at the ForestGEO website under request (https://forestgeo.si.edu/explore-data/ordway-swisher-termsconditionsrequest-forms) and with agreement to the use terms.

The lidar data were collected in June 2019 with the GatorEye Unmanned Flying Laboratory (http://www.gatoreye.org/). This system is composed of a DJI M 600 Pro hexacopter with a Phoenix Scout Ultra core, which has a STIM300 inertial measurement unit (IMU) coupled with a differential GNSS antenna and integrates a Velodyne Ultra Puck 32c, a 24 MP visual camera, and a Headwall Photonics Nano hyperspectral camera. The lidar dense point clouds (∼288 pts/m²) were acquired from three separate flights 80 m apart, following terrain at 80 m above ground level and at a ground speed of 10 m/s (see Table 2). Flightlines were post-processed using the GatorEye multi-scale post-processing (GMSPP) workflow (v. 229 detailed at http://www.gatoreye.org/).

3.1 Canopy height model (CHM) generation

CHM refers to the distance between the ground level and the topmost point of the objects (which are treetops in our case) under consideration and gives the actual height of objects (Figure 4). Traditionally, CHMs can be generated by (i) subtracting the Digital Terrain Model (DTM) (DTM; which gives a measure of the ground elevation) from the DSM (which represents features on the earth’s surface) or (ii) height normalizing the point cloud using the DTM elevation values and deriving the elevation model from the topmost point cloud returns.

Figure 4

DSM, DTM, and CHM generation.

We used the functions available in the lidR package [35] for CHM generation and one sample plot of 900 m² (30 m × 30 m) at random in each study site as examples for processing. The plot polygons are included in the data sources described in Section 2.1 and 2.2 in “.shp” format that you can download (as per instructions provided in Section 2.3) and load into R Studio.

• Step 1 – Loading the data: Use the readLAS function to load the 3D point clouds. This function can read both “.las” and “.laz” files (lines 6 and 8 of S1). Specifications and standards of the “.las” format can be found in http://www.asprs.org/a/society/committees/standards/LAS_1_4_r13.pdf. For files that are excessively large to be loaded at once, the readLAScatalog function can be used (for example, in case if you are loading data for the entire study area all at once). This function is used as a representation of one or multiple .las files and works with several lidR functions for processing the point clouds without loading it into R (details on the LAScatalog can be found here: https://cran.r-project.org/web/packages/lidR/vignettes/lidR-LAScatalog-class.html).

• Step 2 – Clipping the area of interest: Use the clip_roi function to clip a point cloud. This function allows clipping the point cloud based on a given geometry, such as shapefiles. Clipping can help in the visualization of a region of interest (ROI) or for increasing processing efficiency. We clipped one of the example-plots in each site as the ROI. In this example, we created a 3 m buffer (lines 14 and 15 of S1) around the plot before clipping (lines 18 and 19 of S1) with the purpose of avoiding edge effects in the following processes.

• Step 3 – Classifying ground points: Use the classify_ground function to classify the points that represent the ground. This function changes the “Classification” attribute of points that represent the ground to the value of “2” following the.las file formatting standards (see step 1). For lidar point clouds, there is still an option for using only the last returns in this process as those are most likely the ones from the ground. We can assess the quality of the ground classification by plotting the point cloud by its classification values (Figure 3b; lines 31 and 32 of S1) or directly assess the generated DTM in the next step. If improvements are necessary, the parameters of the algorithm can be changed empirically or based on previous publications in similar types of forests. Furthermore, lidR provides other ground filtering algorithms such as the CSF [36]. The points classified as ground returns are used in the next step for generating the DTM.

• Step 4 – Creating a DTM: Use the grid_terrain function for creating the DTM. This function creates a rasterized surface representing the terrain by interpolating the ground points and has several algorithms available. Herein, we used the k-nearest neighbor with inverse distance weighting (lines 37 and 38 of S1). The cell size of the DTM can be defined as a function of the point density [37], empirically, or based on previous works in similar types of forests. Given the low tree density in most of the area and the high density of the point clouds (especially for the lidar dataset), we defined the cell size to be 0.25 m. The quality and consistency should be assessed in this part on a case-by-case basis as the subsequent process relies on it. For instance, you can use the function plot_dtm3d as a reference to visually assess the consistency of the DTM (Figure 5).

• Step 5 – Height normalization of the point cloud: Use the function normalize_height for normalizing the point cloud elevation values. This algorithm subtracts the DTM elevation (z value) from the elevation of all points. After this process, the lowest returns (i.e., ground returns) in the point cloud are set to ∼0, and the point cloud elevation values (z) will represent the true height of the objects (lines 43 and 44 of S1).

• Step 6 – Creating a CHM: Use the grid_canopy function for creating the CHM. This function creates a rasterized surface using the upper returns of the point cloud and has several algorithm options (refer to lidR package description for more details; https://cran.r-project.org/web/packages/lidR/index.html). Herein, we used the method p2r (point to raster) that holds the highest value for a user-defined voxel and interpolated them into a rasterized surface (lines 48 and 49 of S1).

Figure 5

Processing steps for point cloud height normalization using the SfM and lidar-derived point clouds from the study sites in Texas and Florida (a and b). Starting from the raw point cloud with elevation values (a1 and b1), the ground points are classified (a2 and b2). By subtracting the DTM from the point cloud elevation values (a3 and b3), a normalized point cloud with true height values is obtained (a4 and b4).

3.2 Individual tree detection (ITD)

One of the most well-known, effective, and simplest methods for ITD is the LM algorithm, which is incorporated in this study. Within the LM algorithm, the treetops are associated with the high-intensity LM of the imagery, and we can further apply smoothing techniques and height thresholds to get rid of spread-out tree branches and contorted snags which otherwise might create spurious LM [1,2,28,29,38]. For our study, we used a smoothing window size (SWS) of 5 × 5 pixels as advised for open canopies [1,2,5]. Additionally, a fixed tree window size (TWS) of 7 × 7 and 5 × 5 pixels for the Texas and Florida sites, respectively, was employed as it allows us to define the boundary within which the algorithm has to look for treetops. For more information on optimizing SWS and TWS combinations for enhancing tree accuracy based on canopy density, please refer to [5].

• Step 7 – Smoothening the CHM: We used the function CHMsmoothing to smoothen the CHM (see Figure 6). This is an optional but commonly used process that can improve the treetop detection algorithm results as it helps eliminate spurious LM (e.g., caused by branches) [1,2,28,29,38]. This function uses an image convolution kernel allowing the application of mean, median, or gaussian-based effects. We used the “mean” option as often applied using a 5 × 5 window size. Note that the mean filter changes the pixel values, and it may be needed to verify it when assessing the height of the objects (lines 53 and 54 of S1).

• Step 8 – Detecting the treetops: The function FindTreesCHM, which employs an LM filter, was applied to detect the treetops and retrieve their heights (see Figure 6). The LM are found using a user-defined fixed window. It is also necessary to define a minimum height threshold to exclude non-tree features. We defined a 7 × 7 grid cells window for the point clouds from the treetop detection. Furthermore, a 3 m minimum height threshold was set for avoiding shrubs and dead trees (lines 56 and 57 of S1).

Figure 6

Data processing and CHM creation (a1 and b1); smoothened CHMs (a2 and b2); application of LM filter to detect the treetops (a3 and b3). The CHMs were generated from the SfM (a1–3) and li-dar-derived point clouds (b1–3) obtained in UAV surveys in Texas and Florida sites, respectively.

LM filters are a well-known, widely used, simple, and effective method for ITD. It is important to keep in mind that it is a raster-based method and allows the identification of trees that are in the upper canopy strata. For more complex forest environments such as areas with steep slopes or dense canopies, alternative algorithms – such as TM or VF or deep learning – and approaches or their combinations should be tested [19,39,40,41,42,43,44]. Furthermore, it is worth noting that most of the functions used here have optional parameters that might be different depending on the forest structure and must be tested based on inferences drawn from previous studies [2]. The R packages presented here also have several other interesting features – for filtering (see lidR filter_poi function), processing large datasets (see lidR catalog_apply function), delineating crowns (see rLiDAR forestCAS function), segmenting tree crowns (see lidR segment_trees function), and building 3D forest representations (see rLiDAR LiDARForestStand function) – that we encourage the users to explore.

3.3 Accuracy assessment

For validation purposes, we compared the LM-based ITD results with the manual interpretation of tree counts made from the high-resolution imagery for the UAV-SfM dataset. This is the most common method used for regional-level studies as the acquisition of concurrent data is not economically viable for small-scale landowners. For the UAV-lidar data, we used the position of the trees in the field as a reference for comparison. A total of 10 random plots of 900 m² (30 m × 30 m) in each site were considered for this study. Accuracy metrics calculated include true positive (TP, correct detection), false negative (FN, omission error), false positive (FP, commission error), recall (r), precision (p), and F-score (F); equations are listed below Equations (1)–(3) and for more information on the individual metrics, please refer to [45,46,47]. Herein, we can obtain a measure of trees detected from recall; precision gives a measure of correctly detected trees; F-score provides a measure of the test’s accuracy by taking the harmonic mean of recall and precision.

Overall, the F-scores were approximately 0.8, and omission was higher than commission errors (see Table 3). In the Texas site, 49 of the 64 reference trees were correctly detected. The omission and commission errors in this site were 24.6 and 18.4%, respectively. However, in the Florida site, 70 of the 87 reference trees were correctly detected with an omission error of 19.5% and commission error of 18.3%. It should be borne in mind that the tree detection accuracy is highly dependent on the forest structure, topography, and canopy characteristics. Homogeneous forests, plantations, and woodlands usually will present higher accuracies as the same window size will work in a similar way throughout the area [2,5]. It is highly recommended that the user bases their accuracies on previous works in the same type of forests. To improve the accuracy and F-score, it is possible to vary the smoothening and LMF window sizes as well as check the consistency of the DTM.