Estimation of Leaf Area Index in vineyards by analysing projected shadows using UAV imagery

A few decades ago, farmers could precisely monitor their croplands just by walking over the fields, but this task becomes more difficult as farm size increases. Precision viticulture can help better understand the vineyard and measure some key structural parameters, such as the Leaf Area Index (LAI). Remote Sensing is a typical approach to monitoring vegetation which measures the spectral information directly emitted and reflected from vegetation. This study explores a new method for estimating LAI which measures the projected shadows of plants using UAV (unmanned aerial vehicle) imagery. A flight mission over a vineyard was scheduled in the afternoon (15:30 to 16:00 solar time), which is the optimal time for the projection of vine shadows on the ground. Real LAI was measured destructively by removing all the vegetation from the area. Then, the projected shadows in the image were detected using machine learning methods (k-means and random forest) and analysed at pixel level using a customised R code. A strong linear relationship (R2 = 0.76, RMSE = 0.160 m2 m-2 and MAE = 0.139 m2 m-2) was found between the shaded area and the LAI per vine. This is a quick and simple method, which is non-destructive and gives accurate results; moreover, flights can be scheduled during other periods of the day than solar noon, such as in the morning or afternoon, thus enabling pilots to extend their working day. Therefore, it may be a viable option for determining LAI in vineyards trained on Vertical Shoot Positioned (VSP) systems.


INTRODUCTION
According to European Union policies (Zarco-Tejada et al., 2014), Precision Agriculture is a farming management concept based on observing, measuring, and responding to interand intra-field variability in crops, and in which the spatial variability of vineyards plays a vital role. Vineyards are characterised by a strong spatial structure that is generally stable over time (Bramley and Lamb, 2003) and is affected by several factors, some of which are quite constant over time and can be dealt with through differentiated crop management (Bramley and Hamilton, 2004). Although Precision Agriculture can improve yields, the most significant advantage is the reduction of yield variations over time, which leads to more stability and resilience to a changing climate (Yost et al., 2016). On the other hand, until a few decades ago, farmers could precisely monitor their croplands just by walking over the fields, but with increasing farm size this task is becoming more difficult without using technology (Balafoutis et al., 2017). All of this has led to a growing interest in precision viticulture and further research efforts (Santesteban, 2019), but results can vary significantly depending on knowledge of the crop. Therefore, it is necessary to use methods that precisely define key variables to avoid any problems; for example, over-cropping, which will lead to a canopy with an insufficient amount of healthy active leaves which will not be able to produce enough sugar to obtain the desirable level of ripening in all clusters, resulting in grapes lacking in aroma/flavour and/or the desirable phenol compounds (Reynolds, 2010). Key variables to be defined include cluster, flower or berry number, which directly affect yield and ripening, and leaf area, which is an important variable to monitor, because the ripening of grape clusters depends on the leaf area/fruit weight ratio (Keller, 2015;Jackson, 2020). Leaves are organs specialised in intercepting radiation, which is necessary for photosynthesis. Leaf area is commonly measured using the Leaf Area Index, LAI (Lambers and Oliveira, 2019), which is closely related to photosynthetic active radiation assimilation (Pessarakli, 2014) and is defined as the vegetative development of a crop per unit area of land (Watson, 1947).

Leaf Area Index in viticulture
In viticulture, LAI is one of the most used parameters to represent canopy area (Delrot et al., 2010); it is a key indicator, since leaf area directly correlates with other critical parameters, such as transpiration, root development and photosynthetic capacity, thus limiting yield (Keller, 2015). Furthermore, LAI is related to canopy transpiration and water use, and can thus influence irrigation decisions (Netzer et al., 2008;Munitz et al., 2019), and it is affected by management and environmental factors, such as irrigation, nutrient management, and training systems (Oliveira and Santos, 1995).

