Citrus orchards under formation evaluated by UAV-Based RGB Imagery

Few studies have investigated the biometric attributes of citrus orchards under formation that use RGB sensors on board unmanned aerial vehicles (UAV) and the challenges are great. This study aimed to develop and validate a method of using aerial UAV images by automated routines to evaluate the biometric attributes of a crop of ‘Tahiti’ acid lime under formation. We used a multirotor UAV, programmed to capture images at three different map scales, with a frontal and side overlap of 80 %. Geoprocessing was carried out both with and without ground control points on each scale. An automated routine was developed in an opensource environment, consisting of three processing phases: i) Estimation of the plant biometric attributes, ii) Statistical analysis, and iii) Statistical Report Map (SRM). The use of the developed routine allowed to delimit and estimate the crown projection area with an accuracy of more than 95 % as well as identify and quantify the plants with an accuracy of over 97 %. The use of ground control points during the processing stage does not increase accuracy in estimating the biometric attributes under evaluation. On the other hand, map scale is strongly correlated with the quality of the estimates, especially plant height. The results allowed to define a method for the acquisition and analysis of aerophotogrammetric data using a UAV, which can be used to measure the plant biometric attributes under analysis and the method can be easily adapted to


Introduction
Standardizing processes and developing solutions (algorithms) for extracting biometric attributes from plants are essential when a large amount of data from UAV sensors is used (Aasen et al., 2018;Tu et al., 2020), especially in search for low-cost solutions (Ampatzidis et al., 2020). These authors report that different combinations of flight parameters can affect product quality and application. For Osco et al. (2020), no computer vision technique is universally applicable and different approaches should be tested to solve specific problems. Fawcett et al. (2019), Torres-Sánchez et al. (2018), and Zarco-Tejada et al. (2014) worked on optimizing UAV flight parameters to measure plant structure. Tu et al. (2020) conducted a comprehensive study, including a series of flight parameters that could be adopted to measure height of avocado trees, but only when a multispectral sensor was used. Seifert et al. (2019) analyzed several parameters with an RGB sensor for application in the forestry sector.
Despite advances in methods to extract biometric attributes (survival inventory, crown projection area, plant height and volume, planting failures, and average spacing between plants and rows) of citrus plants, no studies have investigated the effects of variations in aerophotogrammetric parameters on the quality of the automatic extraction of such attributes using a low-cost sensor (RGB). This is even more challenging in orchards under formation, where greater spatial variability of such attributes might occur and, if left untreated, orchard yield is certainly compromised. Therefore, estimation of biometric attributes using an on-board UAV sensor proposes various questions: a) Is it necessary to use ground control points in processing the aerophotogrammetric data? b) What is the effect of map scale (or the size of the ground sample distance-GSD) on the estimation of these attributes? c) Is it possible to develop low-cost and intuitive automated computer routines to analyze ortho-photomosaics? Therefore, this study aimed to develop and validate a method to use aerial images (obtained with an RGB sensor of on-board UAV) by automated routines to evaluate biometric attributes of a 'Tahiti' acid lime crop under formation.

Collecting biometric attributes
The experimental site consisted of 24 rows of 39 plants row - level; and crown diameter, in the same direction (Dl), perpendicular to the row (Dr), both obtained using a millimeter tape measure (positioned at half of the height of the plant under evaluation). The plant volume was then estimated according to the formula presented by Rodrigues et al. (2019), Eq. (1): where: V is the volume (m 3 ) of the crown, and Dm is the diameter (m) of the crown, obtained from the arithmetic mean of (Dl) and (Dr).

