Asier Uribeetxebarria
email: uribeetxebarria.asier@gmail.com

J A Martínez Casasnovas

Bruno Tisseyre

Serge Guillaume

Alexandre Escolà

J R Rosell Polo

J Arno
email: jarno@eagrof.udl.cat

José A Martínez-Casasnovas

Joan R Rosell-Polo

Jaume Arnó

Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards

Keywords: Ranked set sampling, Fruit growing, Sampling efficiency, Yield estimation Digital multispectral camera ECa Soil apparent electrical conductivity MTLS UAV Unmanned aerial vehicle

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Fruit tree crops are considered high value products [START_REF] Aggelopoulou | Spatial variation in yield and quality in a small apple orchard[END_REF]), but they demand more tasks from farmers compared to arable crops. Among the tasks required in peach orchards, harvesting is a complex process that, necessarily done by hand, usually also needs to be finished in a short time window. Because of that, a large amount of labor is often needed once the quality standards demanded by the market have been reached (ripening of the fruit). Hence, obtaining a reliable yield estimate in advance would improve the logistics of the entire process [START_REF] Wulfsohn | Multilevel systematic sampling to estimate total fruit number for yield forecasts[END_REF], also contributing to anticipate operational costs more effectively.

In perennial crops, it is known that regular patterns are adopted by planting commercial plots with genetically uniform plant material [START_REF] Miranda | Sampling stratification using aerial imagery to estimate fruit load and hail damage in nectarine trees[END_REF]. Therefore, it is expected to achieve regular growths and certain homogeneity in production and quality. From this point of view, making an estimate of the yield should not raise any issue. By randomly choosing a few trees, the estimation of yield should be quite accurate knowing the total number of trees within the plot. However, perennial plots are not always homogeneous and spatial variability is present POSTPRINT of the article: Uribeetxebarria, A., Martínez-Casasnovas, J.A., Tisseyre, B., Guillaume, S., Escolà, A., Rosell-Polo, J.R., Arnó, J. 2019. Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards. Comput. Electron. Agric. DOI: https://doi.org/10. 1016/j.compag.2019.104931 4 due to environmental factors, as it has been shown in different studies [START_REF] Berman | Water stress and crop load effects on fruit fresh and dry weights in peach (Prunus persica)[END_REF][START_REF] Taylor | A comparison of the spatial variability of vineyard yield in European and Australian production systems[END_REF][START_REF] Arnó | Spatial variability in grape yield and quality influenced by soil and crop nutrition characteristics[END_REF]. Moreover, since many of the factors that affect yield (or fruit load) are spatial dependent, yield spatial distribution within the orchards usually presents patterns and is not random [START_REF] Aggelopoulou | Delineation of management zones in an apple orchard in Greece using a multivariate approach[END_REF]. Under this scenario of structured variability with yield potential varying for each tree and for the different productive areas within the plot, simple random sampling (SRS) strategy may be inefficient for estimating average crop loading per tree. In order to fulfill fruit grower demands, yield forecast should be made by using 5-6 sampling trees and with maximum error of 10%. Therefore, SRS may not be the best option since it does not take into account any spatial organization of the variable to estimate. Moreover, a large sample size may be needed to fulfill growers' constraints [START_REF] Carrillo | Use of multi-spectral airborne imagery to improve yield sampling in viticulture[END_REF]. Only the simplicity to implement and understand the results of SRS would explain why this method is still used in agriculture in forecasting tasks. Therefore, the Precision Agriculture concept allows this issue to be overcome since, as a management strategy, it is about gathers, processes and analyzes temporal, spatial and individual data and combines it with other information to support management decisions according to estimated variability for improved resource use efficiency, productivity, quality and sustainability of agricultural production (ISPA Newsletter, 2019).

Although it is fair to mention that one of the pioneering work on sampling in fruit trees was published in the forties [START_REF] Pearce | Sampling methods for the measurement of fruit crops[END_REF], sampling methods making use of ancillary information provided by high spatial resolution sensors have been proposed in recent years. Previous analysis of the within-field variability of such information is often key when designing the sampling strategy. Even spatial variability on a smaller scale has also been the option used by other authors. So in 1990, Monestiez et al. proposed a geostatistical approach to optimally sample within the trees to improve yield estimates. Wulfshon et al. (2012) discussed using POSTPRINT of the article: Uribeetxebarria, A., Martínez-Casasnovas, J.A., Tisseyre, B., Guillaume, S., Escolà, A., Rosell-Polo, J.R., Arnó, J. 2019. Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards. Comput. Electron. Agric. DOI: https://doi.org/10.1016/j.compag.2019.104931 5 images to optimize sample sizes and opted for a multilevel systematic sampling to estimate the number of fruits reaching errors of just 10%, then extending this same method to estimate quality parameters [START_REF] Martinez | Using multilevel systematic sampling to study apple fruit (Malus domestica Borkh.) quality and its variability at the orchard scale[END_REF]. More recently, stratified sampling based on the use of auxiliary information for the spatial stratification of the sample has also shown promising results. Examples of this research are shown in [START_REF] Arnó | Comparing efficiency of different sampling schemes to estimate yield and quality parameters in fruit orchards[END_REF] and [START_REF] Miranda | Sampling stratification using aerial imagery to estimate fruit load and hail damage in nectarine trees[END_REF][START_REF] Miranda | Sampling stratification using aerial imagery to estimate fruit load in peach tree orchards[END_REF]. Specifically, trees were stratified using clustered NDVI-based aerial images (NDVI-Normalized Difference Vegetation Index), either alone or in combination with other information layers such as the trunk cross-sectional area (TCSA). In fact, because of the increasingly widespread use of remote sensing in fruit growing, there are companies that already offer advice for yield estimation based on the use of images in stratified sampling schemes. This practice manages to improve efficiency compared to random sampling, although sample sizes are still needed above what would be desirable. In addition, obtaining the representative NDVI value for an individual canopy or, even more, a manual measurement of the trees is not always a simple task [START_REF] Fountas | Site-specific management in an olive tree plantation[END_REF].

The use of ancillary information to stratify sampling has been favored by the availability and rapid access to proximal and remote sensors data. Among the sources of information most used for this purpose, it is worth mentioning the NDVI images [START_REF] Meyers | Research Note. Use of Normalized Difference Vegetation Index Images to Optimize Vineyard Sampling Protocols[END_REF], the apparent electrical conductivity (ECa) of the soil [START_REF] Mann | Delineating productivity zones in a citrus grove using citrus production, tree growth and temporally stable soil data[END_REF][START_REF] Arnó | Comparing efficiency of different sampling schemes to estimate yield and quality parameters in fruit orchards[END_REF], and time series of yield maps provided by manual harvesting or monitored by sensors on harvesters [START_REF] Araya-Alman | Using ancillary yield data to improve sampling and grape yield estimation of the current season[END_REF]. Other data are also available while other technologies for acquiring crop data have been fine-tuned. It refers to RGB and multispectral cameras (Ulzii-Orshikh et al.

2017

), mounted on terrestrial platforms or unmanned aerial vehicles (UAV), and mobile terrestrial laser scanners (MTLS) based on 2D LiDAR sensors [START_REF] Rosell | Obtaining the three-dimensional structure of tree orchards from remote 2D terrestrial LIDAR scanning[END_REF][START_REF] Escolà | Mobile terrestrial laser scanner applications in precision fruticulture/horticulture and tools to extract information from canopy point clouds[END_REF] POSTPRINT of the article: Uribeetxebarria, A., Martínez-Casasnovas, J.A., Tisseyre, B., Guillaume, S., Escolà, A., Rosell-Polo, J.R., Arnó, J. 2019. Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards. Comput. Electron. Agric. DOI: https://doi.org/10.1016/j.compag.2019.104931 6 2017), offering new possibilities for capturing canopy information from fruit trees. These higher resolution data allow the trees to be characterized in greater detail, but possible advantages of their use in sampling have not yet been verified.