Current Remote Sensing methods for estimating leaf area index
LAI can be measured using traditional methods, such as the Carbonneau method (Carbonneau, 1976), the adapted Point Quadrat (Smart and Robinson, 1991), the Lopes and Pinto method (Lopes and Pinto, 2005), or by measuring specific parameters related to leaf shape (Williams and Martinson, 2003). However, these methods are inefficient, time-consuming and, in some cases, destructive.
Remote Sensing is a quicker way of estimating LAI via; different technologies can be used, such as field spectroradiometers , multispectral and hyperspectral data (Mananze et al., 2018), satellite imagery (Meyer et al., 2019;Dube et al., 2019) and thermal imagery (Neinavaz et al., 2019). These technologies can vary in accuracy and complexity:  i) Different results have been reported for satellite imagery; for example, Johnson et al. (2003) showed a significant correlation between image-based leaf area and ground-based leaf area (R² > 0.7), while Beeri et al. (2020) found a weak relationship between satellite information and LAI (R² > 0.3).
 ii) Specific tools have been developed exclusively for indirect LAI estimation based on the measurement of radiation extinction through the foliage, such as the Plant Canopy Analyzer (PCA, LAI-2000 and 2200, Li-Cor Inc., Lincoln, NE, USA). Johnson and Pierce (2004) have reported high accuracy using PCA in viticulture, but even though this tool has been specifically developed to measure LAI, it is not easy to employ it correctly, since there are different protocols for its use (White et al., 2018).
 iii) General-purpose sensors such as conventional RGB cameras have also shown great potential for measuring LAI using various methods. For example, Fuentes et al. (2014) used a camera mounted on a pole to acquire downward-looking digital images from the middle of the plant, and De Bei et al. (2016) developed an app (VitiCanopy) to assess canopy architecture parameters in vineyards using upward-looking imagery. Diago et al. (2012) used an RGB camera mounted on a tripod set normal to the canopy plane, 2 m from the row axis and 1.05 m aboveground. Both authors used an RGB camera successfully, reporting high accuracy. Furthermore, these RGB cameras can be mounted on unmanned aerial vehicles (UAVs), thus increasing their possibilities improving the accuracy of the results. Kalisperakis et al. (2015) used a lowcost standard RGB camera to estimate LAI, obtaining good results (R² > 0.7), and they improved the results (R² > 0.8) by mounting a hyperspectral VNIR imaging sensor on a UAV.  iv) Hyperspectral and multispectral sensors are an improvement on RGB cameras. Using multispectral airborne images, Hall et al. (2008) found close relationships between the planimetric canopy area and ground-based measurements of LAI at several phenological stages, whereas no significant relationships were found between NDVI and LAI. Towers et al. (2019) used a multispectral sensor to capture nadir-view images of the canopy, reporting varying accuracies depending on whether soil values were used or not, and showing that soil backscattering can contribute more to the signal than vegetation cover can. In a similar way, but capturing the images from a UAV, Comba et al. (2020) used a multispectral camera combined with point cloud creation and 3D modelling using SfM (Structure from Motion) to estimate LAI, obtaining good results (R² > 0.8). Mathews and Jensen (2013) captured images from nadir and varying oblique angles, obtaining less accurate results than Comba et al. (2020) (R² > 0.5), and showing that increasing the number of differing angles/perspectives with overlap improves the SfM product. Kalisperakis et al. (2015) reported good accuracy (R² > 0.8) by creating a 3D model via 3D triangulation; they used aerial RGB images to generate a dense point cloud by employing dense stereo and multiimage matching algorithms. Therefore, even when the same sensors and technology (RGB cameras) are used, different methods can result in different levels of accuracy.
 v) An alternative approach involves the use of specific tools that can capture the three dimensions (3-D) of the area. For example, Arnó et al. (2013) computed geometric and structural parameters using a tractor-mounted LiDAR system to measure vines (using TAI, tree area index) in a transverse direction along rows with high accuracy (R² > 0.8).
Each method has its advantages and limitations. Some are cheap to implement (e.g., RGB cameras (Fuentes et al., 2014;Diago et al., 2012)), but may have different drawbacks depending on the method; for example, Diago et al. (2012) had problems using RGB cameras in terms of distance to the remnant foliage when the defoliation process was performed over highly dense canopies. Del- Moral-Martínez et al. (2016) reported that LiDAR needs to be very precisely set up, otherwise it can generate incorrect LAI values for vines with poor leaf development, recording zones in the canopy containing a considerable percentage of gaps as effective leaf wall area. Mathews and Jensen (2013) reported that the SfM approach needs more time than other UAV missions due to the higher number of images, angles or overlap. In addition, some methods are more expensive than others due to the technologies involved, such as multispectral (Towers et al., 2019;Comba et al., 2020) and hyperspectral cameras (Kalisperakis et al., 2015) or a tractor for mounting the LiDAR system on. Even a low-cost version of LiDAR (Arnó et al., 2013) can still be much more expensive than a standard RGB camera.
For most of these methods, strong correlations with LAI can generally be found (R² > 0.7), and a key requirement of most of them is to measure data around midday to minimise the effects of shadowing between the vine rows. The development of new methods that will complement current techniques and help expand target areas and extend the work to other periods of the day could be very beneficial for the application of these techniques in practice.

Remote Sensing approach to shadows
In Precision Agriculture, Remote Sensing is typically used for monitoring vegetation by measuring the spectral information directly emitted and reflected from the vegetation; each band is analysed separately, or indices are calculated, such as NDVI (Rouse et al., 1973), which is related to vineyard vegetation (Vélez et al., 2020b). In these approaches, shadows are frequently treated as non-desirable information (Ma et al., 2008;Zhang and Chen, 2010;Wu and Bauer, 2013;Aboutalebi et al., 2018) and are usually removed from datasets (Poblete-Echeverría et al., 2017;Jiang et al., 2018). Nevertheless, some problems can arise when the images are not obtained close to solar noon, as a result of the shaded and sunlit parts of leaves having different reflectance values. In addition, shaded areas can generate noisy data, which significantly affect the estimation of vegetation parameters (Zhang et al., 2015).

