Comparison between satellite and ground data with UAV-based information to analyse vineyard spatio-temporal variability

Currently, the greatest challenge for vine growers is to improve the yield and quality of grapes by minimizing costs and environmental impacts. This goal can be achieved through a better knowledge of vineyard spatial variability. Traditional platforms such as airborne, satellite and unmanned aerial vehicles (UAVs) solutions are useful investigation tools for vineyard site specific management. These remote sensing techniques are mainly exploited to get the Normalized Difference Vegetation Index (NDVI), which is useful for describing the morpho-vegetational characteristics of vineyards. This study was conducted in a vineyard in Tuscany (Italy) during the 2017, 2018 and 2019 seasons. Ground data were acquired to detect some agronomic variables such as yield (kg/vine), total soluble solids (TSS), and pruning weight (kg/vine). Remote sensed multispectral images acquired by UAV and Sentinel 2 (S2) satellite platform were used to assess the analysis of the vegetative variability. The UAV NDVI was extracted using both a mixed pixels approach (both vine and inter-row) and from pure canopy pixels. In addition to these UAV layers, the vine thickness was extracted. The aim of this study was to evaluate both classical Ordinary Least Square (OLS) and spatial statistical methods (Moran Index-MI and BILISA) to assess their performance in a multi-temporal comparison between satellite and ground data with UAV information. Good correlations were detected between S2 NDVI and UAV NDVI mixed pixels through both methods (R2 = 0.80 and MI = 0.75). Regarding ground data, UAV layers showed low and negative association with TSS (MI = 0.34 was the lowest value) whereas better spatial autocorrelations with positive values were detected between UAV layers and both yield (MI ranged from 0.42 to 0.52) and pruning weight (MI ranged from 0.45 to 0.64). The spatial analysis made by MI and BILISA methodologies added more information to this study, clearly showing that both UAV and Sentinel-2 satellite allowed the vigour spatial variability within the vineyard to be detected correctly, overcoming the classical comparison methods by adding the spatial effect. MI and BILISA play a key role in identifying spatial patterns and could be successfully exploited by agricultural stakeholders.