Collecting and processing the geodetic and aerophotogrammetric data
Eleven support points were set up: five ground control points (GCP) for the absolute orientation of the aerophotogrammetric block, and six checkpoints (CP) to control quality (accuracy), both systematically distributed within the area of interest (Figure 1), according to recommendations of Aasen et al. (2018). The respective points were georeferenced by a Geodesic Receiver, with a mean positional accuracy (sigma 95 %) of less than 3 cm (planimetry) and 4 cm (altimetry). Flight planning considered a quadcopter-type UAV (DJI Phantom 4 Adv), equipped with a 20-megapixel RGB sensor, 13.2 mm × 8.8 mm in size, with a focal length of 9.1561 mm. Three flight altitudes [60 m (F60), 90 m (F90) and 120 m (F120)] were defined to get a reference scale, that is, the ground sample distance (GSD).
The flight included a horizontal buffer of 60, 120, and 180 m added to the area surrounding the experimental site for F60, F90, and F120 flights, respectively. This procedure ensured a minimum overlap of nine images point -1 and minimized the rolling shutter effect and variations in the UAV pitch angle (Tu et al., 2020). This procedure can also reduce false matching and increase point cloud accuracy (AGISOT LLC, 2020), improving quality of the ortho-photomosaic (Shi et al., 2016).
The flights were planned along the same orientation of the rows, as recommended by Tu et al. (2020). Finally, the compass was calibrated and the flights carried out on 19 May 2020, between 11h00 and 12h00 to minimize the shading effects on the images, taking into account various flight parameters and RGB sensor settings (Table  1), as recommended by Fawcett et al. (2019), Seifert et al. (2019), Shi et al. (2016), and Wang et al. (2019). During flight F120, the lighting conditions remained unstable (partly cloudy).
The images collected were imported into the Agisoft Metashape® software for geoprocessing (Structure from Motion-SfM) and then the Digital Surface Model (DSM), Digital Terrain Model (DTM) and Classic Orthophotomosaic (COM-RGB) were generated. Two processes were carried out for each flight, one with and one without ground control points (GCP). For the aerophotogrammetric blocks geoprocessed using GCP, the following quality was obtained (RMSE-X,Y, and Z) for the six checkpoints at

Implementing the automated routine
An automated routine (FindCITRUS-V1) was developed to estimate the biometric attributes of the plants, consisting of three processing phases: i) Estimation of Plant Biometric Attributes; ii) Statistical Analysis, and iii) Statistical Report Map (SRM).

Estimating plant biometric attributes -Phase I
Phase I was implemented using the Graphical Modeler of the QGIS 3.10.7-A software and started with the (mandatory) input of COM-RGB, DTM and DSM ( Figure 2). The vegetation index (VI) of the Tract is then obtained using the VARIgreen metric proposed by Gitelson et al. (2002). The Digital Height Model (DHM) is determined based on differences between the DSM and DTM, representing the relative height (rh) of each object (natural and/or artificial elements) above the surface. Pixels with an rh value of less than 45 cm were considered zero, due to the presence of objects, such as boxes and tools, as well as plant cover (between crop rows), The next process is to combine (multiplying) a geometric product (DHM) with a spectral product (VI of the Tract), where the result is used as the input to the r.recode algorithm. The use of r.recode allows to establish a threshold (analyzing the histogram) with a zero value, in which lower values were classified (and vectored) as 'non-crown area' and larger values as 'possibly crown area'; the areas classified as 'possibly crown' were extracted and smoothed. The area (m 2 ) was also calculated for the geometry of each feature. Final extraction of the crown area was only possible using a threshold (which can be changed by the user) to reduce undesirable polygons (noise). In our study, a value of 0.10 m 2 was adopted. Once the crown geometry was defined for each plant, the zonal statistics (arithmetic mean) of VI raster for the Tract (1°) was extracted to generate the Plant VI (3°b).
The next step was to obtain the crown centroid of each plant, followed by a filter that allowed to reduce False Positives (FP) when necessary thereby obtaining the geographical location of the plants. The FP Filter was implemented using the Delaunay algorithm (available from QGIS) to calculate the area, length, and azimuth of each edge of the triangular mesh. Next, a cut-off threshold was established to filter the unwanted edges with the consequent removal of any unwanted points, that is, FP.
Once each plant is located [x,y], Voronoi Polygons are created to give the spatial distribution of the occupied area (m 2 per plant). Another Delaunay Triangulation was then implemented following the logic above. However, edges outside the normal spacing between plants are filtered out by the cut-off threshold, that is, it identifies any lines that represent planting failures, allowing estimating the number of seedlings to be replanted.
Furthermore, location of the geometry that represents the crown area of each plant allowed the application of zonal statistics to the DHM raster, obtaining the maximum value for relative height, considered as plant height (h D HM  , 5°a) derived from the DHM. As an alternative, the routine allows only the DSM to be used to estimate plant height (h DSM  , 5°b). In this case, a mask is defined from the plant location