Proposed approach
The leaf area of the plant is correlated with intercepted light and can be estimated by measuring its shadow (Baeza et al., 2010). Therefore, a newer and less time-consuming approach to estimating LAI could be the study of the shaded soil area within the vineyard using UAVs. Each plant projects its own shadow and image analysis can be used to measure the lateral leaf area of the vines in the shaded inter-row space.
UAV platforms have been extensively used for studying and exploring vineyards, offering valuable technology for estimating numerous vineyard parameters (e.g.; Mathews and Jensen, 2013;Zarco-Tejada et al., 2013;Matese et al., 2016;Santesteban et al., 2017;Weiss and Baret, 2017;Poblete-Echeverría et al., 2017). In general, UAV offers the possibility of obtaining precise and high-resolution multispectral imagery, which is critical for measuring the shape of the vine shadows correctly. Furthermore, it allows the time of the flight to be chosen (Vélez et al., 2020a), which was a crucial factor in this study as the flight needed to be scheduled according to sun elevation, period of the year and desired shadow size (related to plant height). Other authors have hypothesised that shadows can be used as an indirect measure of leaf area and canopy characteristics (Zheng and Moskal, 2009), or have even suggested the possibility of using optical remote Sensing close to solar noon to monitor and map vineyard shaded area (Johnson and Scholasch, 2005). However, to our knowledge, this is the first study to follow and evaluate a field protocol for measuring vine leaf area using the shaded area under real conditions. This study takes advantage of advances in image capturing and UAV technologies, combining previous knowledge in viticulture with remote Sensing and machine learning methods.
In the present study, LAI is estimated by capturing plant shadows using UAV imagery. This method is compatible with previous methods because they are all based on the need to find a quick, non-destructive, and accurate way of measuring leaf area. However, as previously stated, these methods usually require measurements to be carried out around midday to minimise the effect of shadowing between vine rows.

Experimental site
A large-scale field experiment was carried out in a vineyard (cv. Cabernet-Sauvignon) located in 'Zamadueñas estate' (coordinates: 41.7013º N, 4.7088º W, Valladolid, Spain), which belongs to the Agricultural Technology Institute of Castilla y León (ITACyL). According to the WRB classification (FAO), the soil is medium-coarse textured (FLc) Calcaric Fluvisol + (FLe) Eutric Fluvisol (Nafría et al., 2013). The climate is Csb (temperate with dry or temperate summer -Köppen-Geiger Climate Classification) and is characterised by dry summers and mild, wet winters (AEMET, 2011). Table 1 summarises the meteorological data from the study site collected by 'VA101 Finca Zamadueñas' weather station. The information is available at http://www. inforiego.org/. The vineyard comprised a vertical shoot positioned trellis (VSP), and the vines were spur pruned (with bilateral cordon training), with eight spurs per plant, two buds per spur, 2.5 m x 1.2 m row and plant-spacing respectively, and 1.8 m average vine height. The orientation was NE-SW (northeastsouthwest), 35 º to the north. The soil was kept free of any weeds that would have affected image processing (Fountas et al., 2014).

Field determination of LAI
In June 2019, the vines were removed from the vineyard. The vines were geo-positioned ( Figure 1a) using a GPS Triumph-2 JAVAD GNSS model (Triumph-2, JAVAD GNSS Inc, San Jose, California, USA; Figure 1b) with centimetre accuracy to mark the vines in the field to be removed within the study area. The GPS TRIUMPH-2 has 216 channels of dual-frequency GPS and GLONASS and could be connected to the mobile phone via Bluetooth and Wi-Fi to access the local GNSS Reference Station Network. Each grapevine was cut in the lower-middle part ( Figure 1c) and all the material was extracted from the vineyard.
The real LAI (Leaf Area Index) was determined on a sample of 36 vines. The total leaf area of each removed vine was measured using the EasyLeafArea application (Easlon and Bloom, 2014). The EasyLeafArea app automatically calculates leaf area from green leaf and red scale areas.
It is important to note that the shadow area analysed in this study represents PAI (Plan Area Index), as it captures all plant structures and not only the leaves.
However, the term LAI is used in the present manuscript, because, like other indirect methods based on images, the method employed in this study does not distinguish between leaf and nonleaf material in the analysis; therefore, LAI was used for consistency with real leaf area values obtained using the EasyLeafArea application, and from previous literature (Fuentes et al., 2014;Kalisperakis et al., 2015;De Bei et al., 2016;Towers et al., 2019;Beeri et al., 2020). In addition, there is a strong relationship between PAI and LAI from May to October. Doring et al. (2014) observed that estimated PAI did not differ substantially from directly measured LAI in a vineyard, showing a remarkably high correlation between PAI and LAI (R² = 0.93). Moreover, Arnó et al. (2013) found a very high correlation between PAI and LAI (R² = 0.99) using 'TAI' instead of 'PAI' since LiDAR does not distinguish between green and non-green elements.

