High-Resolution UAV Imagery for Field Olive ( Olea europaea L.) Phenotyping

: Remote sensing techniques based on images acquired from unmanned aerial vehicles (UAVs) could represent an effective tool to speed up the data acquisition process in phenotyping trials and, consequently, to reduce the time and cost of the ﬁeld work. In this study, we assessed the ability of a UAV equipped with RGB-NIR cameras in highlighting differences in geometrical and spectral canopy characteristics between eight olive cultivars planted at different planting distances in a hedgerow olive orchard. The relationships between measured and estimated canopy height, projected canopy area and canopy volume were linear regardless of the different cultivars and planting distances (RMSE of 0.12 m, 0.44 m 2 and 0.68 m 3 , respectively). A good relationship (R 2 = 0.95) was found between the pruning mass material weighted on the ground and its volume estimated by aerial images. NDVI measured in February 2019 was related to fruit yield per tree measured in November 2018, whereas no relationships were observed with the fruit yield measured in November 2019 due to abiotic and biotic stresses that occurred before harvest. These results conﬁrm the reliability of UAV imagery and structure from motion techniques in estimating the olive geometrical canopy characteristics and suggest further potential applications of UAVs in early discrimination of yield efﬁciency between different cultivars and in estimating the pruning material volume.


Introduction
Olive groves are widespread in the Mediterranean basin, where 97% of total olive trees are grown [1]. Taking into account the new olive-producing countries (Australia, South Africa, Argentina, Chile and California) outside the Mediterranean Basin, olive trees can be considered one of the most representative tree cultivations in the world [1]. The diffusion of olive growing outside traditional areas of cultivation is mainly due to an increase in olive oil consumption over the last 20 years and to the availability of new planting systems (e.g., high-density and super-high-density) and more efficient orchard management strategies (irrigation and mechanized harvesting and pruning).
Currently, only a few cultivars (Arbequina, Arbosana, Chiquitita and Koroneiki) have the necessary low vigor to be suitable for super high density (SHD) olive groves [2][3][4][5]. In order to prevent possible loss of biodiversity linked to the adoption of SHD systems, the selection of local cultivars with the most desirable traits for these planting systems could represent a forward-looking strategy. However, while plant genotyping is nowadays an efficient process thanks to advances in genomics and biotechnology, phenotypic data are mainly collected by manual or visual methods [6,7].
Crop phenotyping is defined as the monitoring of morphological, physiological and biochemical traits resulting from the interaction between genetic and environmental factors [8,9]. Geometrical canopy characteristics and tree vigor are crucial traits in phenotyping studies in order to assess the cultivar's suitability to be cultivated in specific growing systems, as well as for assessing pruning practices [4,5,10,11]. Tree canopy shape and size are also related to fruit and oil quality as they affect light distribution within the canopy and, in turn, the amount of light intercepted by fruits [12][13][14]. Finally, the combined analysis of tree vigor data, based on both vegetation indices or canopy growth patterns, and tree crop load can be used to classify different olive cultivars according to their fruit yield efficiency.
The traditional manual-based measurements of geometrical canopy characteristics require intensive field work and experienced staff. The same measurements are also used to evaluate the impact of pruning on tree architecture. From the above, the need for time and cost-efficient methods to rapidly identify desirable genotypes is clear. High-throughput phenotyping technologies can play a crucial role in phenotyping trials thanks to their ability to monitor a large number of plants in a short period of time [7,15,16]. Recent studies focused on possible applications of ground and aerial platforms for plant phenotyping under field conditions [16][17][18][19][20][21][22]. Among them, unmanned aerial vehicles (UAVs) have been proposed as sensing platforms for field-based phenotyping due to their efficiency in producing a large amount of field scale information, great flexibility in flight scheduling and low operational costs.
In the few studies focused on the evaluation of UAVs for olive breeding purposes, high-resolution images have been used to infer geometrical canopy characteristics, such as tree height, canopy diameters and canopy volumes, but no information on the relationship between canopy spectral response and fruit yield of different genotypes are available [17,18,23,24]. Another possible application of UAVs in breeding trials concerns the evaluation of pruning effects on the tree canopy. In previous studies, the volume of the pruning material has been estimated as the difference in canopy volume measured before and after pruning [25,26]. To the best of our knowledge, there are no studies where the volume of pruning material has been directly estimated by UAV images.
The overall objective of this study was to assess the ability of UAVs equipped with RGB-NIR cameras in highlighting differences in geometrical and spectral canopy characteristics between eight olive cultivars in a super-high-density olive orchard. The relationship between spectral canopy characteristics and fruit yield was also investigated. Finally, we propose a new approach for direct pruning material estimation using a low-cost RGB camera.

