Use of Unoccupied Aerial Systems to Characterize Woody Vegetation across Silvopastoral Systems in Ecuador

: The trees in pastures are recognized for the beneﬁts they provide to livestock, farmers, and the environment; nevertheless, their study has been restricted to small areas, making it difﬁcult to upscale this information to national levels. For tropical developing countries, it is particularly important to understand the contribution of these systems to national carbon budgets. However, the costs associated with performing ﬁeld measurements might limit the acquisition of this information. The use of unoccupied aerial systems (UAS) for ecological surveys has proved useful for collecting information at larger scales and with signiﬁcantly lower costs. This study proposes a methodology that integrates ﬁeld and UAS surveys to study trees on pasture areas across different terrain conditions. Our overall objective was to test the suitability of UAS surveys to the estimation of aboveground biomass (AGB), relying mainly on open-source software. The tree heights and crown diameters were measured on 0.1-hectare circular plots installed on pasture areas on livestock farms in the Amazon and Coastal regions in Ecuador. An UAS survey was performed on 1-hectare plots containing the circular plots. Field measurements were compared against canopy-height model values and biomass estimates using the two sources of information. Our results demonstrate that UAS surveys can be useful for identifying tree spatial arrangements and provide good estimates of tree height (RMSE values ranged from 0.01 to 3.53 m), crown diameter (RMSE values ranged from 0.04 to 4.47 m), and tree density (density differences ranging from 21.5 to 64.3%), which have a direct impact on biomass estimates. The differences in biomass estimates between the UAS and the ﬁeld-measured values ranged from 25 to 75%, depending on site characteristics, such as slope and tree coverage. The results suggest that UASs are reliable and feasible tools with which to study tree characteristics on pastures, covering larger areas than ﬁeld methods


Introduction
Pasture is one of the most important land uses globally, occupying approximately 25% of total ice-free land and 67% of agricultural land [1]. Pastures host about 1.43 billion cattle and 1.87 billion sheep and goats [2], supporting the livelihoods of millions of people around the world [3]; at the same time, they are major contributors of greenhouse gas emissions to the atmosphere [4][5][6]. When trees and shrubs are present in pasture areas, these ecosystems are classified as silvopastoral systems (SPS) [7], which are known for having a positive effect on the farm economy [8], biodiversity [9], and carbon stocks [10][11][12].
that have not been studied enough, such as spatial arrangements at larger spatial extents. While more sophisticated hardware may provide more accurate results, we used inexpensive and readily available hardware to test whether the method could be used in developing countries, where resources are more limited. In this study, we propose an approach that combines field-level techniques (plot measurements) with remote-sensing information (UAS photogrammetry) to study SPS characteristics. The overall aim of the study was to assess the suitability of UAS surveys to estimate aboveground biomass on SPS. To achieve this, we: (1) compared the best parameter combination for point-cloud processing and classification, by evaluating the canopy-height-model differences against field measurements; (2) assessed the performance of individual tree detection derived from canopy-height models vs. visual identification in ortho-mosaic images; and (3) assessed the tree spatial arrangements and distribution of SPS.

Study Sites
This research was conducted between May and December 2018 on 101 livestock farms in Ecuador, as part of the Climate Smart Livestock project [34]. In total, 43 farms from the Amazon region and 58 in the Coastal region of Ecuador were surveyed. For this study, 5 farms per region were randomly selected. Farms were distributed in an altitude range from 8 to 426 m asl and 451 to 2000 m asl, in the coastal (C1-C5) and Amazon (A1-A5) regions, respectively ( Figure 1). The main activity on the studied farms was cattlelivestock production, but most of the farms also dedicated land to other economic activities, such as perennial crops, annual crops, other livestock, or conservation. To achieve the proposed aims of this study, we divided the work into three stages ( Figure 2).First, combining traditional forest field measurements with UAS surveys, we collected tree characteristic information (a). Secondly, we transformed the 2D data (individual photos) into landscape information in 2D (orthomosaic) and 3D (digital surface model (DSM), point cloud), which is explained in Section b in Figure 2. Finally, To achieve the proposed aims of this study, we divided the work into three stages ( Figure 2).First, combining traditional forest field measurements with UAS surveys, we collected tree characteristic information (a). Secondly, we transformed the 2D data (individual photos) into landscape information in 2D (orthomosaic) and 3D (digital surface model (DSM), point cloud), which is explained in Section b in Figure 2. Finally, we used the products of Section b to estimate tree characteristics and compare them against measured values (Figure 2c). we used the products of Section b to estimate tree characteristics and compare them against measured values (Figure 2c).

Plot Delimitation and Tree Inventory
Farmers were asked to identify grazing areas that were representative of their farms [35]. Once the areas were identified, two 1-hectare square plots were installed. Plot corners were marked with a plastic sheet to make them visible from the UAS survey, and geographical coordinates were collected with a Garmin GPS50 GPS with an accuracy of +/− 5 m. Plots were established with a minimum distance of 200 m between the plot centers. Inside each square plot, a circular plot with radius R = 17.84 m was installed; the center was marked with a plastic sheet and geo-referenced. The circular plots were used to capture the diversity of tree arrangements in the pasture areas. Forest and secondary forest patches were avoided at the time of circular-plot installation, due to the lack of grazing (Figure 2a). In the circular plots, trees with a diameter at breast height (DBH) ≥ 0.10 m and height ≥ 2 m were measured and recorded. Total height, commercial height, and crown diameter was measured, and the trees' geographical coordinates were also recorded.