Statistical analysis -Phase II
This phase consisted of the Descriptive Statistical Analysis of all biometric attributes estimated previously (Phase I). The routine was implemented (by a script) into the Python 3.6 terminal of QGIS, using the following libraries: gdal, math, matplolib, numpy, pandas, seaborn, and scipy. The layers (vector and raster) were read using the QgsVectorLayer and QgsRasterLayer classes. Values were then extracted from the attribute tables (vector layer) using the getFeatures instance. Finally, the DataFrames were built and the descriptive statistics were obtained.

Statistical Report Map (SRM) -Phase III
This phase consisted of the development of a simple and intuitive report, known as the SRM, containing the following information: Registry (Property, Owner, and Technician in charge), Cartographic, and Geospatial Biometric Attribute Maps (associated to the Descriptive Statistics of Phase II). For that purpose, a further script was developed again in the Python 3.6 terminal of QGIS, using the following libraries: Reportlab, Matplotlib, and Geopandas.

Data analysis
The statistical analysis was carried out using the estimated biometric attributes (predicted value) and those measured in the field (ground truth) to obtain a workflow to acquire aerophotogrammetric data and validate the FindCITRUS-V1 routine. The Python 3.6 terminal of QGIS was again used to develop the customized routines, including the following libraries: matplolib, numpy, pandas, sklearn, scipy.stats, seaborn, and statsmodels.

Delimiting and estimating the crown projection area, plant height, and volume
The methodology described by Jing et al. (2012) and Dong et al. (2020) was used to identify and delimit the crown areas. First, plant crowns were delimited (photointerpretation and manual vectoring in QGIS) thereby obtaining a reference geometry (RG). Next, polygons that represent plant crowns, considered here as Extracted Geometry (EG), were extracted using the computer routine. Both were analyzed into six categories (E ma : number of matched tree crowns, E nm : number of near-match tree crowns, E me : number of merged tree crowns, E sp : number of split tree crowns, E mi : number of missed tree crowns and E wr : number of wrong tree crowns). Accuracy (Producer Accuracy -PA) was estimated from Eq. (2), where E ma represents the total number of matched tree crowns in EG, E e represents the EG total number, and E wr represents the wrong EG number.
To validate the size of the Crown Projection Area, only the E ma and E nm classes were considered as estimated values (predicted value), compared with values measured in the field (ground truth), that is, a circular area. The heights measured in the field (hp) were compared with the estimated values, that is, the heights derived from the DHM (h D HM  ) and those derived from the DSM (h DSM  ). The values of 'true volumes' (V) were calculated from Eq. (1) and then compared to the volumes derived from the DHM (V D HM  ) and DSM (V DSM  ). The data were submitted to the Kolmogorov-Smirnov test of normality. The Root Mean Square Error (RMSE) and Mean Absolute Error (MAE) were then determined.

Identifying and counting individuals (plants and planting failures)
The following performance measurements and metrics were applied, according to Oliveira et al. (2020): sensitivity (Sb), algorithm capacity to detect acid lime plants; specificity (Sp), algorithm effectiveness to identify 'non-acid lime plants'; general accuracy (Ac), overall performance measurement of the proposed method; producer accuracy (Pacc), and percentage of identified acid lime plants; used to determine the performance of the method proposed.