Plant Material and Site Characteristics
The experiment was carried out in an irrigated olive orchard located near Sciacca in South-West Sicily (Italy) and planted in 2012 with 1-year-old self-rooted olive plants of eight different olive cultivars: Cerasuola, Nocellara del Belice, Biancolilla, Nocellara Etnea, Minuta, Calatina, Abunara and Koroneiki. The climate at the experimental site is characterized by hot and dry summers and mild winters. The soil was a sandy clay loam, consisting of 60% sand, 22% clay and 18% silt, with pH 7.7 and less than 5% active carbonates. The tree rows were oriented North-East to South-West. Distance between the rows was 4 m while three different spacing (2, 3 and 4 m) were used between trees along the rows. The trees were trained with a single central trunk with three main couples of branches, oriented along the row, starting at 0.7 m height above ground. A manual topping was performed to keep plant height below 3 m.

Field Data
On 27 February 2019, the canopy volume and the projected canopy area (PCA) were measured on trees of cvs. Biancollilla, Calatina and Koroneiki (5 trees for each cultivarplanting distance combination for a total of 45 trees) according to the following formulas: where D1 and D2 are crown diameters, Ht is tree height and Hb is the lowest canopy height from the ground. The same trees were manually pruned on the same day, after the canopy volume measurements. The pruned material from each tree was placed in the middle of the inter-row and weighed using a portable weight scale (Figure 1).