UAS Missions and Image Processing
UAS surveys were performed using a DJI Mavic Pro quadcopter (DJI Inc, Shenzen, China), equipped with a 12-megapixel, 35-millimeter-focal-length RGB camera. A flying area of 1.56 ha (125 × 125 m), covering the square plots, was designed using the Pix4Dcapture application (Pix4Dcapture 4.3; Pix4D SA, Lausanne, Switzerland), using tablet device with an Android operating system. In the Pix4D application, flights were set to 80% forward and 72% side overlap, flying height was set to 60 m, and camera pitch was set to 90° heading forward. Mission-camera parameters were as follows: focal length, 4.695; image height, 3000; image width, 4000; minimum photo interval, 2.0; sensor width, 6.17. The surveyed area was set to cover a bigger area than the square plot (Figure 2a), in

Plot Delimitation and Tree Inventory
Farmers were asked to identify grazing areas that were representative of their farms [35]. Once the areas were identified, two 1-hectare square plots were installed. Plot corners were marked with a plastic sheet to make them visible from the UAS survey, and geographical coordinates were collected with a Garmin GPS50 GPS with an accuracy of +/− 5 m. Plots were established with a minimum distance of 200 m between the plot centers. Inside each square plot, a circular plot with radius R = 17.84 m was installed; the center was marked with a plastic sheet and geo-referenced. The circular plots were used to capture the diversity of tree arrangements in the pasture areas. Forest and secondary forest patches were avoided at the time of circular-plot installation, due to the lack of grazing (Figure 2a). In the circular plots, trees with a diameter at breast height (DBH) ≥ 0.10 m and height ≥ 2 m were measured and recorded. Total height, commercial height, and crown diameter was measured, and the trees' geographical coordinates were also recorded.

