Individual Tree Detection in a Eucalyptus Plantation Using Unmanned Aerial Vehicle (UAV)-LiDAR

: The present study addresses the tree counting of a Eucalyptus plantation, the most widely planted hardwood in the world. Unmanned aerial vehicle (UAV) light detection and ranging (LiDAR) was used for the estimation of Eucalyptus trees. LiDAR-based estimation of Eucalyptus is a challenge due to the irregular shape and multiple trunks. To overcome this di ﬃ culty, the layer of the point cloud containing the stems was automatically classiﬁed and extracted according to the height thresholds, and those points were horizontally projected. Two di ﬀ erent procedures were applied on these points. One is based on creating a bu ﬀ er around each single point and combining the overlapping resulting polygons. The other one consists of a two-dimensional raster calculated from a kernel density estimation with an axis-aligned bivariate quartic kernel. Results were assessed against the manual interpretation of the LiDAR point cloud. Both methods yielded a detection rate (DR) of 103.7% and 113.6%, respectively. Results of the application of the local maxima ﬁlter to the canopy height model (CHM) intensely depends on the algorithm and the CHM pixel size. Additionally, the height of each tree was calculated from the CHM. Estimates of tree height produced from the CHM was sensitive to spatial resolution. A resolution of 2.0 m produced a R 2 and a root mean square error (RMSE) of 0.99 m and 0.34 m, respectively. A ﬁner resolution of 0.5 m produced a more accurate height estimation, with a R 2 and a RMSE of 0.99 and 0.44 m, respectively. The quality of the results is a step toward precision forestry in eucalypt plantations.


Introduction
In order to manage the growth and yield of planted forests, several parameters are required. The basic parameters include structural characteristics, the spatial position of stems and crowns, and species identification [1]. The numeric parameters are calculated by introducing, in allometric equations, the variables of interest, above all: tree species identification, tree height, and diameter at breast height (DBH) [2]. Traditional forest inventory is carried out mainly to measure these variables. This consists of field measurement campaigns that are time-consuming and, thus, costly, and the analysis of these results has low applicability in other plantations.
Light detection and ranging (LiDAR) sensors have the potential to overtake the productivity of traditional forest inventory by collecting information with high accuracy, high resolution, and in a short timeframe [3]. LiDAR-acquired data allow for the estimation of stand attributes over large areas at a much lower cost than traditional field-based methodologies. The most common platforms to obtain LiDAR data are airborne platforms, either manned aircrafts or unmanned aerial vehicles (UAVs) [4]. Airborne laser scanning has proven its utility at the aforementioned forest inventory applications as well as for ecological purposes such as fuel assessment, fire prevention, biodiversity assessment, Voxel space detection and delineation (VDD)/LM based on seeded k-means clustering (CDPD) DR: VDD 99%/CDPD 100.6% OE: VDD 6.4%/CDPD 4.6% CE: VDD 7.3%/CDPD 5.5% 1 RMSE: root mean square error; DR: detection rate (number of detected trees divided by the number of reference trees); OE: omission error (false negative); CE: commission error (false positive).
Most of the aforementioned research works on Eucalyptus stands are based on LiDAR sensors aboard manned aircrafts. The high cost of integrated LiDAR sensors into UAVs makes their regular use only possible for very large projects. Therefore, the use of such systems in environmental monitoring and forestry is still greatly limited [46,47]. UAV-LiDAR allows the production of LiDAR data at a fraction of the cost of airplane LiDAR [48]. However, with UAV, it is possible to reach difficult-to-approach areas such as narrow valleys and to acquire higher resolution point clouds [20,48,49]. For these reasons, the weight of laser scanners was reduced to meet the carrying capacity of UAVs [50]. In 2016, Pilarska et al. [51] mentioned that ultralight laser scanners weigh up to 2-3 kg. In 2018, using a laser scanning with a weight of 1.6 kg, Salach et al. [52] obtained a point density of 180 points/m 2 .
As can be seen in Table 1, Gonçalves-Seco et al. [41] obtained a low detection performance in a Eucalyptus plot. They attributed this result to the high density of trees per hectare (up to 3680), the severe defoliation by the snout beetle (Gonipterus sp.), and to the insufficiency of the point density. Barnes et al. [53] agreed with regard to how defoliation can affect the accuracy at tree segmentation. Wallace et al. [45] confirmed that point density significantly influences the successful detection of four-year-old Eucalyptus plantations. When the ITD approach is combined with an ABA, the exact location of each tree can be lost, as in the methodology presented by Shinzato et al. [9]. This work aims to describe two alternative methodologies based on UAV-LiDAR data to perform ITD on Eucalyptus plantations, which differ from the commonly-used CHM-LM approaches. The accuracy in tree detection was evaluated and compared with the CHM-LM algorithm results. The height of detected trees was automatically estimated and evaluated. The study area was an Eucalyptus plantation in Galicia (Spain). To date, the Galician efficiency in forest inventory remains underdeveloped, according to the public authorities [54]. The purpose of this paper was to provide tools to improve Eucalyptus forest inventory in this region through LiDAR point cloud analysis. Due to the extent of the analyzed Eucalyptus stand, LiDAR was obtained from UAV.