INTRODUCTION
In the last years, the agriculture community has entered the world of technology, which is an increasingly indispensable resource in crop management. The need for ever more precise information aimed at supporting stakeholder's decisions, led to the use of remote sensing (RS) data, which are timely, synoptic, cost-efficient and repetitive information for crop management purposes (Atzberger, 2013;Di Gennaro et al., 2019;. The application of these technologies and methodologies for vineyard site specific management represents the core of precision viticulture (PV) . Improving vineyard yield and grape quality through a proper knowledge of the vineyard spatial variability and thereby reducing costs and environmental impact is one of the current challenges for vine growers (Song et al., 2014;Silvestroni et al., 2018;Khaliq et al., 2019). This crop has a very wide international diffusion since it can be cultivated in very different conditions due to a great number of varieties with different behaviour and numerous training systems. In particular, this crop has a macroscopic ecophysiological response to pedomorphological and climatic conditions that often vary within the same vineyard . These peculiarities make remote sensing a fast and useful solution for spatio-temporal monitoring (Borgogno-Mondino et al., 2018). The most common remote sensing platforms used in viticulture are Unmanned Aerial Vehicles (UAV), airborne sensors and satellites, which can be equipped with different types of optical sensors to analyse plant response in a wide spectral range, such as RGB, multiand hyper-spectral, thermal infrared and LiDAR (Hall et al., 2002;Mathews and Jensen, 2012;Pádua et al., 2019;Sozzi et al., 2020). Among these remote sensing platforms, airborne, but above all UAV, are characterized by low operational costs, excellently timed and flexible surveys and high ground spatial resolution (even a few centimetres) providing detailed information of vine features within the field Jay et al., 2019;Pichon et al., 2019;Sozzi et al., 2020). However, these two systems allow information to be obtained on limited areas and require long and specialized post-processing that can be difficult for the farmers to interpret (Candiago et al., 2015, Sozzi et al., 2020. Although satellites are affected by cloud cover, they are a fundamental tool for long-term and extensive crop monitoring. Nowadays, a large range of satellite sensors regularly provide free datasets, thus promoting satellite technology for many agricultural applications (Atzberger, 2013;Yang et al., 2017;He et al., 2018;Khaliq et al., 2019).
Remote sensing in PV is mainly used for reflectance spectroscopy, an optical technique based on reflectance measurements of the incident electromagnetic radiation at different wavelengths, in particular in the visible, near infrared and thermal infrared regions . This data allows considerable information to be obtained, useful for describing the morpho-vegetational characteristics of vineyards (Johnson, 2003;Santesteban et al., 2013;Matese et al., 2019). This information includes a wide set of spectral vegetation indices based on different electromagnetic bands, useful for describing plant physiology (Magney et al., 2017;Pádua et al., 2019). The well-known Normalized Difference Vegetation Index (NDVI) (Rouse et al., 1973) is one of the most widely used. It makes use of the different responses of vegetation to the visible (red) and near infrared (NIR) spectral bands that are closely related to crop vegetative and productive features (Romboli et al., 2017;Khaliq et al., 2019). The presence of grass cover, bare soil or shadows may strongly affect the computation of spectral indices, leading to a wrong crop status evaluation. In order to solve this issue, some algorithms focused on pixels classification have been developed to perform pure canopy pixels extraction by means of crop-specific masks (Atzberger, 2013;Peña et al., 2013;Pérez-Ortiz et al., 2015;Senthilnath et al., 2017;Cinat et al., 2019;Khaliq et al., 2019).
The high heterogeneity of vineyards can have different impacts on vine physiological response, with direct consequences on grape quality and yield. It is therefore clear that vineyards need specific agronomic management according to the spatial variability within the field (Proffitt et al., 2006). The use of statistical methods, which take into account the spatial structure of vineyards to compare ground measurements with information derived from remotely-sensed data, could be a robust and powerful tool to evaluate the practical implications of vineyard variability, and represents a valuable trigger to increase the spread of PV techniques . Many studies base their data comparison on either visual evaluation or classical Ordinary Least Square (OLS) regression (Baluja et al., Laura Pastonchi et al. 2012;Yang and Everitt, 2012); however, both approaches have several shortcomings, as visual inspection can suffer from subjectivity errors, and the OLS is not well-suited to spatially structured data because OLS scores can be biased and their significance inflated . There is a need to use other methods more suitable for this kind of study. Within the geostatistics solutions, there are two methods that take into account autocorrelation, represented by the Moran Index (MI, Moran, 1950), and the Bivariate Local Indicators of Spatial Autocorrelation based methods (BILISA, Anselin, 1995). These spatial dependence statistics are more objective means of quantifying a pattern in terms of patchiness (Miller et al., 2007). The MI measures the global correlation between data in the comparison process; instead, BILISA focuses on degree of spatial clustering and dispersion between two different variables in nearby locations (Anselin and Rey, 2014). Currently, there are not many studies that apply these statistical approaches to data analysis for agricultural purposes; there is also a lack of research concerning local spatial autocorrelation indices used to evaluate the spatio-temporal patterns in fields .
The purpose of this study is thus to use the aforementioned classical and spatial statistical methods (OLS, MI and BILISA) to evaluate their performance in the comparison between satellite and ground data with information acquired from a UAV platform. The aim is to identify a methodology able to overcome the limits of the classical statistical approach and provide objective means of quantifying the variability within a vineyard. This analysis was conducted for 2017, 2018 and 2019 seasons in an experimental vineyard located in an area of excellence of Tuscan viticulture (Italy): the Chianti Classico. Ground data were acquired during the harvest periods in order to detect some agronomic variables. Through UAV flights and S2 satellite, multispectral images were obtained (S2 images were acquired by selecting those closest to the UAV flight date) and the NDVI was calculated. In addition to these layers, the vine thickness was extracted from UAV data. The statistical framework was thus used to assess, within a multi-temporal context, the degree of similarity of these different datasets, providing a clear demonstration of how remote sensing can be exploited by winegrowers for crop management.