As already mentioned, SRS is a well-known method in agriculture. It is statistically consistent, but not always the most efficient. Moreover, there is a growing interest in developing sampling methods that, in addition to provide accurate and unbiased yield estimates, allow small sample sizes to be used. Stratified sampling [START_REF] Cochran | Sampling Techniques[END_REF]) may be an option, because it fulfills the characteristics described above. However, to stratify the samples based on ancillary variable maps, fruit growers have to define the strata as a preliminary step by choosing between different classifications and zoning techniques. How many strata, how the strata are obtained, and how sampling points should be allocated in each stratum are decisions to be made. Therefore, this combination of subjective factors probably ends up affecting the goodness of yield estimates.

Another option is using the Ranked Set Sampling (RSS) [START_REF] Mcintyre | A method for unbiased selective sampling, using ranked sets[END_REF]. This method is interesting from the operational point of view because it does not need these preliminary steps.

To obtain the sample, once the number of sites to sample is set, RSS only requires applying a ranking mechanism based on the distribution of an auxiliary variable [START_REF] Wolfe | Ranked set sampling[END_REF]. This ranking mechanism allows the trees (sampling points), previously taken randomly in a first iteration from the population, to be ranked from lowest to highest, according to the value that in these trees takes an ancillary variable. More details about the ranking mechanism and the different iterations of the whole process can be found in the corresponding section below (Basics of the sampling methods). With everything, the ancillary information is therefore key in the procedure since it has to present high correlation with the variable to estimate (the fruit POSTPRINT of the article: Uribeetxebarria, A., Martínez-Casasnovas, J.A., Tisseyre, B., Guillaume, S., Escolà, A., Rosell-Polo, J.R., Arnó, J. 2019. Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards. Comput. Electron. Agric. DOI: https://doi.org/10.1016/j.compag.2019.104931 7 load, in the case of this work). In short, finding the best ancillary information to perform RSS in fruit orchards is a pending issue to solve.

Focusing our research on comparing RSS and SRS, the main objective of the present work was to evaluate whether the use of ranked samples significantly improves fruit load estimates in fruit growing. To do this, two main issues were addressed: i) to determine the most appropriate ancillary information to be used in the ranking process, and ii) to assess the efficiency in terms of coefficient of error when comparing RSS to SRS for different sample sizes.

Materials and methods

Study plot

The study was carried out in a commercial peach orchard located in Aitona (45° 94' 71" N, 0° 29' 20" E, ETRS89, 160 m a.s.l., Lleida, Catalonia, Spain). The plot covered an area of 2.24 ha. It was planted in 2012 with white peach (Prunus persica cv. 'Platycarpa') according to a regular plantation pattern of 5x2 m (Fig. 1). The total number of trees was 1816, which were delimited and counted using an airborne image. Peach trees were planted in form of "Catalan" vessel shape, which is the most common training system in Catalonia. This canopy management system produces some visible gaps between trees within the same row, being relatively easy to individualize each tree. Regarding agricultural tasks, the plot was managed like many of the commercial orchards in the region.

Drip irrigation system was used for water and fertilizers supply.

The plot was representative of the so-called Ebro depression area, characterized by the presence of broad flat regions resulting from infilled valley bottoms and residual landforms. Climate is typical of hot semi-arid areas, with strong seasonal temperature variations and an annual rainfall frequently below 400 mm. As a consequence, soils in the region usually have high concentrations of calcium carbonate (CaCO3) and low content of organic matter (OM).

However, as a main feature, within-field soil heterogeneity is prevalent in many plots due to POSTPRINT of the article: Uribeetxebarria, A., Martínez-Casasnovas, J.A., Tisseyre, B., Guillaume, S., Escolà, A., Rosell-Polo, J.R., Arnó, J. 2019. Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards. Comput. Electron. Agric. DOI: https://doi.org/10.1016/j.compag.2019.104931 9 the levelling of original terraces in the eighties to facilitate agricultural mechanization. The plot under study resulted from this type of transformation. As for the spatial variability of trees within the plot, individualized canopy projected area values for each tree were obtained (see section 2.4.2 for more details). The average value was 2.82 m 2 , with a range between 0.62 m 2 and 5.50 m 2 , and a Coefficient of Variation (CV) close to 21%, which can be considered a lowto-moderate variability level for fruit orchards.

Systematic sampling of the orchard

A regular sampling grid of 15x15 m was defined in 2016 to sample the plot systematically (Fig. 1). The trees coinciding with the grid nodes were taken as sampling points. Grid size was previously defined by variographic analysis of the apparent electrical conductivity (ECa) data provided by a Veris 3100 soil sensor (Veris Technologies, Inc., Salinas, KS, USA) from an onthe-go survey of the plot. Once a strong autocorrelation of the ECa data was verified, half of the variogram range (ArcGIS 10.4.1 software) was taken as the optimum distance between sampling points [START_REF] Kerry | Sampling in Precision Agriculture[END_REF]. On the other hand, to avoid border effect, no sample sites were positioned in a buffer zone of 15 meters from the limits that separated other plots (short sides of the rectangular shape of the plot). However, this buffer zone was not applied on the sides alongside with other fields because the cultivated peach variety and the canopy training system were the same in these adjacent plots.

In total, 104 trees (sampling points) were defined (Fig. 1), and the number of fruits was counted manually in each tree four weeks before harvest. As proposed by [START_REF] Miranda | A statistical model to estimate potential yields in peach before bloom[END_REF], yield in fruit orchards is usually measured in terms of fruit load (number of fruits per tree). As POSTPRINT of the article: Uribeetxebarria, A., Martínez-Casasnovas, J.A., Tisseyre, B., Guillaume, S., Escolà, A., Rosell-Polo, J.R., Arnó, J. 2019. Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards. Comput. Electron. Agric. DOI: https://doi.org/10.1016/j.compag.2019.104931 10 yield (kg/tree) basically depends on fruit load, the latter was the variable used in this work.

Covering the entire plot, these 104 trees were taken as a reference population on which to compare two sampling methods.

Basics of the sampling methods

Simple random sampling (SRS)

The size of the sample was initially set at N = 12 trees or sampling points. Randomly selected without replacement, each of the twelve measured fruit load values could be considered representative of the population (or parcel under study), and with identical probability of choice.

Therefore, the mean of the sample allowed the average number of fruits per tree to be estimated without bias. However, there was a possibility that such a mean did not provide a truly representative picture of the orchard mean (Wolfe 2010).

Ranked set sampling (RSS)

Ranked set sampling is not a new method. Initially developed by [START_REF] Mcintyre | A method for unbiased selective sampling, using ranked sets[END_REF], the method was proved to be effective in improving the efficiency of pasture sampling. In this work, making use of adequate ancillary information (necessary for the ranking mechanism as it has been previously mentioned), the method was updated to evaluate its use in fruit growing.

Briefly, the method is based on an iterative sampling process [START_REF] Wolfe | Ranked set sampling[END_REF]. Following the procedure for a sample size N = 12, a first sample of size k = 12 was taken randomly without POSTPRINT of the article: Uribeetxebarria, A., Martínez-Casasnovas, J.A., Tisseyre, B., Guillaume, S., Escolà, A., Rosell-Polo, J.R., Arnó, J. 2019. Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards. Comput. Electron. Agric. DOI: https://doi.org/10.1016/j.compag.2019.104931 11 replacement from the 104 trees initially marked within the plot, and then the 12 selected trees were arranged according to some attribute (ranking mechanism). This process of sorting the sampled trees allowed a ranking to be established, usually from the lowest to the highest presence of the selected attribute. In our case, the attribute used in the ranking mechanism was an ancillary variable (or secondary variable) that, being related to the number of fruits per tree (as a relevant component of yield), should also be easy to measure and at an affordable cost.

The individual (tree) with the lowest value in this ranking was finally included as the first item in the ranked set sample, and the property of interest (fruit load) was formally considered only for this unit and denoted by 𝑌 [?] . The remaining k-1 observations were discarded in this without replacement strategy first sampling step.

In a second iteration, another sample of k = 12 observations was again selected without replacement and ranked following the same procedure as before. The individual (tree) ranked as the second smallest attending the attribute (ancillary variable) was now chosen, providing a second number of fruits value 𝑌 [@] that was added to the sample. The process was repeated until obtaining N = 12 measured observations 𝑌 [?] , 𝑌 [@] , ⋯ , 𝑌 [B] that constitutes a balanced ranked set sample of size N.

The mean of the N values (ranked set sample mean) is, like SRS, an unbiased estimate of the population mean [START_REF] Takahasi | On unbiased estimates of the population mean based on the sample stratified by means of ordering[END_REF]). However, having applied a ranking mechanism, each of the individual sample items in a balanced RSS represented a very distinctly different portion of the underlying population. Figure 2 represents this situation. The histogram represents the actual distribution of the number of fruits per tree (fruit load) in the plot under study (approximately, normal). Obtaining a sample through SRS involves sampling within that POSTPRINT of the article: Uribeetxebarria, A., Martínez-Casasnovas, J.A., Tisseyre, B., Guillaume, S., Escolà, A., Rosell-Polo, J.R., Arnó, J. 2019. Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards. Comput. Electron. Agric. DOI: https://doi.org/10.1016/j.compag.2019.104931 12 normal population (red curve in Fig. 2), and non-representative sample values may be drawn.