Image Processing Methods
The multispectral images were first mosaicked using Autopano Giga 3.5 Software (Kolor SARL, Challes-les-Eaux, France), then georeferenced and orthorectified using the ground-referenced points (ArcGIS software ® , ESRI, Redlands, CA, USA). Vegetative indices were calculated by means of the map algebra technique implemented in ArcGIS software (ESRI, Redlands, CA, USA). Normalized Difference Vegetation Index (NDVI), Green Normalized Difference Vegetation Index (GNDVI), Green Ratio Vegetation Index (GRVI), Radar Vegetation Index (RVI), Green Red Normalized Difference Vegetation Index (GRN-

Multispectral Imagery Acquisitions from the Unmanned Platform
The acquisition campaign was performed using an S1000 UAV octocopter (DJI, Shenzhen, China) able to fly autonomously over a predetermined waypoint course. The S1000 was equipped with a 2-axis stabilized gimbal equipped with a consumer photo camera (RGB) and a multispectral camera (NIR-RG). The RGB camera was a Coolpix P7700 (Nikon, Shinjuku, Japan) embodying a 12.2-megapixel CMOS sensor, whereas the NIR-RG multispectral camera was a Tetracam ADC-lite (Tetracam Inc., Gainesville, FL, USA) equipped with a 3.2-megapixel CMOS sensor. Images were recorded in the visible red (R), green (G) and near-infrared (NIR) domain with nominal wavelengths of 520-600, 630-690 and 760-900 nm, respectively. Images were acquired before (27 February 2019) and after (28 February 2019) pruning at noon under clear sky conditions. The flight altitude was 50 m above ground level (AGL), flying at 2.5 m s −1 speed. A total of 126 and 132 multispectral and RGB images, respectively, were acquired on 27 February, whereas 143 RGB images were acquired on 28 February. The image forward and side overlap (80% and 70%, respectively) guaranteed optimal photogrammetric processing. Geometrical parameters and vegetative indices used for cultivars characterization were calculated by processing the images acquired on 27 February (before pruning) in order to avoid the effects of canopy manipulation.
The three-dimensional canopy volume and the PCA of each tree were calculated following the same procedure reported in Caruso et al. [33]. Briefly, RGB images were processed using Agisoft PhotoScan Professional Edition (Agisoft LLC, St. Petersburg, Russia) for the generation of the 3D point clouds. The digital surface model (DSM) was obtained from the 3D point cloud and then processed in ArcGIS to obtain a digital terrain model (DTM). The normalized DSM, obtained by subtracting the DTM from the DSM, allowed to retrieve the height of each three-dimensional axes of the canopy point above the ground. The net tree canopy volume was calculated by subtracting the volume comprised between the ground level and 0.7 m (the height of the insertion of the first couple of branches) to the total volume of each tree [33]. The very high-resolution of the images and the absence of grass cover facilitated the separation of tree crowns from the background. The PCA was obtained using a height value threshold (0.7 m) which delimited soil and tree components. The canopy height profile of each tree along the row was obtained by extracting the pixel value (height above the ground level) from the normalized DSM at 0.1 m intervals. The mean canopy height profile was calculated as the average of 10 trees for each cultivar-planting distance combination. The volume of the pruned material of Biancolilla, Calatina and Koroneiki trees was estimated following the same procedure used for the canopy volume estimation. The entire volume (from ground level to the top of the pruning heap) was considered representative of the pruning material volume (Figure 1).

Statistical Analysis
The experimental design was a split-plot where cultivar was considered as the main plot and planting distance as the sub-plot. Each combination (cultivar-planting distance) consisted of 10 (canopy volume, PCA and canopy height profile) and 7 (fruit yield and TCSA) replicate trees. Means were separated by honestly significant differences (Tukey's HSD) at p < 0.05 after analysis of variance (ANOVA). Where applicable, data were analyzed by regression using Costat (CoHort Software, Monterey, CA, USA). The root mean square error (RMSE) and the mean absolute percentage error (MAPE) of residuals were calculated along with the regression fit and the squared correlation coefficient between measured and estimated tree parameters. Discriminant analysis was performed using JMP (JMP SAS, Campus Drive Cary, NC, USA) to evaluate the vegetation indices derived by aerial images as a set of independent variables for discriminating the eight cultivars (Cerasuola, Nocellara del Belice, Nocellara Etnea, Minuta, Calatina, Abunara, Koroneiki). All the tests used to evaluate the analysis (Wilks' Lambda, Pillai's Trace, Hotelling-Lawley, Roy's max root) were significant (p < 0.0001). Canonical details on discriminant analysis are provided as supporting materials (Table S1).

Results
The relationships between measured and estimated canopy height, projected canopy area and canopy volume were linear regardless of the different cultivars and planting distances (slope of 0.89, 0.96 and 0.85 for canopy height, projected canopy area and canopy volume, respectively) (Table 1, Figure 2). Good correlation (R 2 equal to 0.79, 0.78 and 0.82 for canopy height, projected canopy area and canopy volume, respectively) and accuracy levels (RMSE equal to 0.12 m, 0.44 m 2 and 0.68 m 3 for canopy height, projected canopy area and canopy volume, respectively) were also observed (Table 1, Figure 2). The geometrical canopy characteristics derived by RGB images showed significant differences between cultivars and between planting distances ( Table 2, Table S2). The biggest canopies were those of cultivar Cerasuola (7.30 m 3 , average of all the sampling distances), whereas the smallest ones were measured in Calatina trees (5.27 m 3 ). The projected canopy area was similar between cultivars (4.07-4.27 m 2 ) except for the cultivar Calatina, which showed significantly lower values (3.29 m 2 ) than the others ( Table 2). ent in canopy height ( Figure 3). Differences in the shape of the canopy profile became more evident in trees planted at 4 × 3 and 4 × 4 m (Figures 4 and 5). The percentage of the vegetative wall filled by tree canopies was calculated by assuming a complete and uniform canopy development below the maximum canopy height profile. Trees planted at row distances of 2 and 3 m showed similar percentages of filled canopy walls between each other, whereas lower values were measured in olives planted at 4 × 4 m (Figures 3-5).        Pruning material mass varied according to tree cultivar (11.8, 7.21 and 4.32 kg pruning material per tree-an average of 15 trees per cultivar-in Koroneiki, Biancoli and Calatina trees, respectively). The combined use of RGB images by UAV and structu from motion techniques allowed to properly estimate the mass removed by pruning fro each tree ( Figure 6). A good relationship (R 2 = 0.95) was found between the pruning ma material weighted on the ground and its volume estimated by UAV ( Figure 6A). Similar the reduction in canopy volume (canopy volume before pruning minus canopy volum Pruning material mass varied according to tree cultivar (11.8, 7.21 and 4.32 kg of pruning material per tree-an average of 15 trees per cultivar-in Koroneiki, Biancolilla and Calatina trees, respectively). The combined use of RGB images by UAV and structure from motion techniques allowed to properly estimate the mass removed by pruning from each tree (Figure 6). A good relationship (R 2 = 0.95) was found between the pruning mass material weighted on the ground and its volume estimated by UAV ( Figure 6A). Similarly, the reduction in canopy volume (canopy volume before pruning minus canopy volume after pruning) estimated by UAV was correlated (R 2 = 0.84) with the pruned mass ( Figure 6B). The impact of different pruning intensities on the projected canopy area was less evident (R 2 = 0.63) ( Figure 6C).
Horticulturae 2021, 7, x FOR PEER REVIEW 11 of 18 after pruning) estimated by UAV was correlated (R 2 = 0.84) with the pruned mass ( Figure  6B). The impact of different pruning intensities on the projected canopy area was less evident (R 2 = 0.63) ( Figure 6C).  A positive linear relationship between NDVI measured on February 2019 and the annual TCSA increment was found regardless of cultivar and planting distances ( Figure 7A-C). Fruit yield per tree ranged between 4.8 kg (Nocellara del Belice 4 × 2 m) and 26.3 kg (Cerasuola 4 × 4 m) in 2018 and between 2.1 kg (Abunara 4 × 2) and 25.1 kg (Biancolilla 4 × 4) in 2019 ( Figure 7D-I). In 2019, fruit yield was strongly affected by both abiotic and biotic stresses. Particularly, hot winds and rain at blooming reduced fruit set, whereas olive fruit fly (Bactrocera oleae Rossi) caused a consistent fruit drop before harvest. The NDVI measured in February 2019 was positively related to fruit yield measured a few months earlier (November 2018) with the only exception of Calatina, which did not follow the general trend of the other cultivars ( Figure 7G-I). In 2019, yields were very highly variable depending on the intensity of biotic and abiotic stresses effects, and, as expected, no relationship between NDVI and tree productivity emerged ( Figure 7D-F).  Figure 7D-I). In 2019, fruit yield was strongly affected by both abiotic and biotic stresses. Particularly, hot winds and rain at blooming reduced fruit set, whereas olive fruit fly (Bactrocera oleae Rossi) caused a consistent fruit drop before harvest. The NDVI measured in February 2019 was positively related to fruit yield measured a few months earlier (November 2018) with the only exception of Calatina, which did not follow the general trend of the other cultivars ( Figure 7G-I). In 2019, yields were very highly variable depending on the intensity of biotic and abiotic stresses effects, and, as expected, no relationship between NDVI and tree productivity emerged ( Figure 7D-F).  Figure 8 shows the results of the discriminant analysis performed using six vegetation indices derived by the green (G), red (R) e near-infrared (NIR) bands from the multispectral images. Discriminant analysis, performed including all planting distances, was able to separate the different cultivars using the six vegetation indices. In particular, GRVI and RVI were able to discriminate Minuta and Calatina toward Canonical 1, while Cerasuola and Nocellara del Belice were mainly discriminated by GNDVI (Figure 8). Koroneiki was separated on Canonical 2 by NDVI and MSAVI indices. annual TCSA increment = 5.7 NDVI -2.5, R 2 = 0.69 (A); annual TCSA increment = 6.4 NDVI -2.7, R 2 = 0.66 (B); annual TCSA increment = 8.7 NDVI -3.8, R 2 = 0.74 (C); fruit yield 2018 = 287 NDVI -136, R 2 = 0.75 (G); fruit yield 2018 = 217 NDVI -98 R 2 = 0.65 (H); fruit yield 2018 = 335 NDVI -159 R 2 = 0.67 (I). Figure 8 shows the results of the discriminant analysis performed using six vegetation indices derived by the green (G), red (R) e near-infrared (NIR) bands from the multispectral images. Discriminant analysis, performed including all planting distances, was able to separate the different cultivars using the six vegetation indices. In particular, GRVI and RVI were able to discriminate Minuta and Calatina toward Canonical 1, while Cerasuola and Nocellara del Belice were mainly discriminated by GNDVI (Figure 8). Koroneiki was separated on Canonical 2 by NDVI and MSAVI indices.