Study Area
The study area is located in Galicia, Spain. Galicia is one of the European regions with higher forest productivity, ranging from 6 to 12 m 3 /ha/year, depending on the particular species [55], and supplies 4.5% of Europe's timber [56]. With only 5.9% of the Spanish area, Galicia provides 54.3% of the national timber cutting [57]. In 2009, 21.9% of the Galician forestland consisted of Eucalyptus spp., primarily globulus and nitens spp. [58,59]. This accounts for circa 60%-70% of Spanish Eucalyptus cuttings [60], and accounted for approximately 50% of the Spanish production of chemical pulp over the period from 2008-2015 [61]. The study was performed in two Eucalyptus plots, the west plot and east plot of 34 ha and 25 ha, respectively ( Figure 1). The elevation ranges from 324 m to 440 m and the slope from 0 • to 69 • . The tree spacing across the whole plantation was 2.0 × 4.0 m. Defoliation by the snout beetle affected both plots. Despite the uniform tree spacing in both plots, a number of paths and forest gaps were present across the scanned area, especially in the west one. This work aims to describe two alternative methodologies based on UAV-LiDAR data to perform ITD on Eucalyptus plantations, which differ from the commonly-used CHM-LM approaches. The accuracy in tree detection was evaluated and compared with the CHM-LM algorithm results. The height of detected trees was automatically estimated and evaluated. The study area was an Eucalyptus plantation in Galicia (Spain). To date, the Galician efficiency in forest inventory remains underdeveloped, according to the public authorities [54]. The purpose of this paper was to provide tools to improve Eucalyptus forest inventory in this region through LiDAR point cloud analysis. Due to the extent of the analyzed Eucalyptus stand, LiDAR was obtained from UAV.

Study Area
The study area is located in Galicia, Spain. Galicia is one of the European regions with higher forest productivity, ranging from 6 to 12 m 3 /ha/year, depending on the particular species [55], and supplies 4.5% of Europe's timber [56]. With only 5.9% of the Spanish area, Galicia provides 54.3% of the national timber cutting [57]. In 2009, 21.9% of the Galician forestland consisted of Eucalyptus spp., primarily globulus and nitens spp. [58,59]. This accounts for circa 60%-70% of Spanish Eucalyptus cuttings [60], and accounted for approximately 50% of the Spanish production of chemical pulp over the period from 2008-2015 [61]. The study was performed in two Eucalyptus plots, the west plot and east plot of 34 ha and 25 ha, respectively ( Figure 1). The elevation ranges from 324 m to 440 m and the slope from 0° to 69°. The tree spacing across the whole plantation was 2.0 × 4.0 m. Defoliation by the snout beetle affected both plots. Despite the uniform tree spacing in both plots, a number of paths and forest gaps were present across the scanned area, especially in the west one.