Instead, using a perfect ranking mechanism (since the ancillary variable that provides the attribute is well correlated with the fruit load), each of the observations included in the RSS behave like mutually independent order statistics, with densities that can be represented by the individual marginal density curves in Figure 2 (in this case, k = 6). Therefore, because of this extra structure provided by the judgment ranking (Wolfe 2004), RSS observations were much more likely to represent the full range of variation to estimate the population mean in a more effective way. This procedure is probably better than or equivalent to stratifying according to the probability distribution of the ancillary variable (stratified sampling based on the distribution of the percentiles). Adopting the previous notation for a sample size N = 12, the mean was calculated using (1),

𝑌 CDD = ? B ∑ 𝑌 [%] B %F?
.

(1) Regarding the variance, Wolfe (2010) provided the formula (2):

𝑉𝑎𝑟"𝑌 CDD ' = 𝑉𝑎𝑟"𝑌 DCD ' - ? B L ∑ "𝜇 [%] * -𝜇' @ B %F?
(2)

where

𝜇 [%] * = 𝐸"𝑌 [%] ', for 𝑗 = 1, … , 𝑁. Since ∑ "𝜇 [%] * -𝜇' @ ≥ 0 B %F?
, it follows from Eq. ( 2) that the variance of the RSS mean (𝑌 CDD ) is, at most, equal to the variance of the SRS mean (𝑌 DCD ).

POSTPRINT of the article: * is a key in improving the efficiency of the RSS compared to the SRS, providing the ancillary variable is highly correlated with the fruit load.

Ancillary variables in the ranking mechanism

In this section, the most relevant characteristics of the ancillary variables used to rank sampled trees are shown. Both terrestrial and remote sensors, whether used or recently introduced in the framework of precision agriculture, were considered to provide these data. Table 1 shows the platforms and sensors tested with additional indication of the ancillary information provided as well as a first assessment of the technical difficulty to obtain and process the data.

POSTPRINT of the article: Germany). The UAV was handled in manual configuration at constant speed of 18 km•h -1 and constant altitude up to 100 m. Pitch and roll movements were minimized using MK HiSight SLR2 gimbal support equipped with two servo motors. The camera was configured to take an image every two seconds during the flight, being the autofocus activated and the exposure fixed.

Nadir images were pre-selected to obtain an orthomosaic image with a spatial resolution of 2 cm per pixel using Agisoft Photoscan Professional software (Agisoft LLC, St. Petersburg, Russia).