Discussion
Previous studies ascertained the ability of UAV and image analysis in geometrical canopy characteristics estimation in orchards and vineyards [22,[33][34][35][36][37]. In this study, geometrical canopy characteristics, such as canopy height, projected canopy area and canopy volume, were satisfactorily estimated (MAPE of 3.26, 10.40 and 9.25%, respectively) following the methodology proposed by Caruso et al. [33]. In a similar study carried out on young (27 months old) olive trees, the relative error between estimated and measured canopy volume accounted for 17.65% [17]. The significant differences in canopy volume and in projected canopy area between cultivars were also observed in previous experiments [17,18,23], confirming their geometrical canopy peculiarities observed using an on-ground measurement approach [38]. The analysis of the canopy height profile provides useful information about the ability and timing of different olive cultivars in filling a predetermined leaf area wall (in our experiment 4.6, 6.9 and 9.2 m 2 , for 4 × 2, 4 × 3 and 4 × 4 m planting systems, respectively). One of the assumptions in our procedure (a complete and uniform canopy development below the maximum canopy height profile, Figures 3-5) represents a simplification that cannot be applied a priori under different experimental conditions. However, our results are in agreement with those observed in previous experiments [39,40], suggesting a possible use of this technique under specific site conditions. Geometrical canopy characteristics are also important in order to optimize the spray volume in pesticides application [41,42]. Different dosing models are proposed to calculate the optimal amount of pesticide to be applied, including the Crown Height model (CHT model), the Surface Orchard model (SO model) and the Tree Row Volume model (TRV model) [43]. In this regard, the ability of UAVs in easily providing these canopy features could positively contribute to reducing the environmental impact of pesticide applications.
In previous studies, the pruning biomass was estimated by comparing the tree canopy volume estimated by UAV before and after pruning [25,26,44,45]. To the best of our knowledge, this is the first study where UAV imagery was used to directly estimate the pruning material volume. In our research, following the same approach proposed for the canopy volume estimation [33], we estimated the volume of the pruning heap placed by pruners on the ground. The relationship between the pruning material mass and its volume estimated by UAV was very good (R 2 = 0.95). These results could represent a substantial improvement in data collection in phenotyping trials where intensive field work is needed to acquire this information at tree level on the ground. It is important to point out that these relationships are "pruners"-specific because the way through which the operator places and arranges the pruning material on the ground could affect its final volume. Moreover, in order to be clearly visible from the UAV, the pruning material should be placed in the middle of the inter-row. The above specifications only slightly reduce the applicability of this methodology because, in olive trials, pruning is usually carried out always by the same pruners. The amount of pruning material was also consistent with the reduction in canopy volume after pruning. Some discrepancies between the pruning mass and the corresponding reduction in canopy volume may be due to the position of the shoots and branches to be removed. Inner shoots and branches have little impact on the canopy shape and, therefore, their removal may have only a slight (or none) effect on the canopy volume. Johansen et al. [26] used multispectral UAV images to assess changes in lychee trees structures (tree crown perimeter, width, height, area and Plant Projective Cover, PPC) induced by pruning. They reported that all measured tree structural parameters collected after a pruning event were lower than those measured before pruning. Similar results were found in a previous study where similar technologies were used to quantify the impact of different pruning strategies (traditional, adapted and mechanical pruning) on olive tree architecture and annual canopy growth [25].
The olive cultivars evaluated in this study were different for canopy spectral response. While morphological markers are widely used for olive cultivars identification [46][47][48], little information is reported about leaf and canopy reflectance [49,50]. The combined use of leaf hyperspectral reflectance and data-mining techniques was successfully applied to discriminate 10 olive cultivars in an experiment carried out in Portugal [50]. In another experiment, Avola et al. [49] used a UAV equipped with a multispectral camera (520-600, 630-690 and 760-900 bands) to discriminate field-grown olive cultivars (Frantoio and Leccino) in an experimental plantation with different scion/rootstock combinations. The same authors were able to effectively discern the two scions, whereas the effect of rootstock on canopy reflectance appeared unclear. Using a similar approach, we were able to highlight differences in spectral response between the eight olive cultivars evaluated.
A positive linear relationship between NDVI and fruit yield was found for all cultivars, with the only exception of Calatina trees. This result is not surprising, because Calatina has been previously described as a low vigor but highly productive cultivar [38]. In this sense, the integration of data derived by remote sensing (tree vigour-NDVI) and field measurements (fruit yield per tree) could improve the characterization of olive cultivars according to their vegetative-to-reproductive balance. This is a crucial information in phenotyping trials which are mainly aimed at finding low vigor but highly productive olive varieties. Recent studies focused on olive yield forecasting based on geometrical and spectral canopy characteristics [51,52]. Sola-Guirado et al. [51], in an experiment carried out in Greece, developed a multiple linear regression based on NDVI, canopy volume and the average slope of the tree location, yielding an R 2 of 0.6. Although good performance of vegetative indices in discriminating olive cultivars and in yield forecasting was found in our study and in previous ones [49][50][51][52], it is important to pay attention to possible limitations in practical transferring of these techniques. For instance, the ability of vegetative indices in cultivar identification has been tested in experimental fields characterized by homogeneous pedo-climatic conditions and agronomic practices within the orchard [49, this experiment]. Under uncontrolled conditions (larger areas or abandoned olive orchards), abiotic and biotic stresses should induce changes in vegetative indices values greater than those induced by genotype. Similarly, forecasting fruit yield based only on NDVI values and canopy geometrical characteristics can be risky because biological parameters affecting production are not taken into account. In our experiment, we experienced this issue in 2019, when we did not find any relationship between the NDVI measured in February and the final fruit yield due to the abiotic and biotic stresses that occurred during the growing season. In general, unfavorable events affecting flowering, fruit-set or ripening (e.g., spring frost, heat waves, olive fruit fly) can dramatically modify the final fruit yield per tree but induce negligible effects on tree canopy volume and vigor (NDVI).
Overall, our results confirm the reliability of UAV imagery and structure from motion techniques in estimating the olive geometrical canopy characteristics such as canopy volume and projected canopy area. The relationship between NDVI and fruit yield suggests the usefulness of the combined use of field and remotely sensed data for the early discrimination of yield efficiency between different cultivars. The possibility to remotely estimate the pruning material volume may dramatically speed up the acquisition process of this information and, consequently, reduces the time and cost of the field work. Finally, in a commercial farm context, the UAV-derived information could provide olive growers the opportunity to take corrective actions to improve orchard management activities such as pruning and pesticide applications.
Author Contributions: G.C.: conceptualization, investigation, methodology, formal analysis, resources, software and writing-original draft. G.P.: investigation, data curation, formal analysis, software and writing-original draft. F.P.M.: formal analysis, funding acquisition, validation and writing-review and editing. T.C.: conceptualization, supervision, validation and writing-review and editing. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.