Results and Discussion
Delimiting and estimating the crown projection area The ground control points (GCP) did not increase accuracy (PA) in the geometric delimitation of the crowns of the plants evaluated (Table 2). No significant differences were found and this can be explained by the SfM processing technique applied to the aerial images to reconstruct the 3D surfaces, later generating the digital elevation models and ortho-photomosaic. Snavely et al. (2008) highlight that the SfM process initially estimates the relative location of objects, being the absolute location (georeferenced) just one option in this process of referencing (Sanz-Ablanedo et al., 2018). Snavely et al. (2008) also state that most problems for a 3D reconstruction can only be solved using the relative position of the objects. Szeliski (2010) presents details of the mathematical modelling of 3D reconstruction. Currently, most software programs available for processing images acquired by UAV are based on the SfM technique (Ferrer-González et al., 2020). This approach, unlike traditional digital photogrammetry, solves the collinearity equations without the need for any control points, providing a sparse point cloud in an arbitrary coordinate system and a complete camera calibration (Agisoft, 2020).
Regarding the map scale (or GSD), lower PA values were found at greater GSDs, especially from F120. In the process of delimiting the crown projection, the geometric product and the spectral product are combined, the latter is hampered by the variation in lighting conditions during flight F120 (Figures 3A and  3B). This helped to reduce delimitation accuracy, since the shading effect can cause noise in the image (Yeom et al., 2019), reduce the contrast between objects (Tu et al., 2020), and consequently, reduce point cloud accuracy (Dandois et al., 2015).
Similarly, the use of GCP did not result in significant differences in the estimating process of the Crown Projection Area. Correlation magnitude remained the same (with or without GCP). However, when abstracting from flight altitude only, a reduction can be seen in the values for RMSE, MAE and r as the map scale is reduced, that is, GSD increases (Table 2).
In general, the average values for RMSE and MAE can be understood from the methods used to estimate the area, for instance, circular geometry (ground truth) and irregular polygon (predicted value) ( Figure 3B). When the map scale is increased, image detail increases (Seifert et al., 2019) and, consequently, geometric description of the actual shape of the crown projection is improved, resulting in greater differences (ground truth -predicted value), that is, higher RMSE and MAE values.
The field measurements (ground truth) to calculate the crown projection area consider only two diameter samples (Dl and Dr), whereas using the proposed methodology, irregular polygons are more consistent with reality and theoretically more accurate than the values obtained manually, simplifying crown geometries as circular in shape.
In general, irrespective of the GCP use, the correlation coefficients were greater than 0.83. Further details on the data distribution and dispersion from  extracting the crown projection area for flights with no ground control points or field measurements can be seen in Figure 4.

Identifying and counting individuals (plants and planting failures)
Identifying the crown center of each plant by the centroid algorithm is directly related to accuracy of the previous phase (delimiting the crown projection area), mainly to the E ma and E nm values, and such accuracy may be degraded by the amount of Esp. However, the FP Filter considerably minimized this situation, reducing double identification of each plant, increasing process accuracy to identify and count individuals. This procedure reduced dependence on the previous step with regard to accuracy and can be used whenever necessary. For the combinations shown, the GCP did not generally increase the Paac value, either 'With GCP' (94.4 %) or 'Without GCP' (94.8 %); however, when analyzing F60 and F90 only, the gains were not significant. Accuracy and quality of cartographic products obtained with UAV are not restricted to the use of ground control points or georeferencing accuracy. Agisoft (2020), Ferrer-González et al. (2020), and Sanz-Ablanedo et al. (2018 state that accuracy and quality are the result of the choice and adjustment of several variables, such as aircraft type, sensor type, flight planning, amount of frontal and side overlap, flight altitude, quality of the image captured, modelling and camera calibration, SfM algorithms, and configuration of the algorithm parameters, et cetera. Another important and relevant factor are the surface characteristics in terms of shade texture and variations. Surfaces with random textures, with no uniform coloring, can facilitate the search for key points in the images, improving the SfM processing quality (Iglhaut et al., 2019;Remondino et al., 2014).
For all combinations, reduction in flight altitude (smaller GSD) increased accuracy to identify and count plants. The best performance was for combination F60 and 'Without GCP' ( Table 3). The best result was 97 % (Sb) for F60, without the use of GCP. However, as expected, the main error to identify plants was due to FN, which are plants with less height and a smaller crown. Data variability (ground truth) is also shown: a) average height 1.70 m and coefficient of variation (CV) of 22.0 %; b) average (circular) crown area 3.40 m 2 and CV 39.1 %. Fawcett et al. (2019) shows the heterogeneity effect when evaluating biometric attributes in younger plantations using an RGB sensor in UAV. The authors estimated accuracy when counting plants of the African oil palm (Elaeis guineensis) aged two, seven, and ten years and found 80.4 %, 98.2 %, and 94.9 % respectively.