Regarding the airplane image, a 4-band digital multispectral camera (DMSC) was used allowing images centred at 450 nm (blue), 550 nm (green), 675 nm (red) and 780 nm (near infrared) to be acquired. The images were pre-processed by SpecTerra (SpecTerra Services Pty In both cases, images were processed by manually delimiting each of the 104 individual trees sampled within the plot. The polygons used to delimit the projected canopies to the ground allowed the area and the perimeter of each tree canopy to be calculated. ancillary information was denoted with the letter H (Table 1). As with the previous maps for the other sensors, ArcMap 10.4.1 (ESRI, Redlands, CA, USA) performed spatial data analysis.

On the other hand, it is worth mentioning that the single information of the NDVI could be inconsistent with the expected fruit load of a tree. The bigger the canopy, the higher the expected fruit load [START_REF] Mann | Delineating productivity zones in a citrus grove using citrus production, tree growth and temporally stable soil data[END_REF] 19

𝑁𝐷𝑉𝐼 Q = 𝑁𝐷𝑉𝐼 × RQS RQS TUV (3)
where NDVI and NDVIC were the initial value and the corrected value of the NDVI, respectively, TCA was the tree canopy projected area obtained from the airplane image, and

TCAmax was the largest tree canopy projected area within the plot.

Apart from the remote sensor, a terrestrial OptRx sensor (Ag Leader Technology, Ames, IA, USA) was also used to acquire canopy reflectance data using an all-terrain vehicle (ATV) as a platform. Since the readings were taken sideways passing through the alleyways of the orchard, data referred to the lateral reflectance of the canopy. Specifically, the sensor was maintained at a distance of approximately 1 m from the canopy and moved at a constant speed of 5 km•h -1 .

With the OptRx, only one side of each row was measured. Three different wavelengths were acquired every second (1 Hz) in the ranges of red (670 nm), red edge (728 nm), and near infrared (775 nm). So, the Normalized Difference Vegetation Index (NDVI) and the Normalized Difference Red Edge (NDRE) were finally obtained and georeferenced by also acquiring the position using a Trimble AgGPS332 GPS with SBAS differential correction (EGNOS system).

An average of 169 data per row was obtained. Then, following the same procedure as before, the NDVI and NDRE data were interpolated by ordinary kriging to obtain two surface maps with a final resolution (grid) of 1 m for each vegetation index. By superimposing the surface maps with the layer of 104 polygons, the average values of NDVI and NDRE were obtained and denoted as variables E and D, respectively (Table 1). The use of a MTLS allowed four different ancillary variables to be obtained. The sensor used was a UTM30-LX-EW time-of-flight LiDAR (HOKUYO, Osaka, Japan) that could perform 40 scans per second (40 Hz). In addition to this 2D LiDAR sensor, the MTLS integrated a GPS1200+ (Leica Geosystems AG, Heerbrugg, Switzerland) RTK-GNSS system (a real-time kinematics global navigation satellite system receiving GPS and Glonass constellation signals).

So, after processing all data, the MTLS provided measurements of the position of the impact points produced between the laser beam and the canopy. These georeferenced impacts formed a point cloud that was then analyzed to calculate the canopy volume. A more detailed description of features, components for acquisition and field use methodology can be found in [START_REF] Escolà | Mobile terrestrial laser scanner applications in precision fruticulture/horticulture and tools to extract information from canopy point clouds[END_REF].

Like the OptRx sensor, the MTLS laterally scanned the canopies of all the rows within the plot, but in this case, each row was measured from both sides to obtain a point cloud representative of the entire canopy. Point cloud visualization and specific computations were performed with CloudCompare (CloudCompare [GPL software] v2.6.1 2015). In short, the 3D point cloud was first classified into canopy and ground points. Points located at a height less than 0.4 m above the ground were discarded to eliminate weeds and ridges. The remaining points were then processed to individualize each tree trying to eliminate possible impacts (outliers) that were clearly outside the canopy (presenting positions more than 2.5 standard deviations away from 6 neighboring points).

POSTPRINT of the article: Uribeetxebarria, A., Martínez-Casasnovas, J.A., Tisseyre, B., Guillaume, S., Escolà, A., Rosell-Polo, J.R., Arnó, J. 2019. Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards. Comput. Electron. Agric. DOI: https://doi.org/10.1016/j.compag.2019.104931 21

Once the cloud of impact points for each tree was delimited (in our case, 104-point clouds), the first variable that was obtained was the number of canopy impacts (variable L, Table 1).

CloudCompare 'octree' allowed the cubical initial volume including the entire canopy for each tree to be recursively divided into smaller cubes until visually adjusting the actual volume of the canopy. Adding the unit volume of each adjusted cube or voxel (VOlumetric piXEL), the final volume of the canopy in m 3 was obtained. Depending on either how the voxel size was set, manually or automatically, two different volumes could be obtained. The variable M (Table 1) referred to manually adjusted voxels, while the variable N (Table 1) estimated volume by adjusting a voxel size given by default parameters of the software.

The last volume (variable O, Table 1) was calculated using the '2.5D Volume' tool in the CloudCompare software. Having projected the point cloud for a specific tree on the XY plane (Z-axis projection direction), this tool allowed the volume between the 2D rasterized cloud and the ground surface (taken as arbitrary plane) to be computed. As multiple points of the canopy could fall inside each cell, the maximum height was taken for calculation.

Obviously, other methods could have been used to calculate the volume from LiDAR point clouds [START_REF] Escolà | Mobile terrestrial laser scanner applications in precision fruticulture/horticulture and tools to extract information from canopy point clouds[END_REF]). However, the special structure of the 'Catalan' vessel training system, open with branches in different directions and large gaps inside the canopy, could make it more convenient to use the voxelization procedure as proposed by [START_REF] Underwood | Mapping almond orchard canopy volume, flowers, fruit and yield using lidar and vision sensors[END_REF] in almond. The readings with MTLS system were made on May 17th, 2016.

POSTPRINT of the article: Uribeetxebarria, A., Martínez-Casasnovas, J.A., Tisseyre, B., Guillaume, S., Escolà, A., Rosell-Polo, J.R., Arnó, J. 2019

Evaluating the efficiency of sampling methods

To compare the SRS versus the RSS, B (B = 1000) realizations of selecting a sample with size N, with N ranging from 4 to 12 were resampled each time and for each method, with the sampling population being the 104 trees systematically sampled within the plot. By calculating the sample mean (average number of fruits per tree) for realization i, 𝑌 < , the variance of the distribution of means for the SRS scheme was computed according to Nane and Kooijman (2018) using expression (4),

𝑣𝑎𝑟"𝑌 DCD ' = ? W ∑ "𝑌 < -𝑌 . ' @ × (1 -𝑓) W <F? (4)
where B = 1000 was the number of realizations in the resampling method, 𝑌 < was the sample mean of the fruit load, and 𝑌 . was the average fruit load for the 104 trees previously sampled covering the entire area of the plot (𝑌 . = 144.7 fruits per tree). The latter was considered the representative value of the plot (or best estimator of the population mean), and it was practically coincident in all the resampling processes with the average of the 1000 means (data not shown).

The finite population correction factor was given by 1 -𝑓, with 𝑓 = 𝑁 104 ⁄ defining the sampling fraction.

POSTPRINT of the article: Uribeetxebarria, A., Martínez-Casasnovas, J.A., Tisseyre, B., Guillaume, S., Escolà, A., Rosell-Polo, J.R., Arnó, J. 2019 The variance of the mean (𝑌 CDD ) of a ranked set sample from a finite population was obtained by the formula (5) [START_REF] Kowalczyk | Ranked set sampling and its applications in finite population studies[END_REF],

𝑣𝑎𝑟"𝑌 CDD ' = ? B Z[1 - ? ?\H ] 𝑆 @ - ? B ∑ "_𝐸"𝑌 [𝑗] '`-𝑌 . ' @ B %F? a (5)
where N and 𝑌 . had the same meaning as before, and 𝑆 @ was the known variance of the total trees (104) sampled within the plot. To estimate each of the expected values, 𝐸"𝑌 [%] ', the mean over B = 1000 bootstrap realizations for each of the corresponding ranked yield values j (𝑗 = 1, … , 𝑁) was applied.

Both sampling methods (SRS and RSS) provided an unbiased estimate of the population mean (in this case, of the average of the 104 trees used as a representative measure of the plot mean).

So, in both cases, the mean squared error (MSE) equals the variance. The inverse of the variance is a measure of the precision, thus for a given sample size N the two sampling methods were compared using the so-called relative efficiency, RE (6),