Data Acquisition and Pre-Processing
The aerial survey was conducted in 2017 above both studied plots. The LiDAR data were collected using the Matrice600-PRO multicopter, which was equipped with a VelodynePuck LITE TM laser scanner. The LiDAR sensor had a measurement range of 100 m, an accuracy range up to 3 cm, and generated up to~600,000 points/second at dual return mode across a 360 • horizontal field of view and with a 30 • vertical field of view. Rotation rate was 5 Hz-20 Hz and laser wavelength was 903 nm. The data were acquired by flying at 60 m above ground level with a programmed speed of 7 m/s. Direct georeferencing was applied using a STONEX S10 GNSS receiver. The final position accuracy was 1.5 cm in x-y and 2 cm in z.
The raw LiDAR data were filtered in order to delete the obtained noise. The point cloud resulted in a density of 55.2 points/m 2 in the west plot and 67.0 points/m 2 in the east plot. The next processing steps were performed using the available tools of LASTools software [62]. First, the ground points were identified using lasground, while the height of each point above the ground was computed using lasheight. Therefore, each point of the point cloud contains its X, Y, and normalized Z coordinates. After that process, it was possible to obtain a CHM of the study area. Two CHMs with different resolutions (0.5 m and 2.0 m) were directly created from the normalized point cloud. They were created using lasgrid tool with the option "-highest", which assigns to each pixel the highest Z coordinate among all LiDAR points contained into the corresponding 0.5 m by 0.5 m (or 2.0 m by 2.0 m) cell.

Tree Detection and Tree Height Determination
This work evaluates two procedures to calculate the number of Eucalyptus trees: Method 1 and Method 2. The results were compared to those obtained through the application of local maxima algorithms. Figure 2 presents the workflow followed.

Data Acquisition and Pre-Processing
The aerial survey was conducted in 2017 above both studied plots. The LiDAR data were collected using the Matrice600-PRO multicopter, which was equipped with a VelodynePuck LITE TM laser scanner. The LiDAR sensor had a measurement range of 100 m, an accuracy range up to 3 cm, and generated up to ~600,000 points/second at dual return mode across a 360° horizontal field of view and with a 30° vertical field of view. Rotation rate was 5 Hz-20 Hz and laser wavelength was 903 nm. The data were acquired by flying at 60 m above ground level with a programmed speed of 7 m/s. Direct georeferencing was applied using a STONEX S10 GNSS receiver. The final position accuracy was 1.5 cm in x-y and 2 cm in z.
The raw LiDAR data were filtered in order to delete the obtained noise. The point cloud resulted in a density of 55.2 points/m 2 in the west plot and 67.0 points/m 2 in the east plot. The next processing steps were performed using the available tools of LASTools software [62]. First, the ground points were identified using lasground, while the height of each point above the ground was computed using lasheight. Therefore, each point of the point cloud contains its X, Y, and normalized Z coordinates. After that process, it was possible to obtain a CHM of the study area. Two CHMs with different resolutions (0.5 m and 2.0 m) were directly created from the normalized point cloud. They were created using lasgrid tool with the option "-highest", which assigns to each pixel the highest Z coordinate among all LiDAR points contained into the corresponding 0.5 m by 0.5 m (or 2.0 m by 2.0 m) cell.