Study area and ground measurements
The research took place in an experimental vineyard planted in 2008 and located in Castellina in Chianti within the Chianti Classico denomination, which is an area of excellence of Tuscan viticulture (Italy) (Figure 1). Sangiovese  (Vitis vinifera L.) vines were trained to a vertical shoot-positioned trellis system with spur-pruned single cordon. Vine spacing is 2.2 m × 0.8 m (inter-row and intra-row, respectively), rows are NW-SE oriented and the vineyard is on a slight southerly slope at 355 m above sea level. Pest, soil and canopy management were performed following the farm practices.
Six flight campaigns (referred as F1-F6, Figure 1) were performed during 2017, 2018 and 2019 seasons, and ground truth data, related to productive and vegetative parameters, was collected at the end of each season through destructive sampling in three vigour zones representative of vineyard vegetative variability ( Figure 1).
A preliminary vineyard vigour assessment was done using Sentinel 2 imagery. NDVI maps obtained from S2 surveys at the beginning of June 2017 allowed to delineate a zone with high vigour (HV-green circle) and another with low vigour (LV-red circle) in the north and middle part of the vineyard respectively, whereas a heterogeneous area with mixed vigour (MVorange circle), due to the co-presence of high and low vigour vines, in the south part of the field. NDVI maps obtained with the first UAV flight (F1) campaign confirmed the presence of those vigour zones. At the end of each season, 18 sampling vines were selected within each vigour zone, for a total of 54 vines. Yield (kg/vine) and total soluble solids (TSS) (°Brix) were measured for each vine using a field scale and a hand-held optical refractometer respectively. Pruning weight in terms of total shoot fresh mass (kg/vine), has been used as a proxy of vine vigour. This was determined in the field for each vine during the dormant period at the end of the growing seasons.

UAV data acquisition and processing
During each flight campaign, multispectral images were acquired through an open-source UAV platform, consisting of a modified multirotor Mikrokopter (HiSystems GmbH, Moomerland, Germany), described in detail in Matese and Di Gennaro (2018). The flight altitude, identified to ensure an optimal characterization of vineyards, is 50 m (AGL), providing 3.5 cm/pixel ground resolution. The UAV was equipped with a Tetracam ADC Snap multispectral camera (Tetracam, Inc., Gainesville, FL, USA), which is able to acquire data in the green (520-600 nm), red (630-690 nm), and near infrared (760-900 nm) bands. Camera settings were set to a fixed exposure with automatic trigger at 2 s frequency. Additional characteristics of this camera are described in Matese et al. (2016). The multispectral images were then mosaicked using Agisoft Metashape Professional v.1.6.0 (Agisoft LLC, St. Petersburg, Russia). This Structure from Motion (SfM) processing software has allowed a dense cloud of the vineyard to be obtained, and then a high-resolution NRG orthomosaic (Nir, Red, Green), a multilayer raster as defined in Cinat et al. NDVI unfiltered (mixed pixels, i.e. both vine and inter-row), NDVI filtered (vine-only pixels) and vine thickness, which is the section width of vine rows. First of all, a filtering procedure of pure canopy pixels was done by removing inter-row and mixed pixels. The basis of this procedure was that vine rows have a greater height from the ground and can easily be discriminated by Otsu's global thresholding, an algorithm that allows two different zones to be detected: vine rows and ground (Matese et al., 2016). This step produced a NRG mask layer with the same size and pixel resolution as DSM (Figure 2c), where each pixel was classified as canopy or not canopy (soil and shadows). The digital numbers (DN) of the NRG orthomosaics were then converted to reflectance through a radiometric calibration with the radiance detected by the sensor (Kelcey and Lucieer, 2012) on three OptoPolymer (OptoPolymer -Werner Sanftenberg, Munich, Germany) homogeneous and Lambertian reference panels (95 %, 50 % and 5 % reflectance). These panels were present on each flight. Subsequently, a rectangular polygon grid was defined as a function of spacing values of the vineyard (Figure 1), ensuring that each polygon included three plants. Therefore, the grid for the experimental site is 2.2 x 2.4 m ( Figure 2d). Finally, for each polygon, the radiometric information of all the mixed pixels was used to obtain the NDVI unfiltered data (Figure 2e). The NDVI filtered layer was instead obtained following the same approach but using a filtered mask to extract only the canopy pixels ( Figure 2f). That mask was also used to extract Laura Pastonchi et al. the third UAV product related to vine thickness, which was computed from the DSM following the 2.5D approach described in a previous paper (Di Gennaro and Matese, 2020) ( Figure 2g). As next step of this workflow, a punctual shapefile for each flight containing the average values of the three parameters calculated for each three vines grid (NDVI unfiltered, NDVI filtered and vine thickness) was created. Lastly, in order to compare the UAV products with ground sampling measurements (yield, TSS and pruning weight), the August UAV data (chosen as available in all the years) was joined to the ground sampling vines data in order to use this punctual shapefile in GeoDa software (Anselin et al., 2010).