Flight campaign/survey
The UAV images were acquired on 27 June 2019. Considering the average vine height and using basic trigonometry (Figure 2a), the sun elevation was fixed at β = 45 º. As a result, the flight was scheduled for 15:30 to 16:00 solar time (18:00 local time) using NOAA Solar Calculator (https://www.esrl.noaa.gov/gmd/grad/solcalc/), and under 1 Okta cloud cover conditions and azimuth α = 265 °. This time was used to maximise the vine shadow projection on the floor.

Unmanned Aerial System (UAS)
Before the UAV flight, a set of 12 ground control points (GCPs) were located in the vineyard and georeferenced using a real-time kinematic (RTK) GPS Triumph-2 JAVAD to improve the geometric accuracy of the image mosaicking process ( Figure  2b). Toffanin (2019) defines a Ground Control Point (GCP) as a position measurement made on the ground that can be set using existing structures like pavement corners or lines on a parking lot. A minimum of five GCPs are enough when they include points near corners and within the study area. In this study, clearly distinguishable field structures such as lampposts, maintenance holes and vineyard posts were used (Figure 2b).
An Unmanned Aerial System (UAS) composed of a 3DR quadcopter platform (3DR SOLO, 3D Robotics, Berkeley, California, USA, Figure 1d) was used to fly autonomously on a previously planned mission using open-source autopilot software: ArduCopter 3.3 (Ardupilot) installed on Pixhawk 2.0 mainboard (3D Robotics, Berkeley, California, USA). The image was acquired using a MAPIR® low-cost RGB+RGN system (Survey3W, MAPIR® Inc, San Diego, California, USA) commanded by a drone's flight controller and composed of two cameras with 12-megapixel rolling shutter CMOS sensors. These sensors had a focal length of 3.37 mm with a horizontal field of view (HFOV) of 87 ° and -1 % extreme low distortion glass lens and dual-band filter made of silicon, which is sensitive in the Visible and Near-Infrared spectrum from about 400 to 1200 nm. The sensor captures Blue (450 nm), Green (550 nm), Red (660 nm) and Near Infrared light (850 nm). Each camera had a f/2.8 aperture and was factory calibrated. The images were automatically geotagged using an attached GNSS system (NEO-m8, u-blox, Thalwil, Switzerland). Moreover, it is crucial to correctly configure the camera settings (shutter speed, ISO and EV) so that no pixels reach the maximum pixel value. If a pixel is higher than the maximum value, the information will be lost. Pixels have a value that ranges from a minimum to a maximum based on the image bitrate. The higher the bitrate, the more information that can be stored in the image. A sensor captures each image in a RAW format and then saves the RAW or converts it to a more standard format (typically by compressing it). Based on the information provided by the manufacturer, the Survey3 cameras capture 16bit RAW photos per RGB/RGN channel, which means that there are 16 bits (65,536) pixels and that a pixel value ranges from 0 to 65,535. In this work, since the reflectance of light was captured to analyse differences between pixels to identify shadows, RAW was the used format. In order to keep the pixels from reaching the maximum value mentioned above, the camera was configured as follows: shutter speed: 1/1000 seconds, ISO:50, and exposure: +0.0. The Ground Sample Distance (GSD) estimated by the software (Mission planner) for the used camera profile (SURVEY 3W) and a flight height of 22 was 1.01 cm/px.

Orthomosaic generation
The 16-bit raw images were first converted to tiff format using MAPIR® Camera Control (MCC) software. Each pixel was corrected using a known reflectance value by capturing a photo of the 'MAPIR® Camera Reflectance Calibration Ground Target Package' just before the survey; the photo contained four targets that were measured along incremental wavelengths by a calibrated spectrometer (based on the information provided by the manufacturer, reflective measurements were made along incremental wavelengths of 1 nm and from 350 nm to 1100 nm using multiple Shimadzu spectrophotometers with an integrating sphere).
The pixel values of the captured target image were then compared with the known reflectance values of the targets. Using this information in MAPIR® Camera Control (MCC), the pixel values were transformed, and thus the survey images were calibrated. Then, the images were imported into image mosaicking software (Agisoft Metashape 1.5.2, Agisoft LLC, St Petersburg, Russia) based on the structure-from-motion method (SfM) algorithm (Westoby et al., 2012) to generate the orthophotos. The Exif metadata of each image from the GNSS was first used to help in the image alignment process. Then the locations of GCPs from the JAVAD GNSS were manually identified and added to the aligned images to optimise camera positions and orientation data and in turn to improve orthophoto accuracy.

Data analysis
The data analysis pipeline is presented in Figure  3. Shadows were detected using remote sensing data via three methods: i) Image reclassification (sensitivity method/manual), ii) K-means, and iii) Random Forest. Each plant area was then cropped and the shadow within the area was measured and correlated with real LAI.