Tree Detection and Tree Height Determination
This work evaluates two procedures to calculate the number of Eucalyptus trees: Method 1 and Method 2. The results were compared to those obtained through the application of local maxima algorithms. Figure 2 presents the workflow followed. Local maxima algorithms were applied using two softwares: Fusion [63] and Quantum Geographic Information System (QGIS) [64]. Fusion tool CanopyMaxima identifies local maxima on a CHM using a variable size evaluation window. Equations used to calculate the window size were specified in the Fusion User's manual [65]. Default parameters were used to apply the algorithm to the study case. QGIS uses the local maxima and minimum algorithm [62]. Both algorithms were applied on CHM with different resolutions of 0.5 m and 2.0 m.
Instead of delineation of individual tree crowns, as in the CHM-LM methods, the proposed methods will attempt to delineate individual tree locations by classifying returns from the tree stems. The first step for these approaches is to extract the layer of points that correspond to the stems. To that end, threshold height values were used to classify the point cloud into vertical strata: ground, shrubs, stems, and canopy. Classification was performed on the normalized LiDAR point cloud using the lasheight tool from LASTools software [66]. Height thresholds are determined through point cloud visual inspection of the shrub stratum height and the crown base height for both analyzed plots. The adopted threshold values are shown in Table 2. They differ between plots as the west plot had taller shrubs and trees than the east plot. Local maxima algorithms were applied using two softwares: Fusion [63] and Quantum Geographic Information System (QGIS) [64]. Fusion tool CanopyMaxima identifies local maxima on a CHM using a variable size evaluation window. Equations used to calculate the window size were specified in the Fusion User's manual [65]. Default parameters were used to apply the algorithm to the study case. QGIS uses the local maxima and minimum algorithm [62]. Both algorithms were applied on CHM with different resolutions of 0.5 m and 2.0 m.
Instead of delineation of individual tree crowns, as in the CHM-LM methods, the proposed methods will attempt to delineate individual tree locations by classifying returns from the tree stems. The first step for these approaches is to extract the layer of points that correspond to the stems. To that end, threshold height values were used to classify the point cloud into vertical strata: ground, shrubs, stems, and canopy. Classification was performed on the normalized LiDAR point cloud using the lasheight tool from LASTools software [66]. Height thresholds are determined through point cloud visual inspection of the shrub stratum height and the crown base height for both analyzed plots. The adopted threshold values are shown in Table 2. They differ between plots as the west plot had taller shrubs and trees than the east plot. Afterward, the stem stratum is subset from the point cloud. Finally, the stem strata point cloud is projected in the horizontal plane. Figure 3 illustrates the described process. Figure 3a is an example of the point cloud classification step through a classified vertical profile of a single tree row. Figure 3b contains a front view of the segmented stem stratum. Finally, Figure 3c contains the horizontal projection of the stem stratum points.
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 17 Afterward, the stem stratum is subset from the point cloud. Finally, the stem strata point cloud is projected in the horizontal plane. Figure 3 illustrates the described process. Figure 3a is an example of the point cloud classification step through a classified vertical profile of a single tree row. Figure  3b contains a front view of the segmented stem stratum. Finally, Figure 3c contains the horizontal projection of the stem stratum points.  Method 1 starts by creating a 2D buffer around each point of the projected stem point cloud. The buffer width must be above the X-Y point spacing and below the tree spacing. Several values were tested (0.30 m, 0.50 m, 1.00 m, 1.25 m). As a result, the point cloud was transformed into a cloud of polygons. Then, the buffers that intersected were combined into a single polygon. Figure 4b shows . As a result, the point cloud was transformed into a cloud of polygons. Then, the buffers that intersected were combined into a single polygon. Figure 4b shows an example of the resulting polygons. Then, the cartographic coordinates X, Y of the centroid of the single polygons were obtained. These centroids approximate the geospatial position of every individual tree.
Method 2 is based on considering that a high point density in the horizontally projected point cloud of stems corresponds to an individual tree. The first step is the creation of a two-dimensional raster of the projected stem point cloud calculated from a kernel density estimation with an axis-aligned bivariate quartic kernel. The spatial resolution selected to create the raster was 0.1 m and the kernel radius was 1 m. In order to optimize high density areas, a standardization of the density grid was performed. Figure 4c provides an example of the standardized grid. Finally, a local maxima algorithm was applied to the resulting standardized grid. The result was a point set that corresponded to individual tree geolocation.
Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 17 an example of the resulting polygons. Then, the cartographic coordinates X, Y of the centroid of the single polygons were obtained. These centroids approximate the geospatial position of every individual tree. Method 2 is based on considering that a high point density in the horizontally projected point cloud of stems corresponds to an individual tree. The first step is the creation of a two-dimensional raster of the projected stem point cloud calculated from a kernel density estimation with an axisaligned bivariate quartic kernel. The spatial resolution selected to create the raster was 0.1 m and the kernel radius was 1 m. In order to optimize high density areas, a standardization of the density grid was performed. Figure 4c provides an example of the standardized grid. Finally, a local maxima algorithm was applied to the resulting standardized grid. The result was a point set that corresponded to individual tree geolocation.  Once the tree detection methods were accomplished, the height analysis was performed. The height of each individual was estimated through the CHM. Thus, the center coordinates of each detected tree was intersected with the CHM to obtain the height value of each tree. Finally, all points under 5 m height were deleted, as they might correspond to shrubs. As a result of the detection methods, the position of the detected trees was obtained in the form of a set of points with cartographic coordinates. Thus, the CHM height value for those coordinates was extracted. This procedure was carried out using the CHM with a 0.5 m resolution and the CHM with a 2.0 m resolution.

Assessment of Methods
The evaluation of the accuracy in the tree detection of the analyzed methods was accomplished through visual inspection of the acquired point cloud before processing. Field work was not performed given the reduced accessibility through the analyzed stands due to the height of thorny shrubs (Ulex sp.). Seven tree rows were randomly selected in the study area, and the point cloud was segmented to isolate every row. The number of trees were manually counted in the front view of the Once the tree detection methods were accomplished, the height analysis was performed. The height of each individual was estimated through the CHM. Thus, the center coordinates of each detected tree was intersected with the CHM to obtain the height value of each tree. Finally, all points under 5 m height were deleted, as they might correspond to shrubs. As a result of the detection methods, the position of the detected trees was obtained in the form of a set of points with cartographic coordinates. Thus, the CHM height value for those coordinates was extracted. This procedure was carried out using the CHM with a 0.5 m resolution and the CHM with a 2.0 m resolution.