Estimating plant height and volume
The ground control points (GCP) did not increase accuracy in the process to estimate plant height. In general, the  Agric. v.79, n.5, e20210052, 2022 values for RMSE and MAE 'Without GCP' were lower than 'With GCP' (Table 4).
MDS derived from the SfM processing is also subject to systematic and random errors introduced by GCP (James et al., 2017). During georeferencing point cloud, similarity transformation occurs in some of the parameters (translation, rotation, and scale) that can only compensate linear models (Agisoft, 2020). The authors point out that possible non-linear deformations can be removed by optimizing the parameters when editing point cloud and by calibrating the camera based on GCP coordinates. Nevertheless, the result was not sufficient to increase accuracy when using GCP.
On the other hand, the map scale had a direct influence on estimates of plant height, that is, lower GSD values provided better accuracy (lower RMSE and MAE) and greater correlation levels. Low-altitude flights significantly increase point cloud density, which improves 3D surface reconstruction (Seifert et al., 2019;Torres-Sánchez et al., 2018). Better results were achieved with F60, mainly for 'Without GCP'. The data derived from the UAV and SfM processing tend to underestimate plant height values ( Figure 5A), as reported by Dandois et al. (2015).
Certain factors contributed to underestimating plant height values: a) irregular crown architecture ( Figure 5B), which is still under formation; b) DSM resolution (GSD) not compatible with the size of some objects, such as branches and leaves, with higher GSD values tending to smooth the surface, reducing crown heterogeneity, and consequently reducing accuracy (Zarco-Tejada et al., 2014); c) point cloud precision (derived from the SfM), which is lower in plant crown compared to the terrain surface (Fawcett et al., 2019); and d) plant movement during acquisition, possibly due to the wind.
For relative height (DHM) and absolute altitude (DSM), better results were obtained with the DSM model, for instance, lower values for the RMSE and MAE. When using DHM with the methodology in this study (DHM = DSM -DTM), DTM is derived from two filters in the dense cloud, followed by classification, which must have smoothed the values for absolute height and consequently DHM and the plant height values (h D HM  ). The process of subtracting DTM from DSM is known and can result in underestimations by not representing the plant apex (Holman et al., 2016). This situation is fairly common in data derived from SfM processing and explains the systematic underestimation of plant height (Castro et al., 2019;Solvin et al., 2020). Based on this tendency, our study presented an alternative procedure using the developed routine, which succeeded in minimizing this effect, that is, when using the DSM only, the results were more satisfactory.
When estimating plant volume, the results above should be considered, particularly delimitation of the     Dong et al. (2020) in terms of the effects of different map scale (three flight altitudes) and lighting conditions (as occurred at flight F120) on estimating the biometric attributes of plants. They also complement recommendations of Tu et al. (2020), that is, testing different combinations of aerophotogrammetric parameters (Table 1) with other sensors (RGB) on different crops (citrus plants).
The GCP use is only necessary to carry out a temporal analysis of the crop of interest, since UAV with a navigation receiver does not have a differential correction system to provide sub-metric corrections. On the other hand, GCP can play an important role in areas with homogeneous surfaces, since they can improve the self-calibration process of the camera during alignment (Agisoft, 2020). Furthermore, 3D SfM models derived from aerial images may initially contain deformations or systematic errors. The use of ground control points can reduce these deformations (Eltner and Schneider, 2015;James and Robson, 2014). However, the user needs to be careful to the positional quality of the GCP, for instance, accuracy of geodetic coordinates, otherwise, when inaccurate coordinates are applied in SfM processing, a more complex error surface is introduced, rather than reducing the initial deformation (Sanz-Ablanedo et al., 2018).
Quantity and distribution of GCP on the ground have already been studied and reported in the literature. However, there are certain specific peculiarities in UAVbased photogrammetry (employing SfM processing) and because it mostly uses non-metric cameras and self-calibration, some results can be inconclusive or even contradictory, as reported Sanz-Ablanedo et al. (2018). Theoretically, when more GCP are used and well-distributed, the results are better for horizontal and vertical accuracy (Agüera-Veja et al., 2017;Ferrer-González et al., 2020;James et al., 2017;James and Robson, 2014). Sanz-Ablanedo et al. (2018) state that new studies reduce GCP dependence for SfM projects. Not using GCP saves time, eliminates georeferencing of photoidentifiable targets in the field, and reduces investments in the purchasing or renting a Geodetic Receiver. For Fawcett et al. (2019), low-cost alternatives to GNSS Receivers and techniques with accuracy less than 10 cm have proven promising, such as the technology Post Processed Kinematic (PPK), which can reduce GCP use (Zhang et al., 2019).
In terms of operational yield in the field, considering only the interval used for recording images during the flights, an effective field capacity of 22.4 ha h -1 (F60), 51.3 ha h -1 (F90) and 81.0 ha h -1 (F120) was obtained. Regarding SfM processing, Torres-Sánchez et al. (2018) state that lower flight altitude (lower GSD) increases the number of captured images and requires greater storage and processing time to generate the aerophotogrammetric products (Castro et al., 2019).