Grid definition
A vectorial grid was developed in order to isolate the shadow of each vine. First, azimuth and row orientations were used to calculate the grid angles, resulting in a parallelepiped with angles of 95 ° (360 -265 = 95 °) and 35 ° ( Figure 4a) and four sides which were 1.2 and 1.8 m long. The next step was to replicate the grid in the whole vineyard (Figure 4b), excluding a 0.7 m width line corresponding to the vegetation.
Once the grid was designed (Figure 4c), the points taken in the field using GPS were used to indicate which grid area corresponds to each vine. Finally, the image was cropped crop using the grid to isolate the area that belonged to each vine.

Pixel value analysis
Once the image corresponding to each vine was isolated, a random sampling strategy was manually carried out at various points where there was shade, obtaining the R, G, B and NIR values.
Subsequently, a sensitivity threshold was defined by the ' a ' coefficient: f(x) = (1 ± a) x whereby the detected shaded area is a function of the band value x modified by the sensitivity coefficient a. As a result, a range was defined by an upper and lower limit, depending on the sensitivity value. Within the range, the value is considered 'shadow positive'. It is important to consider that not all shades have the same intensity due to the differences in vegetative development modifying intercepted radiation (Zheng andMoskal, 2009, Baeza et al., 2010); therefore, in order to establish an optimal sensitivity value, a visual analysis was performed of the threshold effect on the image. The coefficient varied from 0.1 to 1(10 % to 100 %) at intervals of 0.1 (10 levels in total), covering the whole range of shadows: from low to high shadow intensity. This process was carried out for each band, reclassifying the image with the values corresponding to the range of shadows, where '1' corresponds to 'shadow' and '0' indicates 'absence of shadow' (Figure 5a). Subsequently, the bands were combined to create an RGB or RGN product comprising a combination of the corresponding bands ( Figure  5b); the corresponding shaded area was calculated from the total area of the vine and a product with the following shaded area was obtained: ( , , ) = (1 ± ) · (1 ± ) · (1 ± )

-EQUATION 1
where x, y and z are the band values (RGB or RGN) and a, b and c are the sensitivity coefficients for each band (10 levels ranging from 0.1 to 1.0). For each RGB and RGN product, 10 3 = 1000 cases were calculated. This information was validated with the real LAI values in order to find the closest relationship between the detected shadows and the LAI. The best result was used as input for the machine learning algorithms.
Finally, the potential of the method is shown by the classification of the vines depending on shaded area into three LAI levels (high, medium, and low), thus demonstrating that the method could be useful for real applications, such as zoning for differential management, irrigation, or fertilisation.

K-means clustering (Unsupervised machine learning)
K-means clustering is an unsupervised machine learning method that classifies the input data objects into multiple classes based on their inherent   distance from each other, thus dividing the input data set into k clusters previously defined by the user (Hung et al., 2005). The clustering problem can be formulated as an optimisation problem, described by: : minimise ( , ) = 0 0 !" 3 ! , where w ij = 1 implies object x i belongs to cluster G(j), and d(x i , μ j ) denotes the Euclidean distance between x i and μ j for i = 1,…, n and j = 1,…, k. (Pérez-Ortega et al., 2020). The standard version of the algorithm consists of four steps: 1) centroids (k points) are randomly generated in the space, 2) each point is assigned to its closest centroid, according to the distance from all the centroids, 3) new centroids are calculated using the mean value of the objects that belong to each cluster, and 4) the process is repeated from step 2 until equilibrium is reached (i.e. when the number of points remains stable within each cluster).
By using k-means it is possible to classify pixels automatically without training samples. The number of k clusters was set from 2 to 6 to determine the effect of different clusters in the shadow classification process. Subsequently, the cluster with the best visual coincidence with the shadow locations was assigned to the shadow class. The initial location of the centroids was set randomly, and the maximum number of iterations was set to 1000.

Random Forest classification (Supervised machine learning)
Random Forest is a supervised machine learning method that uses a decision tree for classification and prediction. The algorithm fits many classification trees to a dataset and then combines the predictions from all the trees. Each tree can be computed separately from other trees, because each tree is independently constructed using a bootstrap sample of the data set (Kuhn, 2008).
The Random Forest algorithm works in four steps: 1) many bootstrap samples from the data are selected, 2) a classification tree is fit to each bootstrap sample, 3) each tree is used to predict the out-of-bag observations, and 4) the predicted class of an observation is calculated by majority vote of the out-of-bag predictions for that observation (Cutler et al., 2007). The observations in the original dataset that do not occur in a bootstrap sample are the out-of-bag observations. In this study, the number of trees was set to 500, and the training set was divided into four classes: soil, shadows, shaded vegetation, and vegetation.