Assessment of Methods
The evaluation of the accuracy in the tree detection of the analyzed methods was accomplished through visual inspection of the acquired point cloud before processing. Field work was not performed given the reduced accessibility through the analyzed stands due to the height of thorny shrubs (Ulex sp.). Seven tree rows were randomly selected in the study area, and the point cloud was segmented to isolate every row. The number of trees were manually counted in the front view of the point cloud corresponding to every sampled row plantation. This way, 242 trees were identified. Afterward, the DR of each method was calculated as the number of trees detected by the algorithm in the sample rows were divided by the total number of trees that were manually counted in the point cloud (242).
Finally, the verification of the height attribute was performed by measuring the point-to-point distance in the front view of the point cloud of a sample of randomly selected trees. Specifically, the considered sample consisted of 30 trees randomly distributed all over the study area. The measurement was done manually by selecting the apical point and a ground point placed just by the stem. An example of the measuring method is shown in Figure 5. The resulting values were compared to those obtained from the CHM. point cloud corresponding to every sampled row plantation. This way, 242 trees were identified. Afterward, the DR of each method was calculated as the number of trees detected by the algorithm in the sample rows were divided by the total number of trees that were manually counted in the point cloud (242). Finally, the verification of the height attribute was performed by measuring the point-to-point distance in the front view of the point cloud of a sample of randomly selected trees. Specifically, the considered sample consisted of 30 trees randomly distributed all over the study area. The measurement was done manually by selecting the apical point and a ground point placed just by the stem. An example of the measuring method is shown in Figure 5. The resulting values were compared to those obtained from the CHM.

Tree detection
The numerical results of the tree detection and counting and associated predicted density are collected in Table 3. It should be remarked that, with regard to the final buffer value used in Method 1, the best performance was achieved with 0.30 m. With the other tested values, neighboring trees resulted in overlapped polygons. The higher the applied buffer value, the more this overlapping resulted. Table 3. Stem counting and density results on the whole area using the three applied methods. According to the average distance between the trees (2 m) and rows (4 m), the expected density should be 1250 trees/ha. The density provided by the ITD methods is below this value due to the presence of forest gaps (see Figure 6). CHM-LM with a 2.0 m pixel size provided the lowest number of trees, regardless of the algorithm used. An overview of the ITD results obtained by Method 1 and Method 2 is presented in Figure 7.

Tree detection
The numerical results of the tree detection and counting and associated predicted density are collected in Table 3. It should be remarked that, with regard to the final buffer value used in Method 1, the best performance was achieved with 0.30 m. With the other tested values, neighboring trees resulted in overlapped polygons. The higher the applied buffer value, the more this overlapping resulted. According to the average distance between the trees (2 m) and rows (4 m), the expected density should be 1250 trees/ha. The density provided by the ITD methods is below this value due to the presence of forest gaps (see Figure 6). CHM-LM with a 2.0 m pixel size provided the lowest number of trees, regardless of the algorithm used. An overview of the ITD results obtained by Method 1 and Method 2 is presented in Figure 7.   Although the number of detected trees gives an idea of the achieved results, a verification through visual interpretation on the point cloud was performed to gather accuracy metrics. The DR on the sampled rows is presented in Table 4. The DR for the CHM-LM with a 0.5 m pixel size using QGIS was 106.2%. However, Fusion CHM-LM at the same pixel size provided a DR of 21.5%. Both Method 1 and Method 2 overestimated the number of trees and provided a DR of 103.7% and 113.6%, respectively. Thus, false positives surpass false negatives. Table 4. Stem counting on a set of rows of trees.