Performance of the FindCITRUS-V1 automated routine
In terms of computer processing time, using a Laptop with a Linux operating system (Ubuntu 18.04-bionic), i7-4500U CPU (4 core, 1.  The information is grouped into a personalized layout, with registry (Property, Assessed Tract, and Owner) and cartographic (DSM, DTM, and Orthophotomosaic) information, resulting in the Statistical Report Map (SRM). The report is generated in PDF format on a one-plant scale, which provides valuable information to monitor quality of 'Tahiti' acid lime orchards under formation, including statistical metrics associated to biometric attributes, which allow identification of management zones with different characteristics, providing greater assertiveness (Ampatzidis et al., 2019;and Ok and Ozdarici-Ok, 2017). It also allows creating high-resolution time series that can assist in irrigation management (Yeom et al., 2019).
The results are sufficient to define a method to obtain aerophotogrammetric data (using an RGB sensor in a UAV) that can be used to monitor quality of 'Tahiti' acid lime orchards under formation. The use of GCP is not necessary during the SfM geoprocessing stage and the map scale has strong correlation with the estimation quality of the attributes under evaluation. Moreover, our study brings the following contributions: (i) an automated and calibrated routine, using aspects of aerial photogrammetry and computer vision; (ii) routine estimations of various biometric attributes of the plant, analyzing them using descriptive statistics, also presenting them in a report (SRM) that can be used to analyze the quality of citrus groves under formation; (iii) routine requires minimal user intervention (a graphic and intuitive environment), allowing default parameters to be modified when necessary; (iv) routine was developed and could be used in a single open-source piece of software, QGIS; (v) routine presents excellent performance with ortho-photomosaics derived from RGB channels, that is, low-cost sensor that is easy to operate (Yao et al., 2019); (vi) assistance in phenotyping plants in the field; (vii) the method also provides the user (researchers, consultants, and producers) with a new tool to evaluate citrus orchards (in the formation process); and (viii) it is easily adapted to other perennial crops, such as mango, avocado, and coffee, among others.

Conclusion
The use of GCP during the SfM geoprocessing stage does not increase accuracy of estimating the biometric attributes under evaluation. The map scale (GSD size) is strongly correlated to the estimation quality of the attributes under evaluation, especially plant height. The routine developed shows high accuracy and excellent computational performance, in addition to an intuitive environment and is easily adapted to perennial crops.