Model comparison
Calibration (cross-validated, 50 % split) and validation were used to assess and compare the ability of each model to predict actual LAI. The coefficient of determination (R²), root mean square error (RMSE) and mean absolute error (MAE) were used to define the model's performance.
Some authors recommend MAE as the most natural measure of average error magnitude instead of RMSE, since measures of average error (such as RMSE), which are based on the sum of squared errors, are functions of the average error (MAE), the distribution of error magnitudes (or squared errors) and n 1/2 ; therefore, they do not describe average error alone and MAE is less sensitive to the effect of outliers than RMSE as an indicator of model performance (Willmott and Matsuura, 2005). However, RMSE is also indicated, since it is commonly used in Remote Sensing literature All image and data analyses were carried out using AutoCAD® (version 2021 R.47.X Autodesk, Inc. San Rafael, California, USA), QGIS (version 3.14.X, QGIS developer team 2020), customised codes written in R (version 3.6.X, R Core Team 2019) and R packages stat, caret and randomForest obtained from the Comprehensive R Archive Network (CRAN).

Shadow pixel analysis
The sensitivity of the thresholds for the shadow detection of each band was adjusted separately (Figure 6a). Most of the bands have a maximum correlation with LAI at a sensitivity level of around 30 %.
Once the maximum R² value is reached, the correlation decreases progressively. As previously explained, the band combinations were analysed using sensitivity values from 0.1 to 1.0, at intervals of 0.1 (10 levels in total). A strong linear correlation can be observed between the real LAI and the area obtained from the shadows. For the RGB product, using a sensitivity for all coefficients of 0.3, the maximum R² value was 0.57. For the RGN product, the maximum was R² = 0.64 (Table  2), using a sensitivity of 0.4 (Figure 6b).

Shadow pixels classification
RGN was the input selected for the machine learning models, since it achieved the maximum correlation with real LAI. Figure 7 shows k-means clustering results for each cluster number. Since the input data is automatically classified using an unsupervised machine learning method, the products obtained differed substantially; however, for any number of clusters, and beginning from k = 2, k-means was able to identify shadows. Regarding Random Forest, the algorithm was  able to identify the four required classes, and the shadows were appropriately classified (Figure 8).

Observed vs predicted LAI values
First, RGB and RGN datasets were compared, because they had the highest correlation with LAI.
Once the model was defined, the predicted LAI was compared with the observed LAI (Figure 9).

DISCUSSION
Aerial orthophotography provides information about the canopy size in two dimensions, but shadows comprise a projection in the third dimension, providing additional valuable information. The shadows depend not only on the shape of the vegetation but also on the lighting. Therefore, it is imperative to have an optimal light source (sunny days). Shadows vary throughout the day depending on the position of the sun; therefore, it is possible to plan the flight at the optimum time when the sun is in the desired position (e.g., in the afternoon). However, images can be captured at other times of the day, as long as there are shadows, but the accuracy will probably be related to the size and quality of the shadows.
Our results show that it is possible to accurately plan a flight at a given time by calculating the azimuth and elevation of the sun so that the information extracted from the shadows can be maximised. This study shows that, in terms of shadow extraction, the afternoon is a good time to take images; however, the exact time should be adjusted depending on vineyard characteristics. At the optimal time, shadows effectively cover the ground between rows, but they are neither affected by the vegetation in the adjacent line, nor by the lack of precision due to an absence of shade. A balance can thereby be achieved, and the maximum amount of information can be obtained from the shadows generated on the ground. The optimum time for taking the image will be determined by the height of the vines and the distance between the rows of plants.
Overall, the results showed a significant positive relationship between the plant shadows and the LAI. According to the results of the RGB and RGN analysis, a correctly identified shadow resulted in a positive and close significant relationship of up to R² = 0.74. Figure 6 shows a maximum correlation between the shadows and the LAI at around 30 % sensitivity: if the sensitivity is too low, the plant's shadows cannot be easily detected, whereas if it is too high, the algorithm will identify pixels as shadows when they are not ( Figure 11). Therefore, errors can occur by overestimating or by underestimating the shadows. However, the slope is much more stable after reaching the maximum value of R²; thus, in order to detect shadows, it is better to overestimate than to underestimate FIGURE 11. Effect of the sensitivity level.
the value of the pixels. In this study, the models tend to underestimate at high LAI values and overestimate at low LAI values (Figures 9 and  10), which is probably due to an increase in leaf overlap in big canopies.
When comparing RGB and RGN products, the results show RGN to be a better approach, probably because NIR can identify the shaded vegetation, thus helping to discern it from shaded soil.
Regarding machine learning methods, Random Forest improved the RGN accuracy for predicting LAI reaching R² = 0.76; this was matched by k-means with k = 6 (R² = 0.76). The rest of the k-means classifications showed good accuracy, even in the worst case, with R² ranging from 0.62 to 0.75. The best accuracy was obtained for k = 3 and k = 6, and the worst was found for k = 4, showing that in k-means, a higher number of clusters for classification is not strictly related to a better accuracy level.