UAS Missions and Image Processing
UAS surveys were performed using a DJI Mavic Pro quadcopter (DJI Inc, Shenzen, China), equipped with a 12-megapixel, 35-millimeter-focal-length RGB camera. A flying area of 1.56 ha (125 × 125 m), covering the square plots, was designed using the Pix4Dcapture application (Pix4Dcapture 4.3; Pix4D SA, Lausanne, Switzerland), using tablet device with an Android operating system. In the Pix4D application, flights were set to 80% forward and 72% side overlap, flying height was set to 60 m, and camera pitch was set to 90 • heading forward. Mission-camera parameters were as follows: focal length, 4.695; image height, 3000; image width, 4000; minimum photo interval, 2.0; sensor width, 6.17. The surveyed area was set to cover a bigger area than the square plot (Figure 2a), in order to ensure sufficient photographs were taken for modeling of the surveyed area. The durations of the flight missions were approximately 8 min, during which between 80 and 150 pictures were obtained, following a grid flight path, which is illustrated in Figure 2a. Images were automatically geotagged with the UAS-incorporated GPS. As a preliminary step, poor-quality images (featuring blurriness, or a bad photograph ensemble) were manually selected by the operator and discarded. The remaining pictures were included in the photogrammetric analysis, which was performed using Web Open Drone Map software (OpenDroneMap, https: //github.com/OpenDroneMap/ODM/blob/master/README.md, 10 April 2022) [36]. For the point-cloud creation, three processing parameters were selected: Web Open Drone Map default settings (def), and two alternative versions that differed from the default settings, except for the resolution parameters, which were increased-resolution (par1) and highresolution (hres) settings (Table 1). Processing parameters are characterized in Appendix A. A digital-surface model (DSM), an orthomosaic (OM), and an unclassified dense point cloud (PC) were obtained from photogrammetric analysis (Figure 2b). Due to the time spent on processing under hres, the PC was not considered in the point-cloud-classification process.

Point-Cloud Classification and Canopy-Height-Model Building
A canopy-height model (CHM) is an object containing the height of the vegetation aboveground and is calculated as the difference between a DSM and a digital terrain model (DTM). Although the software allows a DTM to be obtained as a primary output of the photogrammetric analysis, site characteristics and image quality can produce PC misclassification, which can have an impact on the CHM and, subsequently, on the estimation of trees and forest characteristics. To reduce the misclassification issues, classification was controlled and performed in CloudCompare v2.11 software [37], where PCs were imported. As a first step, points that were visually identified as outliers (below the PC surface or extremely distant from the PC top) were removed [38], and, later, the PC classification was performed using the Cloth Simulation Filter (CSF) plugin [38], as explained in Figure 3a. The algorithm used in the CSF inverts the PC and simulates a rigid cloth on the inverted surface. The nearest points to the surface drawn by the cloth are classified as ground points, and the remaining points are classified as non-ground points (PCng). The CSF has four sets of parameters that can be selected, according to the project characteristics. The scene parameter allows cloth rigidness to be selected from the following options: very rigid (flat), rigid (relief), and soft (steep slope). The slope-processing parameter (yes-no) allows the consideration of terraced slopes. The underlying process for cloth resolution and classification threshold parameters can be found in Zhang et al. (2016). In this study, we evaluated three cloth resolutions (0.2, 0.5, and 1.0) and three classification threshold sizes (0.2, 0.5, and 1.0). Since the best combination of parameters was unknown for each of the plot characteristics, an exploratory combination of all the parameters was performed, except for scene, where an operator visually identified the plot characteristics and selected one or two of the options. Between 18 and 36 classified PCs were produced for every set of PCs; combinations of CSF classification parameters used for each plot are specified in Table 2.
Trees or other non-ground structures, identified in the original PC, are shown as a lack of points in the ground point cloud (Figure 3a). To fill these gaps, the raster tool was used with the average scalar field interpolation method in the CloudCompare software. A filled-ground PC (PCfg) was generated as the result of the rasterization process, which can be interpreted as a DTM for the study area. Finally, a CHM was computed as the difference between PCng and PCfg, using the 2.5D volume tool and exported as a Geotiff file.
Total 684 Trees or other non-ground structures, identified in the original PC, are shown as a lack of points in the ground point cloud (Figure 3a). To fill these gaps, the raster tool was used with the average scalar field interpolation method in the CloudCompare software. A filled-ground PC (PCfg) was generated as the result of the rasterization process, which can be interpreted as a DTM for the study area. Finally, a CHM was computed as the difference between PCng and PCfg, using the 2.5D volume tool and exported as a Geotiff file.

Individual Tree Detection and Tree Characteristics
CHMs were imported into the R environment [39], where, first, we tested whether a smoothing process could improve tree detection. Smoothing was performed using the CHMsmooting function, from the rLiDAR package [40], using a median filter and three pixel sizes (sm3 = 3 × 3; sm5 = 5 × 5; sm7 = 7 × 7), which were compared against a non-smoothed CHM (sm0) version [24]. For individual tree detection, the local-maxima algorithm, with the variable-window-filter function (vwf) from the ForestTools package [41] was used, with a minimum height of 2 m. Five window sizes (w1 = 1, w2 = 2, w3 = 3, w4 = 4, w5 = 5) were tested to explore the differences in the individual tree detection. Subsequently, the marker-controlled watershed segmentation function, with a minimum height of 1.5 m, was used to calculate crown height, diameter, and area for every identified tree. The entire process can be seen in Figure 3b. As a result, for every combination of smoothing level and window size applied to a CHM, a polygon shapefile with information about individual tree geographic coordinates, height (hm), crown area (tcm), and crown diameter (cdm) was obtained.
To compare the results from the photogrammetric analysis and field measurements, the shapefiles from the previous steps were masked into the circular plot area, containing information about tree characteristics that were measured on the field ( Figure 2a). First, we selected only the models that identified the exact number of trees, as measured on the circular plots. Subsequently, root mean square error (RMSE) was calculated to assess the difference in tree height and crown diameter between field measurements and values obtained from the UAS. For each parameter, the 5 outputs with smallest RMSEs were selected.

Tree Density
Individual tree-detection spatial accuracy was evaluated using a methodology similar to [24], which consists of the comparison between modelled and independent visual assessment of trees on the UAS-derived ortho-mosaic. Ten squares of 400 m 2 were randomly distributed on each farm, and, in these, an operator visually identified individual trees and generated a points shapefile containing tree coordinates for each validation square, using QGIS software (QGIS.ORG, Bern, Switzerland) ( Figure 4).  Tree density was merged with validation squares and later imported to the R environment, where the validation squares were used to crop the tree-points shapefile. A dataframe containing the number of observed trees and modeled trees was obtained. Density difference (Difdensity) was calculated as the average absolute (Abs) difference Tree density was merged with validation squares and later imported to the R environment, where the validation squares were used to crop the tree-points shapefile. A dataframe containing the number of observed trees and modeled trees was obtained. Density difference (Difdensity) was calculated as the average absolute (Abs) difference between observed (Obs) and modeled trees (Mod), expressed in terms of percentage (Equation (1)). The parameter combination that delivered the five lowest Difdensity values was selected. (1)

Cloud Processing, Point-Cloud Classification, and Individual-Tree-Detection Selection Parameters
The ideal parameter combination for tree height, crown diameter, and tree density was identified using frequency tables from the combinations that presented 5 lowest RMSEs for each category. For each cloud-processing level (odm), cloud-classification variable (scene, slope processing, cloth resolution, classification threshold) and individual tree-detection parameter (smoothing, window size), the relative frequency was calculated, with values ranging from 0-1. Values close to 1 showed a higher frequency of a given level in a variable, and if values were closer to 0, then there was not predominance of a given level in the studied variable. Values between 0.76-1 were considered to be in strong agreement and are shown in green, whilst a value between 0.51 and 0.75 was considered to be in medium agreement and is shown in yellow. Values equal to or lower than 0.50 were considered to be in weak agreement and are shown in red.

Tree Spatial Arrangements
Tree spatial arrangements were assessed using the best combination of parameters for tree height, crown diameter and tree density in the previous section. Prior to the spatial analysis, shapefiles containing tree location were cropped into a smaller area, to avoid border effects on individual tree identification. A distance-based approach to cluster identification was used, in which clusters are formed when the distance between at least two trees is less than a given distance d [42]. A distance d = 6 m was used as the criterion for cluster definition [43]. Since spatial coordinates for every modeled tree were known, the buffer function from the raster package in R was used [44] to create a buffer area of radius = d around each tree, as shown in Figure 4. Next, the overlapping buffer areas were merged into a single polygon, allowing the spatial arrangements to be identified. Spatial arrangements were classified as scattered trees (individual trees with distance >6 m to the nearest neighbor) and clumps (two or more trees with a minimum distance ≤6 m between trees) [43]. Clumps were subclassified into small (2-9 trees), medium (10-19 trees), and large clumps (>20 trees). Results are presented in terms of percentage of trees under these category arrangements.

Aboveground Biomass
Aboveground biomass from measured values at the circular plot level was calculated using Equation (2) [45]. This equation requires the wood-density parameter, which was retrieved at genus level from the Global Wood Density Database [46] using the getWoodDensity function from the BIOMASS package in R [47]. When genus information was not available, wood density was estimated as the average density at site level (Chave_wdsite) [48] and family level (Chave_wdfamily). Biomass was estimated from measured values, using Equation (3), proposed by [49], which was developed for remote-sensing imagery; the equation relies on tree-height and crown-diameter parameters. The use of Equation (3) with measured values was completed to estimate the differences in measured values separately from the values obtained from the UAS imagery. Using the tree-height and crown-diameter information derived from the CHM, Equation (3)  pared against values obtained in the field. Finally, the share of above-ground biomass was calculated for each spatial-arrangement category.

Site Characteristics and Tree Inventory Results
The circular plot installation (0.1 ha) and tree measurements on the pasture areas each took 2-5 h, depending on the number of trees, using teams of four people. For the UAS survey, a team of two people and 45-60 min were needed for every 1.56-hectare plot. Since the surveys were performed on areas that farmers use for cattle grazing, open areas covered with grass species were dominant. The terrain slopes varied between the plots, but, in general, irregular terrains with slopes ≥ 20% were found in the Amazon region (plots A1, A2, A3, and A5), whilst in the coastal region, the slopes ranged between 4 and 16% in four plots; only in plot C4 was this value above 20% (Appendix B). In total, 58 trees were surveyed, with a minimum of 1 and a maximum of 17 trees per plot. The tree heights ranged from 2.4 to 21 m, and the crown diameters ranged from 2.5 to 15.6 m.

Cloud-Point Creation
The differences in the point-cloud building settings had an impact on the processing time, point-cloud size, and percentage of discarded points, with par1 settings using more processing time, displaying smaller point-cloud sizes and removing fewer points than the def characteristics (SI2). The processing time for all the plots was 12:02:22 h under default characteristics, 16:49:56 h for par1 settings, and 34:16:45 h for high-resolution settings. For the point-cloud size, the par1 point cloud had, on average, 4,281,057 fewer points than the default characteristics, whilst the high resolution had 32,430,874 more points. The percentage of points identified as outliers and removed varied from 0.32, 0.059, and 0.20% for the default, par1, and high-resolution settings, respectively. Due to the processing time needed for the high resolution, the next steps considered only the def and par1 outputs. From the def and par1 settings, 13,680 outputs with individual tree detection were generated. Of these, only 5078 detected an exact number of trees inside the circular-plot validation, making them suitable for assessing differences among measured values. Due to the differences in site characteristics, the point-cloud classification and individual tree-detection parameter combination were expected to vary between plots, but also to show similarities in plots with similar characteristics. The similarities between the plots are discussed in the Discussion section.

Canopy-Height-Model Validation
For the calculation of the biomass from the field measurements, in which the correct measurement of the tree parameters resulted in more robust estimations, the tree parameters obtained from the UAS surveys needed to be validated against field measurements to prove their usefulness on biomass estimations. The tree heights and crown diameters were compared against the field measurements to identify possible sources of error in the final biomass estimates.

Tree-Height Differences
The RMSEs of the best models ranged from 0.49 to 3.53 m and 0.01 to 1.01 m in the Amazon and the Coastal regions, respectively (Figure 5a). The point-cloud-building, classification, and individual tree-detection parameters did not significantly influence the estimates of tree height when analyzed separately. The relative frequencies for each parameter level for the tree height sand crown diameters are presented in Table 2. The point-cloud building specifications (ODM) resulted in high agreement (>0.75) in all the plots, with 80% of the plots showing lower RMSE under the def settings in the Web Open Drone Map software, whilst in the plots A2 and C4, the best results were found using the par1 specifications. With the point-cloud classification parameters, only the cloth resolution (CLR) and classification threshold (CTR) presented a high agreement in more than 50% of the plots, with a CLR size of 1 for 70% of the plots and CTR of 0.2 for 50% of the plots. Meanwhile, for the scene and slope processing parameters, medium agreement within plots was found in 50% and 60% of the plots, respectively. The flat scenes performed better for the plots in the Coastal region, whereas the steep slope and relief scenes were predominant in the plots located in the Amazon region. In 60% of the farms, the tree heights were estimated better without slope processing, whilst a cloth resolution size of 1 and a classification threshold of 0.2 were the most common settings for at least 50% of the plots. The parameters for individual tree detection reached strong agreement in 40% of the farms for window-filter size, and medium agreement in 50% of the farms for smoothing-window size. Only window-filter sizes 1 and 5 presented strong agreement, and the smoothing-window process presented strong agreement at sizes 0, 3, and 5.

Crown Diameter Differences
The RMSEs for the best models ranged from 0.04 to 1.77 m and from 0.53 to 4.47 m for the Amazon and Coastal regions, respectively ( Figure 5B). Similarly to the height estimation, the ODM parameter presented strong agreement in 90% of the plots, among which default characteristics for point-cloud building were present in 80% (Table 3). For the point-cloud-classification parameters, scene resulted in high agreement in 40% of the farms, whilst the other 60% reached only medium agreement. The steep slope was the most common scene for the plots in the Amazon region and flat was the most common

Crown Diameter Differences
The RMSEs for the best models ranged from 0.04 to 1.77 m and from 0.53 to 4.47 m for the Amazon and Coastal regions, respectively ( Figure 5B). Similarly to the height estimation, the ODM parameter presented strong agreement in 90% of the plots, among which default characteristics for point-cloud building were present in 80% (Table 3). For the point-cloud-classification parameters, scene resulted in high agreement in 40% of the farms, whilst the other 60% reached only medium agreement. The steep slope was the most common scene for the plots in the Amazon region and flat was the most common for the Coastal region. Avoiding slope processing presented the best results in 70% of the farms, with strong agreement in 50% of them. Cloth-resolution sizes of 0.2 and 1 were present in 40 and 50% of the plots, respectively, and showed strong agreement in only three plots. The classification threshold reached strong agreement in 50% of the plots, with a size of 1 present in 40% of the plots, followed by 0.5 in 30% of plots, and 0.2 in only two plots. The window-filter size resulted in strong agreement in only one plot and medium agreement in one, but with a size of 5. The smoothing process resulted in agreement in 50% of the plots, with size 5 present in 30% of the plots, and size 3 in 20%.

Tree-Density Assessment
The tree-density differences ranged from 25 to 48% and from 21.85 to 64.28% in the Amazon and Coastal regions, respectively ( Figure 5C). The point-cloud-building specifications showed strong agreement in 70% of the plots for the default characteristics, and medium agreement was found for three plots, with one def and two par1 characteristics (Table 4). For the cloud-classification parameters, scene resulted in strong agreement in 40% of the plots, with two flat and two s.slop; meanwhile, 50% of the farms resulted in medium agreement, with three plots with a relief scene, and one with a flat and steep slope. For the slope processing, medium agreement was found in 80% of the plots, with 50% of the plots resulting in better tree-density estimates with slope processing and the other 50% without it. Cloth resolution resulted in strong agreement in 50% of the plots, whilst the other 30% showed medium agreement, and 20% showed low agreement. The classification threshold resulted in medium agreement for 50% of the plots, and strong agreement was found in only one plot. The most common cloth-resolution size was 0.5 (3 plots), followed by 0.2, and 1, with one plot each. Strong agreement was absent in the window-size-filter parameter, while medium agreement was found only in four plots, with size 2 present in three of them. No-smoothing showed medium and strong agreement in two plots. Using the best models for tree height, crown diameter, and tree density, the average relative frequency was calculated for an ideal combination of parameters that could help to improve the estimations. The relative-frequency values were averaged for every parameter level. The apparent optimum model combination showed a strong agreement for default point-cloud building on 70% of the plots, whilst in the other 30%, there was medium agreement, with one plot with par1 settings. The point-cloud classification parameters gave medium agreement for scene (70% of plots), slope processing (90% of plots), and cloth resolution (60% of the plots). For the tree-detection parameters, the window-size filter resulted in low agreement for 90% of the plots and medium agreement for one plot. The smoothing-processing window resulted in medium agreement for 50% of the plots and low agreement for the other 50%.

Spatial Arrangements
The methodology proved to be useful for identifying the spatial arrangements of trees in livestock landscapes. Although the number of individual trees detected was influenced by a combination of parameters ( Figure 6B), the average tree distribution for was used to determine the spatial arrangement. On the plots located in the Amazon region, with the exception of A5, the majority of the trees were grouped in small clumps, whilst in the plots located on the Coastal region, the majority of the trees were grouped in large clumps, as seen in Figure 6A. On plot A5, most of the trees were grouped in large clumps. The number of trees identified in the square plots varied depending on the parameter-combination selection criteria; in plots A1, A3, C2, C3, and C5, the parameter-combination selection that focused on tree-density validation resulted in a higher number of trees, whilst in plots A2 and A5, the parameter combination that showed the lowest RMSE for tree density identified a higher number of trees in comparison with the other parameter-combination criteria. For plots A4, C1, and C4, only one parameter-combination criterion identified the exact number of trees inside the circular plots.

Biomass Estimation
With information obtained from the field measurements, the difference in biomass estimates as the result of the equation selection was compared. The biomass calculated using Equation (3)

Biomass Estimation
With information obtained from the field measurements, the difference in biomass estimates as the result of the equation selection was compared. The biomass calculated using Equation (3) showed higher mean values compared with Equation (2) in seven plots (Figure 7), although the standard deviations showed similar values. The average biomass estimated with Equation (3) presented higher values in plots A1, A2, A4, A5, C2, and C5, whilst in A3, C3, and C4, the Equation (2) estimates were higher.
When comparing the biomass estimates from the field measurements against those obtained from the UAS survey, the results and best parameter combinations varied according to the plot. This was expected, since the performance and reliability of UAS surveys is greatly influenced by the site conditions; consequently, the point-cloud classification and individual tree-detection parameters varied between the sites.
The sum of the biomass per plot allows us to compare the differences between parameter combinations and biomass obtained from measured values. When comparing against biomass estimated using Equation (2), the best estimates of biomass for plots A2 (average combination), A5 (lowest RMSE for tree density), C1 (crown diameter), and C5 (tree density) were 127.96, 3324.99, 131.28, and 929.09 kg, respectively, with these values differing from the measured values by less than 25% (Figure 8). In these plots, the biomass estimated using Equation (2) was lower than that estimated with Equation (3). In plot A1, with the parameter combination for tree height, a biomass estimate of 2446.34 kg was obtained, which was 36% lower than that estimated using Equation (2). In plots A3, A4, C2, C3, and C4, the biomass values estimated from the UAS imagery were 88.46, 490.40, 2.26, 219.76, and 0.028289 kg, respectively. These values differed by more than 75% from the values estimated from the measured tree characteristics on the field using Equation (2). The parameter combination used in these plots was crown diameter (A3, C2), the lowest RMSE for tree height (A4, C3), and tree density (C4).
When comparing the biomass estimates derived from the UAS metrics against those obtained using the measured values and Equation (3), for plots A5, C1, and C5, the parameter combination for tree height, the lowest RMSE for crown diameter, and the lowest RMSE for tree density resulted in differences of less than 25% compared to the measured values. For the rest of the plots, the differences were greater than 50%.

Biomass Estimation
With information obtained from the field measurements, the difference in biomass estimates as the result of the equation selection was compared. The biomass calculated using Equation (3) showed higher mean values compared with Equation (2) in seven plots (Figure 7), although the standard deviations showed similar values. The average biomass estimated with Equation (3) presented higher values in plots A1, A2, A4, A5, C2, and C5, whilst in A3, C3, and C4, the Equation (2) estimates were higher.  When comparing the biomass estimates from the field measurements against those obtained from the UAS survey, the results and best parameter combinations varied according to the plot. This was expected, since the performance and reliability of UAS surveys is greatly influenced by the site conditions; consequently, the point-cloud classification and individual tree-detection parameters varied between the sites.
The sum of the biomass per plot allows us to compare the differences between parameter combinations and biomass obtained from measured values. When comparing against biomass estimated using Equation (2), the best estimates of biomass for plots A2 (average combination), A5 (lowest RMSE for tree density), C1 (crown diameter), and C5 (tree density) were 127.96, 3324.99, 131.28, and 929.09 kg, respectively, with these values differing from the measured values by less than 25% (Figure 8). In these plots, the biomass estimated using Equation (2) was lower than that estimated with Equation (3). In plot A1, with the parameter combination for tree height, a biomass estimate of 2446.34 kg was obtained, which was 36% lower than that estimated using Equation (2). In plots A3, A4, C2, C3, and C4, the biomass values estimated from the UAS imagery were 88.46, 490.40, 2.26, 219.76, and 0.028289 kg, respectively. These values differed by more than 75% from the values estimated from the measured tree characteristics on the field using Equation (2). The parameter combination used in these plots was crown diameter (A3, C2), the lowest RMSE for tree height (A4, C3), and tree density (C4).
When comparing the biomass estimates derived from the UAS metrics against those obtained using the measured values and Equation (3), for plots A5, C1, and C5, the parameter combination for tree height, the lowest RMSE for crown diameter, and the lowest RMSE for tree density resulted in differences of less than 25% compared to the measured values. For the rest of the plots, the differences were greater than 50%. Figure 8. Comparison of total aboveground biomass for plots obtained from different equations and parameter-combination selection criteria. Field_Chave = biomass estimated from field measurements using Equation (2); Field_Jucker = biomass estimated from field measurements using Equation (3); Mod_combination = PC classification and individual tree-detection parameters described in Table 5; Mod_density = PC classification and individual tree-detection parameter combination for tree density; Mod_diameter = PC classification and individual tree-detection parameter combination for crown diameter; Mod_height = PC classification and individual treedetection parameter combination for tree height; Mod_RMSE_density = PC classification and individual tree-detection parameter combination that presented the lowest density difference for tree density; Mod_RMSE_diameter = PC classification and individual tree-detection parameter Figure 8. Comparison of total aboveground biomass for plots obtained from different equations and parameter-combination selection criteria. Field_Chave = biomass estimated from field measurements using Equation (2); Field_Jucker = biomass estimated from field measurements using Equation (3); Mod_combination = PC classification and individual tree-detection parameters described in Table 5; Mod_density = PC classification and individual tree-detection parameter combination for tree density; Mod_diameter = PC classification and individual tree-detection parameter combination for crown diameter; Mod_height = PC classification and individual tree-detection parameter combination for tree height; Mod_RMSE_density = PC classification and individual tree-detection parameter combination that presented the lowest density difference for tree density; Mod_RMSE_diameter = PC classification and individual tree-detection parameter combination that presented the lowest RMSE for crown diameter; Mod_RMSE_height = PC classification and individual tree-detection parameter combination that presented the lowest RMSE for tree height.

Discussion
Traditionally, the study of trees on pasture areas has been performed using measurement plots, the results from which are then extrapolated to the entire pasture area on the studied farms or region. At larger scales, other approaches have used tree coverage information (in terms of percentage), which are superimposed onto the pasture categories on land-cover maps [50]. While both approaches have improved the understanding of the contribution of trees on pasture areas to carbon budgets, neither adequately captures the heterogeneity within and between farms or landscapes. The use of UASs can be used to complement traditional field measurements and, possibly, to provide an intermediate step between ground-and satellite-based remote sensing.
In the past, only a few studies explored the use of UASs for the study of tree characteristics in livestock landscapes [29,30]. In this paper, a new methodology is proposed for using UAS surveys for the assessment of woody vegetation in silvopastoral systems. The methodology relies mainly on open-source (OS) software, making it more accessible in developing countries, where research funding can be scarce [18,51]. Although commercial software dominates the field, an increasing number of studies use open source software [52][53][54][55]. In terms of the processing steps required, in general, increasing the resolution of the photogrammetric analysis resulted in an increase in total processing time of 04:47:34 h from the def to the par1 and 22:14:23 h from the def to the hres settings. Although the proportion of removed points was lower under par1 and high resolution, the best results were found under default characteristics (Tables 4 and 5), making it unnecessary to invest more time, since no improvement could be obtained in the desired outputs. When assessing the performance of different photogrammetry OS software, the authors of [55] found that WODM default settings performed better in comparison with other OS software, but were less accurate than commercial settings. The plots included in this study differed in tree density and terrain characteristics (Appendix B), which ultimately had an effect on the ability to correctly identify tree characteristics and biomass estimates (Figure 8). For example, in plot C4, field measurements recorded only one tree 3.1 m in height, probably causing a low point density in the point cloud and, subsequently, an inability to identify it from the UAV survey data.

Tree Parameters and Biomass Estimation
Most of the work using UAS for the study of forest characteristics relies either on high-resolution DTM from LiDAR surveys [26,56] or high-accuracy GPS [28,57] to calculate CHM from DTM. This approach relies entirely on the performance of point-cloudclassification algorithms to correctly differentiate ground from non-ground structures. The authors of [23] tested the use of UAS-independent metrics to predict forest-growing stocks in flat forest areas in Norway and Italy. Their findings suggest that it is possible to rely only on UAS-DTM-independent metrics to estimate ground-stand volume, with RMSE values similar to those obtained when using a high-resolution DTM. Site conditions, such as terrain slope and tree cover, have an influence on point-cloud-classification performance, although the authors of [30] suggest that tree coverage is more important for accurate DTM creation. In this study, three of the four plots that had a RMSE between 1.14 and 2.5 m (≥10% relative RMSE) were located in areas with a slope of between 20 and 29 • . With the exception of plot A1, the results of this research include lower RMSE values than those reported by the authors of [28] (1.45 m), who used high-accuracy GPS in a study on restoration areas with terrain slopes between 15 and 30 • in Costa Rica, and the values reported by [26,56] on forest areas with closed canopies in Japan and Sumatra, respectively. The error associated with tree-height estimation is also influenced by the height of the studied trees; the authors of [29] reported an RMSE of 3.8 m in trees with heights equal to or lower than 1 m, with a consistent decrease until 1.5 m for trees 2 m in height, and a 0.9-meter RMSE for trees with a height of 8 m. Since this study only considered trees with heights equal to or above 2 m, with the exception of plot A1, the results showed a similar performance. When analyzing the performance of UAS-derived CHM models in estimating crown diameter, the worst performance was found for plot C2, with an average RMSE of 4.43 m, followed by plots A1, C3, and C5, with 1.74, 1.66, and 1.05 m, respectively. In the rest of the plots, the RMSE values varied between 0.1 and 1 m. The authors of [58] reported RMSE values of 0.73 m for crown width derived from photogrammetry analysis in a eucalyptus plantation in Australia.
Most allometric equations use tree height and DBH as the main parameters to estimate tree biomass [45]. Nevertheless, it is difficult to obtain DBH from traditional UAS surveys. Whilst previous work has assessed the suitability of using UAS for a full reconstruction of trees [53], this could require different flying parameters, requiring more time and making the UAS survey less efficient. The authors of [48] proposed a specific equation for using remote-sensing data in tropical forests, which was proven to result in more reliable biomass estimates when compared with field measurements [48]. Although several papers have reported a relationship between crown width and DBH [48,55], this relationship might be different in trees in agricultural landscapes, because growing conditions differ from those in natural forests, and/or because management practices, such as pruning and foraging, modify crown development. Hence, the crown-DBH relationship was altered, which may have exerted an impact on the final biomass estimates. For example, the authors of [30] reported a low R2 for equations using crown diameter as a predictor of DBH on Tectona grandis trees in SPS in Costa Rica. Although biomass estimates present differences in comparison to biomass estimates derived from field measurements, by assessing the error magnitude, these estimates can be useful for reporting carbon stocks in SPS at larger spatial extents.

Tree-Density Assessment
The preliminary tests of the method for individual tree detection showed that smaller window sizes produced more trees per plot, in some cases splitting large trees into several small trees. Since the nature of this work was exploratory, the range of the detected trees varied across the parameter combinations. Nevertheless, the relative frequency values presented medium agreement for only four plots, with window sizes of 5 for plot A1 and window sizes of 2 for plots A5, C3, and C5. The results suggest that while window size can be an important predictor for correct tree-density assessment, the point-cloud-building and classification parameters played a larger role in the tree-density estimation. The authors of [29] reported that the proportion of trees found was dependent on the tree height, with trees with heights below 4 m displaying identification rates below 75%, which increased to 90% for trees taller than 4 m. The authors of [24], using a watershed algorithm, identified 93.8% of the trees in areas with a density of 150 to 325 trees/ha. Due to the quantity of output from the parameter combination (13,680), these results did not consider omission and commission errors separately for every plot; instead, these errors were integrated into a comprehensible and easy-to-calculate index, which represents the percentage difference with visually identified trees, independent of the number of trees in each subplot.

Spatial Arrangement Assessment
Spatial arrangements have been identified from field measurements by other authors [59], but the problems with these are the same as those reported for tree measurements, i.e., the elevated costs associated with identifying all of the arrangements. One of the negative effects of trees in pastures, as perceived by farmers [60], is the shades that trees cause, and their detrimental effect on grass biomass production [61]. UAS surveys in livestock landscapes could allow areas with excess shade to be identified and, further, allow the distribution of trees to be improved, mitigating the negative effects and increasing the positive effects. The methodology proved useful for identifying tree arrangements on pasture areas, although, in this study, we only classified the arrangements into four categories. Analyzing the clump characteristics (length and width) can improve the classification and contribute to the inference of the functionality of spatial arrangements, depending on their locations on farms.

Conclusions
The methodology described in this study provides a tractable and economical approach to studying the presence of trees, as well as their arrangements, spatial distributions, and individual characteristics (such as height and crown diameter) in livestock landscapes, which can be used to estimate ecosystem services, such as carbon storage. Point-cloud classification remains a challenge for improved DTM generation, especially in areas with highly irregular terrain. In this study, tree species were identified with field observations at circular plot level, but this information was not considered at UAS survey level. Thus, the biomass estimates used a generic equation for the entire plot. Accounting for all the species at landscape level includes the same problems as measuring trees at that spatial extent. To overcome this difficulty, one option is the use of orthomosaics for species identification [33,62]. Although uncertainty over UAS-derived biomass estimates is still high, the possibility of allocating those estimates spatially, across spatial arrangements and landscapes, can be seen as an improvement in understanding the factors influencing carbon stocks in livestock or agricultural landscapes. The generation of orthomosaics from UAS surveys can also be useful for communicating the importance of tree presence in livestock landscapes to farmers and other stakeholders. The visual assessment of the lack/abundance of trees in the landscape might open up opportunities for reflection and discussion about the role of trees on farms.
A methodology that is suitable for use in developing countries was proposed. It is acknowledged that the use of RTK/PPKGPS, NDVI, or LiDAR sensors may lead to an improvement in the results, but high-quality science should also be affordable as shown using the methods described in this paper, and made more accessible to developing regions, where a lack of resources limits the availability of more sophisticated and expensive sensors and commercial tools. Acknowledgments: The authors are grateful for the help from the field teams of the Climate Smart Livestock Project (CSLP), Ecuador, especially to Jhonny Delva and Juan Merino from CSLP for their support on the coordination and field data collection. Since this research is part of a PhD thesis, we are thankful to the University of Aberdeen, UK, for the Elphinstone scholarship and Universidad de Cuenca, Ecuador, for financial support.

Conflicts of Interest:
The authors declare no conflict of interest.