Method
Number The visual analysis of the results allowed us to interpret error sources for Method 1 and Method 2. Some examples of commission errors are shown in Figure 8. In the case of Figure 8a, a single canopy return was identified as individual trees because the buffer distance was not large enough to reach the overlap with the corresponding tree. In Figure 8b, a group of points caused an extra local maximum in the density grid. These errors may be induced either by the multi-branch structure of Eucalyptus, or by the penetration of canopy in the stem layer. Similarly, Figure 8c,d correspond to false positives due to the incorporation of shrub returns to the stem layer or the presence of fallen snag trees. Again, Method 2 seems to be more prone to these types of errors.
the located individual; (c) standardized density raster on the horizontal projection of the points and the located individual.
Although the number of detected trees gives an idea of the achieved results, a verification through visual interpretation on the point cloud was performed to gather accuracy metrics. The DR on the sampled rows is presented in Table 4. The DR for the CHM-LM with a 0.5 m pixel size using QGIS was 106.2%. However, Fusion CHM-LM at the same pixel size provided a DR of 21.5%. Both Method 1 and Method 2 overestimated the number of trees and provided a DR of 103.7% and 113.6%, respectively. Thus, false positives surpass false negatives. Table 4. Stem counting on a set of rows of trees.

Method
Number The visual analysis of the results allowed us to interpret error sources for Method 1 and Method 2. Some examples of commission errors are shown in Figure 8. In the case of Figure 8a, a single canopy return was identified as individual trees because the buffer distance was not large enough to reach the overlap with the corresponding tree. In Figure 8b, a group of points caused an extra local maximum in the density grid. These errors may be induced either by the multi-branch structure of Eucalyptus, or by the penetration of canopy in the stem layer. Similarly, Figure 8-c,d correspond to false positives due to the incorporation of shrub returns to the stem layer or the presence of fallen snag trees. Again, Method 2 seems to be more prone to these types of errors. Omission errors were also checked through visual inspection. As an example, Figure 9 represents two stems that were detected as a single polygon in Method 1. In this particular case, the threshold values considered for the classification of the vegetation stratum did not allow for the complete separation of shrub and stems. Some points corresponding to shrub vegetation remained in the stem layer. As a result, the buffer distance dissolved two trees into a single polygon. This mechanism makes Method 1 more prone to false negatives. Once the stems were identified, their location coordinates could be obtained. Method 1 and Method 2 generated close coordinates, as can be appreciated in Figure 10. Conversely, the positions provided by the CHM-LM detection differ from them. Omission errors were also checked through visual inspection. As an example, Figure 9 represents two stems that were detected as a single polygon in Method 1. In this particular case, the threshold values considered for the classification of the vegetation stratum did not allow for the complete separation of shrub and stems. Some points corresponding to shrub vegetation remained in the stem layer. As a result, the buffer distance dissolved two trees into a single polygon. This mechanism makes Method 1 more prone to false negatives. Omission errors were also checked through visual inspection. As an example, Figure 9 represents two stems that were detected as a single polygon in Method 1. In this particular case, the threshold values considered for the classification of the vegetation stratum did not allow for the complete separation of shrub and stems. Some points corresponding to shrub vegetation remained in the stem layer. As a result, the buffer distance dissolved two trees into a single polygon. This mechanism makes Method 1 more prone to false negatives. Once the stems were identified, their location coordinates could be obtained. Method 1 and Method 2 generated close coordinates, as can be appreciated in Figure 10. Conversely, the positions provided by the CHM-LM detection differ from them. Once the stems were identified, their location coordinates could be obtained. Method 1 and Method 2 generated close coordinates, as can be appreciated in Figure 10. Conversely, the positions provided by the CHM-LM detection differ from them. Once the stems were identified, their location coordinates could be obtained. Method 1 and Method 2 generated close coordinates, as can be appreciated in Figure 10. Conversely, the positions provided by the CHM-LM detection differ from them.

Tree Height Determination
Finally, for each tree, its location coordinates were used to extract the corresponding height on the CHM. The CHM resolutions of 0.5 m and 2.0 m were also considered at the calculation of tree height. The height estimation was evaluated on a sample of three plantation rows. The ground truth values of the height were estimated directly in the original point cloud. The correlating parameters are collected in Figure 11.

Tree Height Determination
Finally, for each tree, its location coordinates were used to extract the corresponding height on the CHM. The CHM resolutions of 0.5 m and 2.0 m were also considered at the calculation of tree height. The height estimation was evaluated on a sample of three plantation rows. The ground truth values of the height were estimated directly in the original point cloud. The correlating parameters are collected in Figure 11.