Method considerations
In order to ensure that this method works optimally, several essential factors must be taken into consideration, for example:  i) Sun position, which influences the shape of the shadows; it is, therefore, essential to consider azimuth and elevation.
 ii) Grapevine management. Trimming the vegetation will modify the vineyard shadows. Depending on the phenological stage, new leaves may cover the gaps. Moreover, errors in shadow detection can occur if the vegetation is not within the trellis for mainly two reasons: → a. The vegetation covers the shade, which therefore goes undetected (Figure 12a) → b. The vegetation itself creates shaded areas onto itself (Figure 12b). NIR can help identify them since the values of shade in NIR are different from that in RGB, thus explaining why it is better to use NIR for the detection of vegetation.
 iii) Weed management. Weeds can add noise to the image and render shade detection difficult.
 iv) Vineyard orientation. The optimum orientation for late-day image capture in the study zone is close to SSE-NW, since it is the orthogonal plane regarding the path of the sun. The applicability of this method decreases as the difference between the real plane and the optimum plane increases. However, this should not be an issue in viticulture, since rows are generally oriented close to north-south to maximise light interception on both sides of the canopy for a part of the day (Jackson, 2020).

Limitations of the method
There are some limitations and important points to raise related to the vineyard characteristics and the method itself:  i) The vineyard vegetation growing in the assigned area extends along the whole length of the vertical trellis, so that it is possible for the vegetation of one grapevine to mix with the adjacent one and thus be wrongly assigned to the other vine. However, this is inherent to the farming system, and it is not a source of error specific to this method.
 ii) If the vineyard dimensions are not exact, by creating an equal and regular mesh of parallelepipeds at certain points, it is possible to under-or overestimate the vegetation due to the error introduced by the designated area. This could be overcome by delineating a boundary (box) around and above each vine trunk using a centroid approach (with prior knowledge of vine spacing).
 iii) Proper overlap is extremely important in order to obtain good quality data and avoid defects such as blotchy artefacts or errors in image alignment.
 iv) The effect of the slope on the shadows could be a problem; however, the effects of topography can be addressed by proper flight/ mission planning and the use of DEMs/DCMs. The terrain slope effect was ignored in this work because the vineyard was on flat terrain.
 v) Light intensity. This study was carried out under good light conditions, and according to our results, a clear sky is needed. If the light intensity varies because of clouds or diffuse radiation, the shadows will not be defined, and the results will be affected.
 vi) Structures, such as VSP wires, poles or irrigation pipes project shadows which could get mixed up with the shadows of the vines.
 vii) The method might not work in highdensity vineyards/orchards due to the mutual shading of vines/trees.

Comparison with other Remote Sensing methods to estimate LAI
As previously discussed, each method has its advantages and disadvantages. Table 3 briefly compares some previously described methods.
The comparison considered similar factors to those considered by Jin et al. (2021), including 'economic cost' in an attempt to represent the cost of the operation, including all the equipment involved and 'time investment', which is the estimated time required for data collection in the field.
In terms of time investment, our proposed method is much quicker than methods that employ specific instrumentation such as PCAs or ceptometers, which are carried on foot to survey plants one at a time (Johnson and Pierce, 2004;López-Lozano and Casterad, 2013;White et al., 2018). LiDAR must be mounted on a ground vehicle (Arnó et al., 2013), as is generally the case for equipment in any field survey method (Diago et al., 2012;Fuentes et al., 2014;Towers et al., 2019). Compared to other methods that use similar technologies (UAV/Airborne + camera), our method involves similar input time (Hall et al., 2008;Kalisperakis et al., 2015;Comba et al., 2020) and is quicker than SfM methods, which require more images (Mathews and Jensen, 2013). However, it is still necessary to go to the vineyard to collect the data; therefore, satellite methods are even quicker (Johnson, 2003;Beeri et al., 2020).
The proposed method is low-cost as basic UAV and normal RGB+RGN cameras can be used. The best performance was obtained with the RGN camera, which is slightly more expensive than an RGB camera, such as that used in Diago et al. (2012) and Fuentes et al. (2014); however, it is cheaper than the multispectral sensors used by Towers et al. (2019) and Comba et al. (2020) or the hyperspectral sensor used by Kalisperakis et al. (2015). Specific sensors such as PCA or ceptometer are also more expensive. Moreover, when considering the cost of the platform, the UAV used in the study is much cheaper than a ground vehicle (Arnó et al., 2013).
As regards accuracy, our method is in line with other suitable methods in studies reporting similar correlation values. Furthermore, to our knowledge, this is the first time a method of this kind has been developed and will therefore likely be improved in further experiments, which has been in the case in other studies; for example, Mathews and Jensen (2013) reported R² > 0.5 accuracy, but in subsequent studies, Kalisperakis et al. (2015) and Comba et al. (2020) reported higher accuracy (R² > 0.8).
In summary, the method described in this article is cheaper and less time-consuming than other methods for calculating LAI and similar levels of accuracy can be obtained (Table 3). Lastly, other methods usually require the experiment to be carried out close to solar noon, or they try to at least avoid shadows; the present method thus allows the flight to be scheduled during other periods of the day than solar noon, such as morning or afternoon, enabling pilots to extend their working day.
It should be noted that this comparison only aims to highlight the potential and accuracy of existing methods versus the usefulness of the method proposed in this paper. The results presented in the literature are not fully comparable, because they have been obtained from studies carried out under different experimental conditions. The reference LAI used was not obtained in the same way in all the experiments. In the present study, real LAI was used, which was measured destructively to validate the estimated results, but other authors have used LAI estimated by other methods, such as with a Li-COR PCA (Johnson, 2003;Fuentes et al., 2014) or an AccuPAR LP-80 ceptometer (Mathews and Jensen, 2013).