Sentinel-2 data acquisition and processing
The Sentinel-2 multispectral images were downloaded from the official Copernicus Open Access Hub (www. scihub.copernicus.eu). Within this online database, cloud-free images were chosen according to their proximity to the UAV flight dates. These completely free data are 100 x 100 km 2 images, with a spatial resolution of 10 m, constituted by orthorectified Bottom-Of-Atmosphere (BOA) reflectance. The spectral bands B4 (635-695 nm) and B8 (727-957 nm), related to red and near infrared regions respectively, were extracted to obtain NDVI raster images with a 10 m pixel resolution ( Figure 2h). As an operational choice, the S2 spatial resolution was chosen as reference resolution for the satellite and UAV datasets. Thus, the three UAV products (NDVI unfiltered, NDVI filtered and vine thickness) were joined to the 10 m S2 NDVI map through QGIS software (Quantum GIS, http://www.qgis.org/), using the mean values per 10 m grid. Such an averaging methodology is commonly implemented for this type of analysis Di Gennaro et al., 2019;Sozzi et al., 2020). A final shapefile containing the four aforementioned values was thus obtained for each of the six surveys. Those shapefiles were then cut according to the vineyard profile, taking care to remove mixed border pixels.

Statistical analysis
For each flight campaign, the three UAV maps (NDVI unfiltered, NDVI filtered and vine thickness) were compared to S2 NDVI values using either classical or spatial statistical approaches, whereas the comparison between UAV layers and ground data was made only with the spatial analysis. This latter, unlike the classical method, considers the structure and geospatial variability present within the vineyard. All these comparisons were made using the GeoDa software. OLS regression was the classical method used to compare the obtained maps in terms of R 2 , and F-statistic was used to check the accuracy of the regression model. The adopted spatial methods are the statistical index I, developed by Moran (Moran, 1950), Moran's Index (MI), and Local Indicators of Spatial Association (LISA), developed by Anselin (Anselin, 1995), and based on the MI. Moran's Index is a measure of global spatial autocorrelation of the data within the whole analyzed area. Its value ranges from −1 to +1, where a positive autocorrelation means that there are similar values in nearby areas (clustered pattern); a negative spatial autocorrelation implies dissimilar values at nearby locations (disperse pattern) and the zero value (MI = 0) indicates the absence of spatial autocorrelation (random pattern). Univariate MI is the correlation between a variable X and the "spatial lag" of X, formed by averaging all values of X for the neighbouring polygons. Bivariate MI is, instead, a correlation between a variable X and a different variable in nearby areas. The fundamental difference is that, in the bivariate case, the spatial lag pertains to a different variable.
Both the Univariate and Bivariate Moran's Indexes have been used in this work. Local Indicators of Spatial Association (LISA) measures local spatial autocorrelation: the statistic is calculated for each areal unit present in the whole area, based on neighbouring units with which it shares a border. The result provides maps of local clusters with similar behaviour in the spatial arrangement of a given variable. The Bivariate LISA (BILISA) focuses on spatial clustering and spatial dispersion between features of a variable and another different variable in nearby locations (Anselin and Rey, 2014).
To perform these analyses, first of all, spatial weights were necessary because they allow the "nearness" or "proximity" between observations to be measured. GeoDa software gives the possibility to manage the spatial weights and, in this paper, the queen contiguity approach with the first-order of neighbours in a 3×3 matrix was chosen (Figure 3a). For this work, a spatial randomization approach that consists of from 99 random permutations allows to test whether the spatial correlation was significantly different from what would be expected in the case of spatial randomness. Significant spatial correlation, in this case, is detected between variables obtaining a Bivariate Moran's Index > 0 or < 0, and the null hypothesis of no autocorrelation is rejected due to Z score > +2.576 or < -2.576. Laura Pastonchi et al.

Classical statistical analysis-OLS
The R 2 values of the OLS regressions obtained for remote sensing datasets indicate generally good correlations between S2 NDVI and the three UAV layers (NDVI unfiltered, NDVI filtered and thickness, Table 1). Indeed, R 2 ranges between the minimum value detected for S2 NDVI vs UAV NDVI filtered of 2018 -F4 dataset (R 2 = 0.53), and the maximum correlation obtained for S2 NDVI vs UAV NDVI both unfiltered and filtered of the same F5 fight in the 2019 season (R 2 = 0.80). Except in this latter case, where two equal values were found, in Table 1 it is shown that the R 2 values associated of UAV NDVI unfiltered are always greater than those of the UAV NDVI filtered. The F-statistic used to check the validity of regression models always showed significant results (p-value < 0.001).

Univariate and Bivariate MI
The GeoDa software allowed definition of the correlation of both a single variable with a neighbourhood of the same variable (Univariate MI), and between two different variables in nearby areas (Bivariate MI). Through the first analysis (Table 2), the most clustered pattern (MI = 0.78) was detected in the 2018 season (F3 flight) for the UAV NDVI unfiltered, whereas the lowest MI value, indicative of a more random pattern, was identified for the UAV thickness of 2017-F2 flight dataset (MI = 0.39).
The comparisons between S2 and UAV data obtained using the Bivariate MI (Table 3) showed good correlations with positive trends; indeed, the lowest value was MI = 0.43 identified for 2017 season -F2 through the comparison of S2 NDVI with UAV thickness. The best correlation, indicative of a high spatial autocorrelation, appeared in 2018-F3 between S2 NDVI and UAV NDVI unfiltered (MI = 0.75). All these comparisons were significant because the Z-score value was always > +2.576 (for positive correlations). Regarding the comparison between S2 NDVI and UAV NDVI both unfiltered and filtered, it is possible to observe that in some cases the MI of S2 NDVI vs UAV NDVI filtered was greater than that obtained through the comparison between S2 NDVI and UAV NDVI unfiltered (Table 3).
The other comparisons performed through the Bivariate MI among ground parameters and    NDVI unfiltered (2018-F4), and between yield and UAV NDVI thickness (2019-F5). The best spatial autocorrelations were instead detected through the comparisons carried out between pruning weight and UAV dataset. The highest Bivariate Moran Index, MI = 0.64 (Table 4), was identified between pruning weight and UAV NDVI unfiltered, calculated for the August 2017 flight (F2). All the aforementioned comparisons are significant because the Z-score value is always > +2.576 for positive correlations and < -2.576 for negative correlations.

BILISA cluster map comparison
The spatial association cluster maps obtained using BILISA method both for a comparison between S2 NDVI and the three UAV layers of 2017, 2018 and 2019 seasons, and between ground information obtained within the three vigour zones (see Figure 1) and UAV data of August for the three years are shown in Figure 4 and Figures 5a, 5b, 5c. The obtained patterns illustrate the relationship between S2 NDVI at a given location and the average value at neighbouring locations of the three UAV datasets considered one at a time ( Figure 4). The same procedure has been then applied between ground data and UAV datasets (Figures 5a, 5b, 5c). The different colours correspond to the four types of local spatial autocorrelation, as already indicated in Figure 3 (Red: High-High; Dark Blue: Low-Low; Light Blue: Low-High; Pink: High-Low).
Regarding the first comparison (Figure 4), a general trend over the years can be seen,   characterized by an H-H association in the north and south-eastern part of the field, while an L-L association is always present in the central and western part. This trend is not so defined in the 2017-F1 dataset (Figure 4 a, b, c) and 2018-F4 comparisons (Figure 4 l, m, n), because more H-L (pink colour) and L-H (light blue colour) associations appeared between the analyzed variables. It should be noted that, for flight n° 3 of the 2018 season, the transition between H-H in the north and L-L in the south is clearer (Figure 4 g, h, i).
When referring to the cluster maps between ground data and UAV datasets, it was observed that each agronomic variable has a similar trend when compared with both UAV NDVI unfiltered, UAV NDVI filtered and UAV thickness of the three analyzed years (Figure 5a, 5b and 5c). A similar trend was also identified for the comparisons related to yield and pruning weight of the three years (Figure 5b and 5c). In particular, they are characterized by H-H correlations mainly within the high vigour area of the field (north part), and L-L correlations present in the other two parts (low vigour-central area and medium vigour-south area). It is also highlighted that in these last two areas, in some sampling points, H-L, L-H and also a few H-H correlations were found. Regarding the TSS variable vs UAV layers (Figure 5a), H-L correlation in the high vigour area, L-H in the central part and mixed spatial association in the south part of the field are present. For the 2018-F4 dataset it is possible to observe that a fair number of sampling points also showed H-H correlations.

DISCUSSION
The purpose of this study was to evaluate both classical and spatial statistical methods in a multi-temporal comparison between satellite and ground agronomic data with information acquired from a UAV platform, identifying the degree of similarity of these different datasets. The UAV NDVI information was also divided into unfiltered and filtered, this latter obtained using a filtered mask to extract only the information corresponding to the rows. This choice allowed a comparison of the spectral response of pure canopy pixels (filtered) and mixed pixels (both vine and inter row pixels, unfiltered) and the spectral information obtained from the S2 platform. The different colours correspond to four types of local Through the classical Ordinary Least Square (OLS) regression, the S2 NDVI maps compared with UAV NDVI unfiltered, UAV NDVI filtered, and UAV thickness showed generally good correlation with a positive trend in all the investigated seasons (2017,2018,2019); indeed, R 2 ranges from 0.53 to 0.80 (Table 1). The R² values obtained in this study were comparable to other works Sozzi et al., 2020), and the statistic used to check the validity of regression models always showed significant results (p-value < 0.001). The OLS analysis also highlighted that the correlation between S2 NDVI and UAV NDVI unfiltered is almost always higher than the one obtained for S2 NDVI vs UAV NDVI filtered, as reported by Di . This due to the same mixed pixel composition from which the NDVI index was calculated. Conversely, Sozzi et al.
(2020) reported a higher R 2 value for S2 NDVI vs UAV NDVI pure vine index but, observing the regression model between S2 NDVI and UAV NDVI mixed pixels, it showed a trend very close to a 1:1 line; this is indicative of the fact that these data are similar.
The use of a spatial statistical approach through the Univariate Moran Index (Table 2) allowed the spatial pattern of remotely-sensed variables to be identified. This analysis showed that the highest value, indicative of a clustered pattern, was detected on June 26 th , 2018 (F3 flight).
Every year, the analyzed vineyard is subjected to vine trimming to make it more uniform. The F3 flight, unlike others flights, was made only a week after this practice: this made it possible to identify, through the statistical analysis, a strongly clustered pattern.
The Bivariate Moran Index was instead used to identify the spatial autocorrelations between both Sentinel-2 NDVI and ground-truth agronomic parameters (TSS, yield and pruning weight) and the three UAV layers. The highest correlation was detected for the June 2018 dataset between UAV NDVI unfiltered and Sentinel-2 NDVI (MI = 0.75, Table 3). Unlike the OLS analysis, the Bivariate MI showed, in some cases, that the spatial autocorrelation detected for S2 NDVI vs UAV NDVI filtered is greater than that obtained through the comparison between S2 NDVI and UAV NDVI unfiltered. Regarding ground data, UAV layers showed low and negative association with TSS values (Table 4), and the lowest Bivariate MI of all the comparisons was detected between TSS and UAV NDVI unfiltered of the August 2018 dataset (MI = -0.34). This kind of correlation between TSS and vine vigour can be found in a previous study conducted by Fiorillo et al. (2012). Better spatial autocorrelations were detected through the comparisons carried out between yield and UAV layers (MI ranged from 0.42 to 0.52), and between pruning weight and UAV dataset (MI ranged from 0.45 to 0.64). All these comparisons showed positive autocorrelation, as highlighted by Romboli et al. (2017). Analysis of the BILISA cluster maps added more information to this study. Regarding the maps obtained for S2 NDVI vs UAV layers (Figure 4), a general trend was observed over the years. High UAV NDVI (both unfiltered and filtered) and UAV thickness were associated to high S2 NDVI (H-H correlation) in the southeastern and northern parts of the field; this latter zone is where high vigour (HV) was detected through ground measurements ( Figure 1). The areas characterized by low satellite NDVI were instead associated to low UAV layers values (L-L correlation) both in the central and western parts of the field, which are the two zones where low vigour (LV) and mixed vigour (MV) were identified. The cluster maps obtained for agronomic variables vs UAV data (Figures 5a,  5b, 5c) showed a similar trend for yield and pruning weight. In particular, they were characterized by H-H correlations with UAV data mainly within the high vigour area of the field whereas L-L correlations were present, as expected, in the low and medium vigour parts. Regarding the TSS variable vs UAV layers (Figure 5a), opposite associations prevail (L-H and H-L correlations) because, as shown in Table 4, this parameter has a negative association with UAV layers. These latter results clearly highlighted that both UAV and Sentinel-2 can correctly detect the vigour spatial variability within the vineyard. Through this study, therefore, it has been demonstrated that geostatistical analysis allows understanding the vineyard zoning to better characterize the areas with different qualitative and quantitative production, which are some of the key points for the terroir zoning. In fact, zoning enables the conservation, when necessary, of homogeneous areas or the differentiated management to create or recreate homogeneous zones lost due to unfavorable weather conditions or plant-soil interactions that have changed over time. The spatial analysis made by MI and BILISA methodologies showed the strong spatial heterogeneity overcoming the classical comparison methods by adding the spatial effect .

CONCLUSION
One of the focal points for vine growers is the need for specific agronomic management according to a proper knowledge of spatial variability and thereby reducing costs and environmental impact. In this study, the use of statistical methods that consider the spatial structure of vineyards to compare satellite and ground agronomic data with UAV data, has proven to be a robust and powerful tool to evaluate the practical implications of vineyard variability. MI and BILISA play a key role in identifying spatial patterns and could be successfully applied to other agricultural and environmental studies. Both Sentinel-2 images and UAV derived information were shown to be a relevant decision support for the stakeholders to manage vineyard vigour by defining withinfield management zones. In particular, good correlations were detected between S2 NDVI and UAV NDVI unfiltered. However, it is necessary to consider that the 10 m resolution of S2 cannot be used for any vine-specific application but only for larger scale purposes. Finally, it is necessary to focus attention on the current lack of studies concerning spatial statistical approaches to data analysis for agricultural purposes; it is on this point that research should concentrate its future developments.