Discussion
Method 1 and Method 2 provide accurate estimations of individual tree detection on Eucalyptus stands, independently of defoliation-related problems. The resulting DR with these methods is similar to the DR obtained in other studies that are focused on tree detection in Eucalyptus healthy stands: 92% [13], 97% [42], and 99% [43]. In the case of Wallace et al. [45], the DR resulted in an overestimation of 100.6% using a LM-based on a seeded k-means clustering approach. Method 1 and Method 2 also produced overestimated results (103.7% and 113.6%, respectively). The overview of the misinterpreted returns illustrates that the obtained errors were mostly linked to tall shrubs and snag trees. In addition, the crown base height of small trees was lower than the taller ones. This fact implies that tree branches of small trees may have penetrated into the generated stem layer, resulting in the misidentification of tree branches as single trees. In summary, most of the analyzed errors were related to the threshold values of height that were used to define the stem layer. These might be adapted to the vegetation of every analyzed plot in order to minimize commission errors and omission errors.
Finally, the accuracy of both methods depends on additional key inputs: the buffer distance and the search kernel radius used in generating the raster. Adequate values depend, above all, on two

Discussion
Method 1 and Method 2 provide accurate estimations of individual tree detection on Eucalyptus stands, independently of defoliation-related problems. The resulting DR with these methods is similar to the DR obtained in other studies that are focused on tree detection in Eucalyptus healthy stands: 92% [13], 97% [42], and 99% [43]. In the case of Wallace et al. [45], the DR resulted in an overestimation of 100.6% using a LM-based on a seeded k-means clustering approach. Method 1 and Method 2 also produced overestimated results (103.7% and 113.6%, respectively). The overview of the misinterpreted returns illustrates that the obtained errors were mostly linked to tall shrubs and snag trees. In addition, the crown base height of small trees was lower than the taller ones. This fact implies that tree branches of small trees may have penetrated into the generated stem layer, resulting in the misidentification of tree branches as single trees. In summary, most of the analyzed errors were related to the threshold values of height that were used to define the stem layer. These might be adapted to the vegetation of every analyzed plot in order to minimize commission errors and omission errors.
Finally, the accuracy of both methods depends on additional key inputs: the buffer distance and the search kernel radius used in generating the raster. Adequate values depend, above all, on two variables: tree spacing and point density. The presented results suggest that a higher point cloud density will allow a lower buffer distance and a lower search radius; consequently, the commission error will be minimized. Furthermore, the more spaced the trees are, the easier the tree isolation will be. Consequently, Method 1 is not expected to be valid in stands with no uniform spacing or low density stands.
Obtained results on tree height estimation show that the extraction of this attribute from both the CHM with a 0.5 × 0.5 m pixel size and CHM with a 2.0 × 2.0 m pixel size provided accurate estimations (R 2 0.99 for both of them, and RMSE of 0.44 and 0.34, respectively). This might indicate that the stem is vertically aligned with the apical point of trees, which is reasonable for a high density plantation like the analyzed one. The obtained RMSE values were similar to those obtained in other studies.

Conclusions
This work presents two methods to individually detect the stems of a Eucalyptus plantation by using UAV-LiDAR data. Both methods were based on the segmentation and processing of stem layer returns. The first (Method 1) is based on creating a buffer around stem returns. The second (Method 2) is based on also creating a kernel density raster on stem returns. Both methods allowed us to address an important barrier reported by other works: tree counting on broadleaves [67] stands and defoliated trees [41].
The results indicate that the proposed methods allow a high-accurate estimation of the number of trees to be obtained: there was a detection rate of 103.7% using Method 1 and 113.6% using Method 2. The accuracy of the proposed methods relies on the morphological uniformity of the analyzed stand: straight and unbranched stems, together with uniform age, height, and crown base height of individuals will provide better results. Vertical structure of the stand might also compromise the results. Shrub heights might be under crown base heights in absolute terms all over the stand. Further research could be oriented to assess the transferability of these methods on Eucalyptus stands with different structures and on different tree species. The presented methods may be automated to perform an individualized characterization in terms of height for every tree in the stand.
The presented study addressed stem counting and tree georeferencing in one of the most widely planted hardwoods in the world. The automated extraction of forest structural parameters is required to move toward precision forestry. Additionally, it could help to cover the lack of management tools for Eucalyptus stands reported by Crecente-Campo et al. [40]. In short, the proposed methodology provides the required accuracy to switch from area-based work to individual tree-based management.