Potential applications
This method can be applied to other woody crops, such as olive or walnut trees, where mutual shading of consecutive trees has a lower effect, depending on the planting distance. However, in these cases (bigger canopies), the internal leaf overlap is higher, an effect which would need to be considered in the application/calibration of the proposed method. In addition, given that the vegetation of one plant can mix with the adjacent in vertical-trellis trained crops, this method may be even more effective when applied to non-trellis trained crops.
Images could also be taken several times a year to monitor the growth of the vegetation (for plant vigour). From the image history of the crop, it would be possible to create a temporal model of the evolution of the shadows, and therefore of the LAI, and thereby adjust vineyard management accordingly (e.g., irrigation, treatments, and leaf removal). However, the number of woody structures covered by the leaves will certainly affect the accuracy of the method.
Finally, to illustrate the real applicability of the method, a cluster analysis was developed: the dataset was split into three groups, with each plant being classified according to shaded area on three LAI levels (high, medium, and low), and 0.88, 0.75 and 0.63 m 2 of soil shaded area as centroids and 1.79, 1.54 and 1.31 for high, medium, and low LAI respectively. Figure 13a shows three different The R 2 values are the values reported by the authors. The 'Economic cost' and 'Time invested' values are based on the methods described by authors compared to the present study method. In 'Economic cost': 'low' total cost up to 3,000; 'medium' up to 6,000 $; 'high' more than 6,000 $. clusters in which the vines are grouped according to LAI; k-means classification was used for this, but other machine learning methods could be applied. Figure 13b shows the result of the classification on a map, showing that this method could be used for, for example, zoning the vineyard according to different irrigation or fertilisation management.
The previously mentioned factors prove that shadows can estimate the LAI. Furthermore, the results show that, by using simple equations and machine learning methods and by controlling the limitation factors, leaf surface area can be estimated from the shadow of the plants. Therefore, all advantages and limitations considered, the shadow captured from UAV images is a good estimator of LAI, provided that the method used to obtain the shadows is the appropriate one.
Further research is needed to assess the temporal stability of the relationship between plant shadows and LAI, and the applicability of the method to other woody crops and farming systems or vineyards planted with cover crops. Moreover, it would be interesting to explore whether this method can be used to estimate other vineyard parameters or to assess the possibility of extracting more information about canopy structure and vegetation density; for example, Figure 11 shows gaps in the detected vegetation using different sensitivity levels (a, b and c) due to differing shadow intensities for each pixel, indicating that it may be possible to use this method for extracting more information about canopy structure and vegetation density. Other authors (Zheng and Moskal, 2009;Baeza et al., 2010) have studied the relationships between the incident radiation intercepted by the vine, LAI, and aboveground biomass, concluding that the measurement of the gap fraction is a way of analysing canopy structure and that it can be parameterised by LAI and leaf angle distribution. Additional research is required in this area.

CONCLUSIONS
This study shows a simple and effective method for estimating LAI in vineyards using machine learning tools and projected shadows captured by a UAV. It is a combination of methods known for centuries and modern image recognition techniques.
A strong positive relationship was observed between the shaded area and the leaf area of the vines, thus resulting in a high level of accuracy for the LAI estimation. (R² = 0.76, RMSE = 0.160 m 2 m -2 and MAE = 0.139 m 2 m -2 ). However, to accurately estimate LAI, shadow areas must be correctly identified, otherwise errors can occur due to over-or underestimating the shadows. In addition, it is essential to manage the influence of parameters on the method (e.g., the position of the sun, vineyard management or weed management) to avoid noise in the image.
The results of the study are useful because LAI is an essential parameter for vineyard management and zoning. Moreover, this method is fully compatible with other Remote Sensing methods and can be incorporated into a flexible working day.