𝑅𝐸 = cdd<e<fgeh "+ jkk ' cdd<e<fgeh "+ kjk ' = lmn"+ kjk ' lmn"+ jkk ' (6)
once obtained 𝑣𝑎𝑟"𝑌 DCD ' and 𝑣𝑎𝑟"𝑌 CDD ' by ( 4) and ( 5). RE values greater than 1 would be expected [START_REF] Chen | Ranked Set Sampling. Theory and Applications[END_REF], showing that RSS would be more precise (or more efficient) than SRS. At worst [START_REF] Webster | Field sampling for environmental science and management[END_REF], 𝑅𝐸 = 1 when the ranking process is completely imperfect, that is, there is no correlation between the fruit load and the ancillary variable used to rank the realizations. Also, as [START_REF] Wolfe | Ranked set sampling[END_REF] for each method was obtained using equation ( 7 × 100, being 𝑌 < the mean of the realization and 𝑌 . the aforementioned reference mean for the plot, the last indicator allowed the number of realizations of a given sampling design and sample size N with sampling error greater than 10%, 𝑛_𝜀"𝑌 < ' > 10%`, to be computed in percentage terms over the total of the 1000 realizations (8),

𝑃𝑅 =

g_v"+ , 'w?\%Ẁ × 100.

Note that this procedure was also used to test the optimal ancillary information to run the RSS.

For this, the resampling of B = 1000 realizations of size N = 12 using each auxiliary variable shown in Table 1 was used. Indices defined by equations ( 7) and ( 8) were used to compare the different auxiliary variables used in the ranking mechanism to provide an estimate of the fruit load. In all cases, resampling was performed by programming in R software, version 3.3.2.

POSTPRINT of the article: Uribeetxebarria, A., Martínez-Casasnovas, J.A., Tisseyre, B., Guillaume, S., Escolà, A., Rosell-Polo, J.R., Arnó, J. 2019. Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards. Comput. Electron. Agric. DOI: https://doi.org/10.1016/j.compag.2019.104931 25

Results and discussion

Relationship between fruit load and the acquired ancillary variables

Figure 3 shows the scatter plot defined by the first two factors of a principal component analysis (PCA) showing the correlation between the ancillary variables altogether and fruit load. By the two factors, 65.9% of variability is explained. In summary, many of the ancillary variables were well summarized by the two first components of the PCA allowing two sets of variables to be discriminated.

The variables with the highest correlation to the first component, apart from the fruit load (variable A), were those obtained using remote sensing (UAV and airplane using RGB and multispectral cameras) and MTLS. So, the first component was basically related to the geometry and vigor of the trees. The soil apparent electrical conductivity (shallow and deep values) was correlated to component 2 instead of component 1. In short, component 2 should be interpreted as related to soil variability. Finally, NDVI and NDRE vegetation indices obtained from the OptRx proximal sensor were not clearly included in any of the previous groups (therefore, they showed a weak correlation with both components). As has been mentioned before, "Catalan" vessel is the adopted training system which favours open canopies to improve solar exposure in detriment of the lateral side. In short, a vegetation wall is not easily measurable from the lateral side of the canopy. Because of that, the field of view of the canopy varies between sensors (lateral from terrestrial platform and nadir from aerial platform), and NDVI readings are surely affected. Probably, the best correlation with nadir vision is due to better cover the full size of the canopy. Relevant information about the relationship between fruit load and all other variables was contrasted visually by the scatterplot of a PCA (Fig. 3). The surprising weak correlation between the number of fruits (variable A) and component 2 made it unwise to use the Veris sensor and the ECa as ancillary information in RSS schemes. This weak correlation was not expected when systematic sampling was designed, according to the findings of other authors in fruit growing [START_REF] Käthner | Interaction of 3D soil electrical conductivity and generative growth in Prunus domestica L[END_REF]. In spite of this, variability at plot scale was also captured by covering the whole area of the orchard. The weak correlation could be due to drip fertigation allowing roots to grow in a controlled environment with little final soil influence (taking ECa readings in the alleyways outside the wet bulb) on canopy and fruit load [START_REF] Uribeetxebarria | Spatial variability in orchards after land transformation: Consequences for precision agriculture practices[END_REF]. There being little correlation between ECa and fruit load, the ranking process may not be optimal to obtain a sample with the additional information needed for a more efficient estimation of the mean. Figure 3 in RSS was rejected. However, the variables correlated to component 1 could be of interest. In short, the options for ancillary data to be used were limited to three data sources, i) RGB cameras mounted on remote acquisition platforms (UAV and airplane, variables F, G, J and K),

ii) airplane images supplied by DMSC (variables H and I), and iii) MTLS point clouds (variables L, M, N and O). The close spatial arrangement within component 1 (Fig. 3), including fruit load, helped in this interpretation.

Going into detail (Fig. 3), variable H (airplane-based NDVI) was the ancillary information with apparently the weakest correlation to fruit load among the variables of component 1. While variable H was only based on reflectance, the other variables were mainly based on geometric parameters. Even so, all the variables that were grouped as mostly correlated with component 1 were selected in order to analyze their convenience as ancillary variables in the next section.

With respect to the volume of the canopy, only the data calculated with the tool 'Volume 2.5D' (variable O) were considered to avoid redundant information

Ancillary information to be used in RSS to improve fruit load estimates

Table 2 shows the results of linear correlation between the number of fruits per tree (variable A) and the most relevant ancillary variables within component 1. As expected, the best correlation was obtained using the UAV-based canopy projected area (variable G) with a value of the linear correlation coefficient of 0.85. The worst result (coefficient of 0.21) corresponded to the airplane NDVI (variable H). The rest of the variables ranging between these two extreme values (Table 2).

POSTPRINT of the article: Uribeetxebarria, A., Martínez-Casasnovas, J.A., Tisseyre, B., Guillaume, S., Escolà, A., Rosell-Polo, J.R., Arnó, J. 2019 As expected, Table 2 shows that MSE values follow exactly the same order as the correlation coefficient showing that the higher the correlation, the lower the error in fruit load prediction with the RSS. The trend of the Coefficient of Error (𝐶𝐸) was similar, going from 5.74% when the ancillary variable was the UAV-based canopy projected area to 8.93% using the airplane NDVI to rank items within the samples. Since a CE of 10% may be acceptable in the fruit production sector, any of the ancillary variables could be candidates to be used in RSS schemes.

However, it is also true that the sample size used in the comparison (N = 12) is not usual in practice, and has undoubtedly contributed to improve the efficiency in fruit load estimates above what is expected. Hence, it is necessary to check how the efficiency varies by reducing the sample size closer to the usual sizes used by fruit growers.

Table 2 also shows the percentage of times the sampling error exceeded the 10% threshold considering 1000 realizations with a sample size N = 12. This probability is a way of knowing POSTPRINT of the article: Uribeetxebarria, A., Martínez-Casasnovas, J.A., Tisseyre, B., Guillaume, S., Escolà, A., Rosell-Polo, J.R., Arnó, J. 2019. Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards. Comput. Electron. Agric. DOI: https://doi.org/10.1016/j.compag.2019.104931 29 the risk assumed by fruit growers since they usually perform a single sampling. Values ranged from 10.1% (variable G) to 26.9% (variable H), and corroborate the results of MSE and linear correlation. The maximum difference between ancillary variables was 166% for this parameter (PR), being 55% for the CE.

On the other hand, the results shown in Table 2 were in accordance with what was noticed by [START_REF] Nahhas | Ranked set sampling: cost and optimal set size[END_REF] and [START_REF] Stokes | Ranked set sampling with concomitant variables[END_REF] when applying RSS. The better the correlation between the variable to be estimated (number of fruits per tree, in our case) and the concomitant variable (or ancillary variable) used in the ranking process, the better the sampling efficiency.

Apart from good efficiency results, the ancillary variable to be used must be easily measurable.

RGB cameras embedded in UAV meet this requirement. The use of UAV in precision agriculture is a booming and low-cost technology as suggested by [START_REF] Lelong | Assessment of unmanned aerial vehicles imagery for quantitative monitoring of wheat crop in small plots[END_REF]. For this reason, the UAV-based canopy projected area (variable G, Tables 1 and2) was finally chosen as the optimal ancillary variable to use in RSS schemes in fruit growing. The advantage of this variable is the high resolution of the images compared to other remote data sources.

Resolution of ancillary variables in RSS: the key factor

In peach orchards, fruit production and quality are clearly influenced by canopy lighting conditions [START_REF] Tang | Light interception efficiency analysis based on three-dimensional peach canopy models[END_REF][START_REF] Minas | Environmental and orchard bases of peach fruit quality[END_REF]. Seeking to enhance floral induction and fruit growth, fruit growers adopt canopy open-center training systems (such as the typical 'Catalan' peach vessel) to maximize exposure to solar radiation. In this way, peach trees with larger canopy projected area usually have a greater number of larger fruits [START_REF] Marini | Heading fruiting shoots before bloom is equally effective as blossom removal in peach crop load management[END_REF]. In addition, POSTPRINT of the article: Uribeetxebarria, A., Martínez-Casasnovas, J.A., Tisseyre, B., Guillaume, S., Escolà, A., Rosell-Polo, J.R., Arnó, J. 2019. Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards. Comput. Electron. Agric. DOI: https://doi.org/10.1016/j.compag.2019.104931 30 the incidence of unproductive interior branches is lower, avoiding what usually happens in canopies formed into traditional or spherical vessels [START_REF] Marini | Tree form and heading height at planting affect peach tree yield and crop value[END_REF]. In RSS, an ancillary information system that manages to individualize the canopy of trees in a precise manner will be essential, i) to ensure an optimal and accurate ranking process of the auxiliary variable, and ii) to cover all the variability of the objective variable (number of fruits) (Fig 2).

Continuing the thread of the previous paragraphs, variables that best meet previously mentioned characteristics are those based on geometrical properties of the canopy (Table 2). UAV-based canopy projected area (variable G) was the best rated ancillary variable. In contrast, LiDARbased canopy volume (O) or impacts (L) and airplane-based canopy projected area (K) obtained less satisfactory sampling efficiency results. What could be the reason for this difference?

Probably, the lack of resolution (57528 impacts per tree) is not the main reason for a lower correlation of LiDAR derived variables with fruit load. The lowest correlation comes from the different orientation of each sensor. While the RGB camera mounted on UAV (variable G) provides a top view of the canopy, MTLS provides a lateral view of the trees. The top view allows the area exposed to sunshine to be more precisely delimited. The efficiency when intercepting the sunrays is the key factor for transpiration and photosynthesis, therefore it is strongly related to the fruit load [START_REF] Da Silva | Light interception efficiency of apple trees: a multiscale computational study based on MappleT[END_REF].

Figure 4A shows the canopy of a peach tree as it is captured with an aerial image from airplane (0.25 m per pixel), and the same canopy when it was captured using UAV and resolution of 0.02 m (Fig. 4B). The difference is obvious. While the delimitation of the canopy projected area was more difficult using the airplane image, the UAV-based canopy projected area allowed the actual canopy projected area to be better defined. This would explain why the use of higher resolution images should lead to more accurate tree ranking processes and more efficient fruit load estimates. To provide an explanatory example, 45 pixels were needed on average to characterize the canopy projected area using the airplane image instead of the 6500 pixels for the UAV image. In other words, the average area of 2.60 m 2 per tree using the higher resolution image became an area of 2.81 m 2 when a larger pixel size was used. The difference was therefore 8%, although large deviations for some trees could be obtained (maximum differences of up to 53% between UAV and airplane were found). Because of that, the use of images of poor resolution can alter the ranking process. In fact, ranking of individuals varied (data not shown) according to the resolution of the ancillary variable that was used (for example, G or K). The moderate Spearman correlation coefficient (0.71) between both variables allowed this property to be concluded. POSTPRINT of the article: Uribeetxebarria, A., Martínez-Casasnovas, J.A., Tisseyre, B., Guillaume, S., Escolà, A., Rosell-Polo, J.R., Arnó, J. 2019. Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards. Comput. Electron. Agric. DOI: https://doi.org/10.1016/j.compag.2019.104931 32 MTLS are another different technology that allows canopy of trees to be characterized with high operative resolution. Figure 4C shows the point cloud obtained from a lateral LiDARbased scan of the same tree. Then, by 3D processing of all the data, a canopy volume expression can be obtained (variable O, Table 2). However, this parameter did not provide results as good as those obtained using the simplest UAV-based tree canopy projected area (variable G). Finding an explanation is not easy. Probably, sunlight penetration is optimal in trees with more open canopy resulting in a larger projected area. Conversely, bulky trees with smaller projected areas could have grown in height instead of laterally and more openly. This could hinder the interior lighting of the canopy and negatively affect the amount of viable fruits and yield. In fact, the correlation between canopy projected area (variable G) and fruit load is higher (R = 0.85) than the correlation between canopy volume (variable O) and fruit load (R = 0.68). POSTPRINT of the article: Uribeetxebarria, A., Martínez-Casasnovas, J.A., Tisseyre, B., Guillaume, S., Escolà, A., Rosell-Polo, J.R., Arnó, J. 2019. Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards. Comput. Electron. Agric. DOI: https://doi.org/10.1016/j.compag.2019.104931

33

The use of NDVI showed acceptable results when adjusting the value of the vegetation index according to the tree canopy projected area. However, the exclusive use of the airplane NDVI (variable H, Table 2) was the least recommended option. The resolution of the airplane images was lower than that of UAV images or MTLS point clouds and, in addition, the reflectance alone does not account for any geometrical parameters of the canopy. In any case, the weak correlation between fruit load and NDVI is also noted in other sampling studies [START_REF] Arnó | Comparing efficiency of different sampling schemes to estimate yield and quality parameters in fruit orchards[END_REF][START_REF] Miranda | Sampling stratification using aerial imagery to estimate fruit load in peach tree orchards[END_REF].

Concluding, resolution of acquired ancillary data was the main constraint in reliably delimiting canopy of trees (Fig. 4). UAV-based RGB images offered high resolution at competitive costs, and good correlation with the fruit load (Table 2). Therefore, RSS was only tested with UAVbased tree canopy projected area as the ancillary information. Comparison of sampling efficiency against SRS is addressed in the next section. The pending issue is to test the method in other orchards trained under systems that are more intensive forming a fruiting wall structure.

Probably in this case, sections of wall vegetation rather than individualized trees should be delimited using an appropriate algorithm and a more automated method.

RSS versus SRS: sampling efficiency and sample size

Table 3 shows the mean squared error (MSE), and the relative efficiency (RE) of RSS versus SRS for different sample sizes (4 to 12 points or trees sampled within the plot). A resampling of 1000 realizations allowed such statistics to be calculated in order to compare the sampling methods. RSS was run in all cases using the UAV-based tree canopy projected area as ancillary variable to drive the ranking mechanism.

POSTPRINT of the article: Uribeetxebarria, A., Martínez-Casasnovas, J.A., Tisseyre, B., Guillaume, S., Escolà, A., Rosell-Polo, J.R., Arnó, J. 2019 In both cases, SRS and RSS, the MSE decreased, as the sample size increased. As expected, the sampling methods became more efficient with an increasing sample size. However, comparatively, RSS always showed lower variance of the mean (lower MSE) than SRS. In short, RSS was more efficient than SRS for any of the sample sizes evaluated. Hence the RE was higher than one (Table 3), ranging from 1.55 (N = 4) to 2.15 (N = 12). The fact of obtaining higher values of RE for the larger sample sizes was possibly due to the RSS procedure.

Takahashi and Wakimoto (1968) reported a loss of efficiency by decreasing the number of individuals checked in the ranking process for the same final sample size. This would explain why the RE of the RSS increases with increasing sample size because a higher number of trees are ranked. Thus, by increasing the sample size (N = 7, 8 ... 12), RSS required progressively increasing the number of trees to rank to obtain the final sample. For example, a sample of size N = 9 would suppose choosing 81 trees within the plot instead of the 25 (5 2 ) needed for N = 5.

In this way, additional information on fruit load structure was more easily achievable by having

Fig 1

 1 Fig 1 Location of the study area and orthophoto (2016) showing the 104 trees sampled within

Fig 2

 2 Fig 2 Distribution of the fruit load (histogram of the 104 trees sampled), normal density curve

  Ltd, Leederville, Western Australia) to correct lens aberration, and adjust scene brightness by the bidirectional reflectance distribution function (BRDF). At the moment of photographing the plot, the airplane was flying at an altitude of 2000 m. The spatial resolution of the pre-processed image was 25 cm per pixel. The aircraft used was a CESSNA 1725 Sky Hawk operated by RS (RS Servicios de Teledetección SL, Lleida, Spain).

  2.4.3. Tree canopy reflectanceThe aerial image provided by the multispectral camera(NDVI image, Rouse Jr et al. 1974) was post-processed according to the following procedure. First, and having eliminated the NDVI pixels with values lower than 0.45 to discard other land covers (i.e. bare soil) that were not tree canopy, all the trees (or canopies) within the plot were then easily delimited by using automatically defined polygons. Specifically, 1816 trees were individualized this way. Then, in a second step, NDVI pixels within each polygon (or canopy) were extracted by overlapping the NDVI image and the polygon layer. The average value of NDVI was then calculated and assigned to the centroid of each polygon. The third step consisted in creating a continuous map (or surface map) by interpolation (ordinary kriging) of the average NDVI values of the centroids. By superimposing this new NDVI map (1 m pixel resolution) with the polygons corresponding to the 104 trees, an average value of NDVI for each tree was obtained. This

  p 𝑣𝑎𝑟(𝑌), the standard deviation of the mean. Having formulated the sampling error of one realization as 𝜀"𝑌 <

Fig 3

 3 Fig 3 Scatterplot showing the correlation between the ancillary variables and fruit load (variable

Fig 4

 4 Fig 4 Comparison of tree canopy projected area delimited for the same peach tree, A) using an

Figure 5

 5 Figure 5 helps to understand this phenomenon. While the two trees shown (number 104 and

Fig 5

 5 Fig 5 Aerial projection of tree 104 (left) and 45 (right) obtained from the LiDAR-based MTLS

  Uribeetxebarria, A., Martínez-Casasnovas, J.A., Tisseyre, B., Guillaume, S., Escolà, A., Rosell-Polo, J.R., Arnó, J. 2019. Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards. Comput. Electron. Agric. DOI: https://doi.org/10.1016/j.compag.2019.104931 14 Therefore, the RSS mean is, theoretically, a more precise estimator of the population mean 𝜇 than the SRS mean based on the same number of measured observations. The reliability of the judgment ranking process in separating order statistic expectations 𝜇[%]

Table 1

 1 Uribeetxebarria, A., Martínez-Casasnovas, J.A., Tisseyre, B., Guillaume, S., Escolà, A., Rosell-Polo, J.R., Arnó, J. 2019. Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards. Comput. Electron. Agric. DOI: https://doi.org/10.1016/j.compag.2019.104931 15 Sensors and ancillary information that were evaluated for the ranking of sampled trees in the RSS method.POSTPRINT of the article: Uribeetxebarria, A., Martínez-Casasnovas, J.A., Tisseyre, B., Guillaume, S., Escolà, A., Rosell-Polo, J.R., Arnó, J. 2019. Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards. Comput. Electron. Agric. DOI: https://doi.org/10.1016/j.compag.2019.104931 POSTPRINT of the article:Uribeetxebarria, A., Martínez-Casasnovas, J.A., Tisseyre, B., Guillaume, S., Escolà, A., Rosell-Polo, J.R., Arnó, J. 2019. Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards. Comput. Electron. Agric. DOI: https://doi.org/10.1016/j.compag.2019.104931 17 GX 7 camera (Panasonic Corporation, Osaka, Japan) with a 20 mm "pancake" lens coupled on a Mikrokopter Oktokopter 6S12 XL eight rotor UAV(HiSystems GmbH, Moomerland,

	Acquisition	Sensor	View of the	Ancillary	Geometry-	Reflectance-	Obtaining
	platform		canopy	information	based	based	process
	ATV	Veris 3100	Soil sensor	(B) shallow ECa			ME

Two ECa raster maps (shallow and deep) were then obtained by ordinary kriging interpolation on a 1 m grid. By superimposing the layer of 104 polygons corresponding to the 104 trees within the plot, the respective values of shallow ECa (variable B, Table

1

) and deep ECa (variable C, Table

1

) were extracted. The mean ECa value was calculated for each tree at each surveying depth.

2.4.2. Tree canopy projected area and tree canopy perimeter

Tree canopy projected area and tree canopy perimeter were obtained by processing two images acquired from two different platforms: airplane and unmanned aerial vehicle (UAV). Both images were acquired on May 16, 2016, approximately at mid-day and under clear sky conditions.

Each platform had its own camera. For the UAV, images were grabbed with a 16 MP Panasonic

  POSTPRINT of the article:Uribeetxebarria, A., Martínez-Casasnovas, J.A., Tisseyre, B., Guillaume, S., Escolà, A., Rosell-Polo, J.R., Arnó, J. 2019. Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards. Comput. Electron. Agric. DOI: https://doi.org/10.1016/j.compag.2019.104931 18 (K), UAV-based canopy perimeter (F), and UAV-based canopy projected area (G).

Table

1

shows the four variables (symbolized with a capital letter) that were used as ancillary information on the canopy's geometry: airplane-based canopy perimeter (J), airplane-based canopy projected area

  . However, average NDVI could take similar values in trees of different sizes. With this in mind, formula (3) was used to correct values of the NDVI to be obtained for each tree (variable I, Table1), POSTPRINT of the article: Uribeetxebarria, A., Martínez-Casasnovas, J.A., Tisseyre, B., Guillaume, S.,Escolà, A., Rosell-Polo, J.R., Arnó, J. 2019. Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards. Comput. Electron. Agric. DOI: https://doi.org/10.1016/j.compag.2019.104931

  POSTPRINT of the article:Uribeetxebarria, A., Martínez-Casasnovas, J.A., Tisseyre, B., Guillaume, S., Escolà, A., Rosell-Polo, J.R., Arnó, J. 2019. Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards. Comput. Electron. Agric. DOI: https://doi.org/10.1016/j.compag.2019.104931

	2.4.4. Tree volume
	The date of the terrestrial sensor
	measurements was May 17th, 2016.
	20

  . Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards. Comput. Electron. Agric. DOI: https://doi.org/10.1016/j.compag.2019.104931 22

  . Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards. Comput. Electron. Agric. DOI: https://doi.org/10.1016/j.compag.2019.104931 23

  points out, good concomitant POSTPRINT of the article: Uribeetxebarria, A., Martínez-Casasnovas, J.A., Tisseyre, B., Guillaume, S., Escolà, A., Rosell-Polo, J.R., Arnó, J. 2019. Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards. Comput. Electron. Agric. DOI: https://doi.org/10.1016/j.compag.2019.104931 24 information is necessary to avoid error-prone ranking.

	Two additional efficiency indicators were computed. First, a predicted Coefficient of Error (CE)

  also shows a nearly zero relationship between fruit load and variables D and E obtained by the OptRx sensor. For this reason, its use POSTPRINT of the article: Uribeetxebarria, A., Martínez-Casasnovas, J.A., Tisseyre, B., Guillaume, S., Escolà, A., Rosell-Polo, J.R.,Arnó, J. 2019. Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards. Comput. Electron. Agric. DOI: https://doi.org/10.1016/j.compag.2019.104931 27

Table 2

 2 . Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards. Comput. Electron. Agric. DOI: https://doi.org/10.1016/j.compag.2019.104931 28 Assessment of RSS efficiency according to the ancillary variables used in the ranking mechanism. Variables are ordered according to their correlation with fruit load.

	R	𝑴𝑺𝑬"𝒀' 𝑪𝑬"𝒀'	PR

Table 3

 3 . Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards. Comput. Electron. Agric. DOI: https://doi.org/10.1016/j.compag.2019.104931 34 Mean squared error (MSE) and relative efficiency (RE) of the sampling methods for different sample sizes.

	Mean squared error (MSE)

(C) deep ECa ME ATV OptRx Lateral view (D) NDRE X ME (E) NDVI X ME UAV RGB camera Nadir (F) Tree canopy perimeter X ME (G) Tree canopy projected area X ME Airplane DMSC Nadir (H) NDVI X ME (I) NDVIC X X HG RGB camera (J) Tree canopy perimeter (K) Tree canopy projected area

Acknowledgments:

The Spanish Ministry of Economy and Competitiveness through the AgVANCE Project (AGL2013-48297-C2-2-R) funded this work. We are also grateful to Frutas Hermanos Espax SL for the possibility to carry out the research in their farm.

POSTPRINT of the article: Uribeetxebarria, A., Martínez-Casasnovas, J.A., Tisseyre, B., Guillaume, S., Escolà, A., Rosell-Polo, J.R., Arnó, J. 2019. Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards. Comput. Electron. Agric. DOI: https://doi.org/10. 1016/j.compag.2019.104931

to select and assess a larger number of trees than the one obtained in small samples. As a result, the RE of RSS increased compared to SRS strategy.

Considering the UAV-based canopy projected area as ancillary information, RSS allowed more accurate fruit load estimates to be obtained. The extent to which this precision was improved was quantified by the CE (Eq. 7) as shown graphically in Fig. 6. Additionally, Figure 6 shows the effect of the sample size on the PR indicator for both sampling schemes. Both measures are important for fruit growers because they complement each other. Farmers collect samples at specific times during the life of the crop. Therefore, knowing the probability of exceeding the threshold and the average CE for each sample size is important. Both curves (CE and PR, Fig. 6) show a similar nonlinear trend. For all sample sizes, RSS has a lower CE than SRS. Specifically, while in SRS the CE ranged from 15.4% (N = 4) to 8.4% (N = 12), in RSS errors ranged from 12.1% (N = 4) to 5.7% (N = 12). From a practical point of view, SRS should be based on minimum sizes of N = 10 while RSS can reduce sample sizes up to N = 5 so as not to exceed a CE threshold of 10% (Fig. 6). This meant that, using sample sizes N = 10 in SRS or N = 5 in RSS, a sampling error of less than 10% could be expected in practically 70% of the samplings (PR close to 30% in both cases, Fig. 6). Assuming normality, 68.3% of the samplings should provide errors of less than 10% (as it has been obtained approximately), although the RSS strategy only needed to sample 5 trees compared to the 10 trees needed in SRS. Likewise, SRS schemes exceed the error threshold of 10% with a higher percentage (PR, Fig. 6) compared to RSS for all the sample sizes that were analyzed. While for SRS, PR values ranged from 53.5% (N = 4) to 25.5% (N = 12), in RSS the same parameter reached significantly lower scores, varying between 44.4% (N = 4) and 10.1% (N = 12). Therefore, using SRS practically doubles (for sample sizes greater than N = 8) the probability of obtaining errors of estimation above the POSTPRINT of the article: Uribeetxebarria, A., Martínez-Casasnovas, J.A., Tisseyre, B., Guillaume, S., Escolà, A., Rosell-Polo, J.R., Arnó, J. 2019. Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards. Comput. Electron. Agric. DOI: https://doi.org/10.1016/j.compag.2019.104931 36 10% threshold than using RSS. However, for smaller sample sizes (N = 4), the difference between SRS and RSS was not so obvious. SRS (simple random sampling). RSS (ranked set sampling). Resampling of 1000 realizations for each of the sample sizes was made, using the UAV-based tree canopy projected area as ancillary variable in the RSS.

As mentioned before, if a CE of no more than 10% was assumed, RSS allowed the sample size to be reduced to N = 5 while the SRS required twice the sample (N = 10) (Fig. 6). Moreover, a sample size N = 4 in RSS could even be justified by reaching an error around 12%.

Undoubtedly, being able to reduce the sampling requirements, RSS could report favorable operative as well as economic consequences in fruit sampling as those obtained by [START_REF] Carrillo | Use of multi-spectral airborne imagery to improve yield sampling in viticulture[END_REF] in vineyard. Although RSS requires random sampling in each iteration to generate the final ranked sample, subjectivity (and, therefore, biased decisions) should not be applied during the ranking mechanism favoring a more automated (without farmers rating) and reliable sampling process. On the other hand, the ancillary variable being available, the method is even POSTPRINT of the article: Uribeetxebarria, A., Martínez-Casasnovas, J.A., Tisseyre, B., Guillaume, S., Escolà, A., Rosell-Polo, J.R., Arnó, J. 2019. Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards. Comput. Electron. Agric. DOI: https://doi.org/10.1016/j.compag.2019.104931 37 applicable in small plots like fruit orchards [START_REF] Zude-Sasse | Applications of precision agriculture in horticultural crops[END_REF], avoiding to perform a stratified sampling scheme.

Looking for the reasons that make the RSS a strategy more adapted to non-homogeneous orchards, it was verified (data not shown) that, in average, the within-sample variability (measured as Coefficient of Variation) in the RSS realizations was always higher than that found in SRS of the same size. This would explain why the RSS manages to capture the heterogeneity of the orchards more reliably providing greater precision in yield estimates.

Finally, it is necessary to point out that these values obtained for this orchard assume that the sampling variance on estimating mean fruit per tree due to a (in effect) selection of 104 trees as the first step can be considered negligible. In addition, these recommended sample sizes apply to this particular orchard with its spatial variability of fruit. They provide a starting point for considering sample sizes, and being prepared to be more conservative for orchards that are more heterogeneous.

Conclusions

Ranked set sampling (RSS) is an interesting method to improve sampling in fruit growing.

Improvements in the efficiency of peach fruit load estimates were proven by applying RSS compared to simple random sampling (SRS), at least for sample sizes from N = 4 to N = 12. However, to be optimal, RSS requires a relevant auxiliary variable, these data being the key factor in the procedure. Among different data sources currently available in agriculture, either commercialized or under development (such as MTLS), the UAV-based tree canopy projected POSTPRINT of the article: Uribeetxebarria, A., Martínez-Casasnovas, J.A., Tisseyre, B., Guillaume, S., Escolà, A., Rosell-Polo, J.R., Arnó, J. 2019. Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards. Comput. Electron. Agric. DOI: https://doi.org/10.1016/j.compag.2019.104931 38 area was identified in this study of a peach orchard as the best option to use in the ranking mechanism aiming at estimate fruit number in peach production. This information can be obtained from RGB cameras that providing high-resolution images allowing the projected area of the canopy (geometric parameter) to be accurately delimited. The close correlation between this area and the fruit load per tree made this the best choice as an auxiliary variable in RSS schemes for the particular management by the producer (open vessel training systems).

In terms of the Coefficient of Error of the estimated mean, RSS using samples of size N = 5 trees allows acceptable errors below 10% to be achieved with about 70% probability. In contrast, SRS requires practically doubling the sample needed for a similar error. Finally, it should be noted that the work has been carried out only on one plot, although it provides relevant information as the ancillary variable to be used in RSS for future research on a larger spatial scale. In fact, UAV-based tree canopy projected area is an interesting (cheap and easy to get) ancillary information. However, since the values of the variable are obtained manually, future research is required to automatically obtain this information by developing the appropriate algorithm.