From Space to Field: Combining Satellite, UAV and Agronomic Data in an Open-Source Methodology for the Validation of NDVI Maps in Precision Viticulture

Abstract

Recent GIS technologies are shaping the direction of Precision Agriculture and Viticulture. Sentinel-2 satellites and UAVs are key resources for multi-spectral analyses of vegetation. Despite being extensively adopted in numerous applications and scenarios, the pros and cons of both platforms are still debated. Researchers have currently investigated different aspects of these sources, mainly comparing different vegetation indexes and exploring potential relationships with agronomic variables. However, due to the costs and limitations of such an approach, a standardized methodology for agronomic purposes is still missing. This study aims to fill such a methodology gap by overcoming the potential flaws or shortages of previous works. To achieve this, an image acquisition campaign covering 6 months and over 17 hectares was carried out, followed by an NDVI comparison between Sentinel-2 and UAV to eventually explore relationships with agronomic variables. Comparative analyses were performed by using both classical (Ordinary Least Squares regression and Pearson Correlation) and spatial (Moran’s Index) statistical approaches: here, 90% of cases show r and MI scores above 0.6 for plain images, with these scores expectedly lowering to 72% and 52% when considering segmented images. Moreover, NDVI thematic maps were classified into clusters and validated by the Chi-squared test. Finally, the relationship and distribution of agronomic variables within NDVI and clustered maps were consistently validated through the ANOVA test. The proposed open-source pipeline allows to strengthen existing UAV and satellite applications in Precision Agriculture by integrating more agronomic variables.

1. Introduction

Precision Agriculture (PA) has been shaping modern agriculture in the last decades. PA is characterized by a wide number of techniques and technologies aimed at optimizing crop management and resources [1]. What made the rise of PA possible was the emergence of advanced modeling by Geographic Information Systems (GIS) and remote sensing technologies combined; they can provide affordable and reliable instruments, such as satellite platforms and devices for the acquisition of remote and field-based geospatial data [2]. In particular, Precision Viticulture (PV), the specialization of PA techniques for viticulture, started to develop significantly later than other crops (later in the middle of the 2000s decade), while experiencing rapid growth [3]. At present, the most recent technologies make it possible to carry out monitoring operations surveys through Unmanned Aerial Vehicles, UAVs. The availability of multi-spectral sensors on a commercial scale, ready to work on board these vehicles, makes it possible to perform spatial analyses by using different vegetation indexes at a sub-meter pixel size resolution. However, in the current fragmentation of the Italian wine-growing regions, being characterized by a mean size of about 1 ha [4], the use of UAV technologies can still be expensive for the agricultural business, both in terms of the monetary aspect and of the active effort required for the data collection [5]. Furthermore, the effects of climate change at regional and local scales make it more complex to collect historical data series that can be consistent over time and suitable to perform predictive analysis [6]. An alternative to vehicle-based monitoring in PV is represented by the European Program for Earth Observation represented by Copernicus. The dedicated satellites of this program among the Sentinel families, (Sentinel-2), active since 2015, are designed to collect multi-spectral images. This type of passive sensor brings numerous advantages for PV: the very first is the effortless collection of open big spatial data at a negligible monetary cost for the final users, able to perform territorial analyses over time. In fact, Sentinel-2 continuously scans the Earth’s surface, capturing multi-spectral images with a temporal frequency between 2 and 10 days depending on the latitudes, delivering post-produced products with a spatial resolution ranging from 10 to 30 m pixel size. Moreover, the Copernicus program makes historical satellite data available, allowing for both spatial and temporal analyses. The main downside of such a tool in PV applications is its low spatial resolution, especially when compared with the centimeter resolution provided by UAV platforms. Although more than sufficient for monitoring operations on a global scale [7], a decametric spatial resolution can represent a relevant issue for surveying operations on a scale such as that of single vineyards. Vineyards are indeed characterized by a regular canopy distribution but with the disadvantage of not being particularly dense and with inter-rows often covered by weeds or cover crops. This means that a single Sentinel-2 pixel can contain all these elements: vines, inter-row bare soil, weeds, or cover crops, with the latest representing a non-negligible source of noise when focusing on vine vigor assessment.

So far in the scientific literature, the issue has generally been investigated by comparing the Normalized Difference Vegetation Index (NDVI) values produced from both Satellite and UAV images with the aim of measuring their relationship. In fact, NDVI is the most widely adopted vegetation index in PV, and it can be useful for several purposes as a tool at the farmers’ disposal [8,9]. Multiple studies have investigated and analyzed this issue, with most of them focusing on different prominent aspects. This study takes into account a sample of these studies in order to elaborate a synthesis of their best practices and methodologies: Sozzi et al. (2020) considered a large area of about 42 hectares while collecting UAV data only once for spatial analyses comparison by focusing on the impact of border pixels and comparing data from both sources using Pearson Correlation and an Ordinary Least Squares regression (OLS) model [10]. With the condition of removing mixed pixels and edges, the authors concluded Sentinel-2 images to be comparable to UAV images as a source for broad canopy and vine management decisions within and between vineyard blocks, though currently unable to provide vine-specific information. Nonni et al. (2018) surveyed two blocks for a total of two hectares using classical OLS models for the comparison [11], noting how preliminary results show the potential for using open-source Sentinel platforms to monitor vineyards. However, in this case, the authors declared that further research was necessary in order to evaluate the possibility of using Sentinel-2 data for the long-term monitoring of vineyards. Also, Di Gennaro et al. (2019) preferred OLS methods for the comparison, while gathering insight about the relationship with agronomic variables such as biomass, yield, and grape composition sampled within different vigor zones at harvest time [12]. Here, the authors confirm the continuity of NDVI maps at a visual inspection but refine the analysis by filtering out ground pixels from UAV maps. The latter study also shows good correlation between NDVI maps elaborated from both platforms, especially considering filtered UAV data and agronomic key parameters. Khaliq et al. (2019) performed four surveys over 2.5 hectares, performing comparisons between maps using the Pearson Correlation before dividing them into different vigor clusters and testing their variance with the ANOVA test [13]. This latter study also considers vines-only maps to be most effective in describing the observed vineyard vigor. In fact, the NDVI maps derived from the satellite imagery were found not to be in accordance with the in-field crop vigor assessment. Instead, NDVI maps from UAV imagery, generated by considering vines-only pixels, were observed to be effective in describing the vigor of a vineyard. Pastonchi et al. (2020) collected images on a single 1.4-hectare block over three different years for a total of six surveys while tracking data about yield and pruning weight, putting emphasis on the need to move from classical statistical methods to spatially savvy techniques [14]. Matese et al. (2019) surveyed a single block of 7.5 hectares once per two years, using both OLS and spatial statistical methods for their comparison [15], underlining how the use of the latter provides a more robust and powerful tool to evaluate the implications of variability within vineyards. All covered studies tend to agree on the necessity to overcome classical statistical techniques in favor of spatial statistical methods since the latter proved to help prevent the overestimation of correlation when spatial patterns are not present and are considered to constitute a more robust evaluation tool. However, each study covered different portions of time and space, whereas studies that covered a larger number of hectares were being held for fewer surveys, and vice versa. Moreover, analysis of the relationship between agronomic variables and NDVI maps built from both aerial platforms was not always included. At last, most studies point out the need to expand the conducted analysis by covering both longer periods of time and larger crop extensions.

This study aims to harmonize what was conducted in previous works by taking advantage of the virtuous methodologies and good practices they employed and expanding the analysis on both a temporal (six surveys over six months) and a spatial (17 hectares) dimension. To achieve this, a multi spatio-temporal comparison between Satellite and UAV NDVI values has been carried out along with the integration of ground-truthed agronomic variables to further validate the correspondence between the vegetation index and actual geometrical measures taken on a number of monitored vines. To do so, after pre-processing both UAV and satellite images, NDVI maps were calculated using data from both platforms; UAV-derived maps were then filtered using a segmentation algorithm in order to remove ground-related pixels and to produce NDVI maps including only pixels related to the vines. UAV images, segmented and nonsegmented, were then co-registered to the corresponding Sentinel-2 images passing from a 4 cm/pixel to a 10 m/pixel resolution by averaging their original values. Once produced, the internal consistency of all sets of images was checked, introducing the Moran’s Index (MI) as a first spatial statistical method. Later, a comparison between Satellite and UAV NDVI values was performed by using more classical approaches and parameters, such as the $R^2$ of an OLS model and the Pearson Correlation coefficient (r) along with spatial correlation measures such as the Bivariate MI as suggested, for example, by Matese et al. (2019) [15]. NDVI maps from both sources were then classified into three vigor clusters based on their distribution terciles as seen in the works of Matese et al. (2019) and Di Gennaro et al. (2019), and their relationship was tested using the Chi-squared test [12,15]. At last, the relationship and distribution of agronomic variables with NDVI and clustered maps were analyzed through Analysis of Variance (ANOVA) [12]. In this study, since the correlation between qualitative grape parameters and NDVI values was deemed weak [12], two quantitative variables were monitored to verify their relationship with NDVI maps: the Average Length of Shoots and the Leaf Wall Area (LWA). These are two relevant parameters for wineries’ technical directors. The Average Length of Shoots is a parameter taken into account to plan for the beginning of phytosanitary treatments, while the LWA is directly related to the health status of vines and is often considered a relevant parameter for the dosage of phytosanitary treatments [16]. Both variables are not necessarily easy to monitor from flying platforms, hence the necessity of testing their relationship with NDVI values registered from above. The relationship and distribution of the values of the two agronomic variables within NDVI and clustered maps were consistently validated through the ANOVA test. Finally, the obtained results are compared with those shown in previous research, and the feasibility and interchangeability of both platforms for different applications in Precision Viticulture is discussed.

2. Materials and Methods

As stated in the introduction, this study aims to develop a standardized methodology for the multi spatio-temporal validation of NDVI maps. To achieve this goal, the proposed statistical analysis was built to include the best practices found in the covered previous research and synthesize them in the workflow depicted in Figure 1. The workflow is composed of the following steps, each of which is explored and detailed later within this section: the first block of the workflow is dedicated to data collection and pre-processing. Satellite images are downloaded from dedicated services, while UAV images and agronomic data are collected on the field. Satellite images are then cleaned from the presence of cloudy pixels, while UAV images are co-registered to match the resolution of the firsts. The second block of the workflow is dedicated to the calculation of NDVI maps and their division into clusters based on their tercile distribution. The final block comprehends all the statistical analyses carried out. NDVI maps are evaluated using both classical and spatial statistics; NDVI clusters are first tested using the Chi-squared test, and then the distribution of the agronomic variables over the clusters is evaluated with the ANOVA test.

2.1. Monitoring Sites, Plant Material and Weather Data

The area of study was that of the Consorzio di Suvereto e Val di Cornia Wine located in Alta Maremma, on the coast of Tuscany in central Italy. The Consorzio comprehends the municipalities of Suvereto, Campiglia Marittima, Piombino, San Vincenzo, Sassetta, and Monteverdi, as depicted in Figure 2. The area is particularly suited to the production of high-quality wines, given its strategic position close to the Tyrrhenian Sea, the minerality of the soils, and a non-negligible centuries-old tradition.

Out of more than thirty wineries, four were selected as study sites, with the intent to represent the diversity available in the area in terms of agronomic management, exposure, and macro composition of the soil, namely, Tenuta Casadei (Casadei), Società Agricola Petra (Petra), Azienda Vinicola Rigoli (Rigoli), and Tua Rita Società Semplice Agricola (Tua Rita) as shown in Figure 3, with their charachteristics summarized in

Table 1. Summary of the characteristics of the four wineries.

Winery	Soil Texture	Soil Management	Layout (m)	Cultivar	Extension (ha)	Training Method
Casadei	Silty clay loam	Alternated cover crop	1.80 × 0.80	Sangiovese	0.45	Double spurred cordon
				Merlot	1.21	Spurred Cordon
				Cabernet Sauvignon	1.11	Spurred Cordon
				Cabernet Franc	1.24	Double Spurred Cordon
Petra	Clay loam	Alternated cover crop	1.60 × 0.80	Sangiovese	/	/
				Merlot	2.82	Guyot
				Cabernet Sauvignon	3.39	Guyot
				Cabernet Franc	0.47	Spurred Cordon
Rigoli	Sandy clay loam	Spontaneous grassing	2.20 × 0.80	Sangiovese	0.88	Spurred Cordon
				Merlot	0.76	Spurred Cordon
				Cabernet Sauvignon	/	/
				Cabernet Franc	0.62	Guyot
Tua Rita	Loam	Alternated cover crop	1.40 × 0.80	Sangiovese	2.22	Double Spurred Cordon
				Merlot	1.49	Double Spurred Cordon
				Cabernet Sauvignon	0.41	Double Spurred Cordon
				Cabernet Franc	0.27	Double Spurred Cordon

.

To account for temporal variability, the presented analysis was repeated once per month during the vegetative season. Six surveys were planned to match UAV flights with the passage of Sentinel-2. The dates, shown in

Table 2. Dates of the six flights performed with UAV platforms.

Flight	Date
F1	24 March 2022
F2	28 April 2022
F3	23 May 2022
F4	27 June 2022
F5	22 July 2022
F6	26 August 2022

, were picked by checking the Sentinel-2 acquisition calendar at Spectator [17].

Grape varieties can represent another source of variability in terms of expressed vigor. To account for this, vineyards from the four wineries were also selected to be, as much as possible, equally distributed among the four most relevant Vitis vinifera grape varieties present in the area: Sangiovese, Merlot, Cabernet Sauvignon, and Cabernet Franc. In total, the study covered an area of about 17 hectares spatially distributed as displayed in Figure 3 and divided among different cultivars as shown in Table 1.

The 2022 season was characterized by one of the worst periods of drought in recent years [18], with practically no precipitation in the most critical months for grape fructification, during June and July (Figure 4). Relevant rainfall was finally observed during the month of August, when recorded precipitations were sufficient to safeguard the quality of the wine products. The exceptional drought, along with the high positive temperature anomalies registered during the summer, not only influenced the seasonal vegetative vigor, yield, and quality but also prevented the growth of cover crops and weeds in the vineyard alleys during the warmest months of the season (June, July and August), with a couple of exceptions arising during the second part of August as witnessed in

Table 3. Cover crops status at different wineries as found in different surveys.

	F1	F2	F3	F4	F5	F6
Casadei	On alternated rows	Just mowed	Absent	Absent	Absent	Absent
Petra	On alternated rows, uneven	On alternated rows, uneven	Absent	Absent	Absent	Partially present
Rigoli	Spontaneous	Just mowed	Just mowed	Absent	Absent	Partially present
Tua Rita	On alternated rows	Just mowed	Absent	Absent	Absent	Absent

.

Since weeds and cover crops are expected to have an impact on the statistical comparison of NDVI values, the four wineries were selected to account for such variability. At the Rigoli study area, weeds were left free to grow until late April before being mowed and controlled for the rest of the season. At the Casadei and Tua Rita study areas, cover crops were largely used (mostly Vicia faba minor L.) as a natural fertilizer and mowed by the end of April. At Petra, the study area presented a mixture of cover crops used on alternate rows. No other techniques to manage the vigor of the vines, and that could affect the recorded NDVI values, were used by the producers during the monitored period.

2.2. Satellite Images

Sentinel-2 products can be freely downloaded from the Copernicus Open Access Hub [20]. For this study, only Level-2A products were used, as they were available for all mentioned dates. Among all available products, only two spectral bands, band 4 (665 nm) and band 8 (842 nm), and the Scene Classification (SC) were considered. The two downloaded bands, respectively RED and NIR, were those needed to compute the NDVI spatial analysis according to Rouse et al. (1974) [8]:

$$NDVI=\frac{NIR-RED}{NIR+RED}$$

(1)

Once computed, Satellite NDVI output was cleaned using the SC layer as a mask. The SC comes at a resolution of 20 m spatial resolution, where pixel values range from 0 to 1 as in Figure 5. This product was extremely useful for masking out unwanted pixels and unnecessary noise from Sentinel-2 multi-spectral images. In this case, all pixels from the Satellite NDVI not corresponding to values 4 (Vegetation) and 5 (Not-vegetated) of the SC were removed. This process was performed mostly to avoid noisy pixels or pixels classified as clouds, which may loom over the studied vineyards, being included in the NDVI calculation.

2.3. Unmanned Aerial Vehicles Images

The UAV platform used for the present study was projected and realized by MAVTech Srl, a spin-off company from the Politecnico di Torino, with the collaboration of LIBRA SRL. This custom UAV platform, namely, the MAV-4QL-01/APR, was equipped with a Micasense Rededge M multi-spectral sensor characterized as described in

Table 4. Micasense Rededge M multi-spectral sensor specifics.

Focal distance	5.5 mm
GSD	8 cm/pixel (per band) at 120 m AGL
Field of View	47.2° HFOV
Resolution	1280 × 960
	B (475 nm ± 20 nm)
	G (560 nm ± 20 nm)
Spectral bands	R (668 nm ± 10 nm)
	Red Edge (717 nm ± 10 nm)
	NIR (840 nm ± 40 nm)

.

All surveys were conducted in the central hours of the day, between 10:00 and 14:00. To maintain a level of coherence within the collected data, the visiting order to the four wineries was always maintained the same for all three dates: Rigoli, Tua Rita, Casadei and finally Petra. Aerial images were acquired at a nominal resolution of 1280 × 960 pixels, from an average flight height of 60 m above the ground, proceeding at an average speed of 2 m/s and with an 80% forward and vertical overlap of the images. Since the drone had no RTK module, all images were georeferenced, exploiting the permanent elements of the landscape already present in the sites of the study as Ground Control Points (GCPs). For each site, at least 10 GCPs at the border of the vineyards and at least 5 GCPs within the perimeter of the vineyards were identified. This allowed for the coordinates of the GCPs to be recorded only during the first survey by making sure the selected GCPs were always going to be visible during the period of the study. The procedure was carried out with a differential GNSS system, namely, the Leica GPS1200, with an accuracy of 0.03 m. The collected pictures were then processed with photogrammetry techniques by the software Agisoft Metashape 1.8.3© to produce the required orthophotos. A camera geometric calibration procedure was performed before the image alignment task; moreover, a radiometric calibration was applied to the image blocks by using the reference images of the Micasense calibrated reflectance panel acquired before and after each UAV flight. Once multi-spectral images were collected, NDVI maps from UAV were calculated with Equation (1).

It must be noted that a calibration issue occurred at F1 at Tua Rita, at F2 at Casadei, and at F3 at Petra. All data corresponding to these surveys were removed from the study.

2.4. Images Alignment, Segmentation and Consistency Analysis

NDVI maps from both platforms were processed to be matched. In particular, UAV images needed to be re-sampled to match the Sentinel-2 images’ pixel size. To achieve this, Coregistration, an image-to-image automatic co-registration processing QGIS plugin, was employed [21,22]. UAV images were therefore brought from a 4 cm/pixel to a 10 m/pixel resolution by averaging the original values, and co-registered to the Sentinel-2 images. The vineyards’ polygons were then dissolved and given a 10 m internal buffer. This was performed since the inclusion of vineyard borders is well known to decrease the correlation between registered Sentinel and UAV data [10]. Though other studies have already applied masks for border effects, a statistical justification is still missing [23]. At this point, the previously computed NDVI maps were clipped using the buffered vineyards polygons, allowing maps derived from both platforms to perfectly match as shown in Figure 6.

As mentioned, vineyards are row-structured crops and, particularly at the early phenological stages when the LWA is still small as compared to the inter-alley distance, vine canopies are not the element of a vineyard that covers the highest percentage of pixels in an aerial image. Soil, shadows, weeds, and cover crops represent a prominent part of orthomosaics. In order to assess the relationship between satellite images and agronomic variables with pure-vine NDVI maps derivable from UAV images, a segmentation process was carried out to remove all non-vine pixels from images.

Correctly classifying UAV images to recognize vines, soil, and weeds/cover crops is one of the most important steps to ensure the high quality of the produced maps. For this purpose, many advanced techniques have been developed in the last few years, with Deep Learning algorithms being seemingly the most promising ones [24]. However, these techniques are still experimental and need more refining. At present, unsupervised algorithms have proved to be accessible and solid solutions for this task. Among the many available, for this study, the K-Means algorithm for RGB image segmentation, as described in Cinat et al. (2019) [25], was implemented in Python 3.10. Basically, the algorithm takes an RGB image of a vineyard as input and converts it into CIE-L*a*b* format to augment the performance of the K-Means algorithm. In the L*a*b* color space, colors are seen in a spatial distribution line of Luminance (L*), green-red line (a*), and yellow-blue distribution line (b*) [25]. The conversion from RGB to L*a*b* was performed with the skimage library, using the d65 reference white that simulates the not exposed colors. K-Means was then set to identify k = 5 clusters. The cluster with the lower value of a*, the strongest intensity of green, was identified as vines, while the shadows cluster was selected as the one with a positive value of a* and a lower value of b*. The remaining three clusters were identified as soil. The algorithm provided good results but resulted to be memory expensive. This was not a major issue for this work but it is still worth noting. Overall, the algorithm was quite capable of correctly classify vines, shadows, and soil as visible in Figure 7.

Operationally, to manage the memory workload, UAV RGB images were clipped by the vineyard polygons and fed to the algorithm. The returned segmented scenes were then remapped, giving a value of 1 to the vines cluster and a value of 0 to all other clusters, and merged back together. The mask built with this procedure allowed to remove non-vine pixels from UAV NDVI maps at their original resolution. At this point, the co-registration process described before was replicated, but only pixels corresponding to those belonging to the vines cluster were averaged to match the Sentinel resolution.

The described processing led to the production of three datasets, with each dataset composed of 24 images:

• NDVI S2: NDVI maps derived from Sentinel-2 multi-spectral images;
• NDVI UAV_plain: NDVI maps derived from full UAV multi-spectral images;
• NDVI UAV_vines: NDVI maps derived segmenting vines-only pixels from UAV_plain images.

Once produced, all images of the three datasets were tested for their internal consistency. To achieve this, the Moran’s Index (MI), a measure of spatial autocorrelation, was used [26]. MI is a measure of global spatial autocorrelation of the data within an analyzed area. Its value ranges from $-1$ to +1: for an observation at position i, a positive autocorrelation means that values in the neighborhood are similar to zi, while a negative spatial autocorrelation implies dissimilar values at nearby locations, and the zero value (MI = 0) indicates the absence of spatial autocorrelation. Univariate MI returns information about the correlation between a variable X and its spatial lag, formed by averaging all values of X for the neighboring polygons [26]. The analysis was processed using GeoDa 1.20.0.20 [27]. At first, the calculation of spatial weights was necessary since they are used to measure nearness and proximity between observations. GeoDa gives the possibility to manage the spatial weights.

As expected, all results shown in

Table 5. Univariate MI for the three datasets.

	Site	F1	F2	F3	F4	F5	F6
Univariate MI for NDVI S2	Casadei	0.848	/	0.769	0.802	0.857	0.857
	Petra	0.901	0.868	/	0.719	0.758	0.785
	Rigoli	0.662	0.851	0.934	0.882	0.887	0.867
	Tua Rita	/	0.743	0.753	0.673	0.855	0.863
Univariate MI for NDVI UAVplain	Casadei	0.624	/	0.718	0.859	0.774	0.756
	Petra	0.845	0.812	/	0.640	0.706	0.699
	Rigoli	0.626	0.700	0.856	0.759	0.852	0.825
	Tua Rita	/	0.673	0.785	0.718	0.736	0.433
Univariate MI for NDVI UAVvines	Casadei	0.695	/	0.566	0.831	0.777	0.856
	Petra	0.817	0.833	/	0.684	0.717	0.778
	Rigoli	0.793	0.898	0.883	0.545	0.960	0.675
	Tua Rita	/	0.682	0.637	0.823	0.726	0.717

were positive and show high MI scores. In fact, only two results scored a MI value lower than 0.6, indicating high spatial consistency for all produced images. p-values and z-scores were not reported, as they were all respectively <0.001 and >2.58.

2.5. Comparative Spatial Analyses

The classical approach in comparative studies is to use OLS models and analyze correlations measured as their coefficient of determination ($R^2$) values and Pearson Correlation between the interested variables. However, this kind of comparison has been deemed possibly biased, as common regression methods can only analyze the relationship between dependent and independent drivers but cannot consider spatial dependence [28,29]. To overcome such issues, a spatial statistical approach was integrated to provide higher reliability to the comparisons. In fact, Bivariate MI helps measure the spatial autocorrelation between a variable X and the spatial lag of a second variable Y, formed by averaging all values of Y for the neighboring polygons. For each flight campaign, NDVI maps from UAV images were compared to NDVI S2 values using either classical or spatial statistical approaches. The latter, unlike the classical method, considered the structure and geospatial variability present within the vineyard. All these comparisons were performed using the GeoDa software Version 1.22.

Continuous NDVI maps are indeed a valuable tool for wineries’ technical directors. These thematic maps can be further elaborated to produce more intelligible insights. Continuous values can, in fact, be divided into clusters to determine areas of Low, Medium, and High vigor as exemplified in Figure 8. These clusters, available in Appendix A, can give technical directors a valuable piece of information that can allow them more efficient vineyard management practices. In order to assess the quality of the information provided by such clustering, the values from all NDVI maps were separated on the terciles of their distributions, as seen in the works of Matese et al. (2019) and Di Gennaro et al. (2019), and tested using a Chi test, used to determine whether there is a significant association between two categorical variables in data [12,15]. Clustering by terciles allows to place the thresholds of the three clusters relative to the NDVI values of each parcel and to represent the internal variability of the latter.

2.6. In-Field Validation and Ground Data

To further evaluate Satellite and UAV data, 107 Regions of Interest (ROIs) were identified in order to monitor vigor-related parameters, such as the Average Length of Shoots and the LWA at each monitoring date during the season. Every ROI encompassed five vines representative of the area. At each plant, measures for both the considered agronomic variables were taken manually using a measuring tape. In particular, for the Average Length of Shoots, the average length of the four shoots more representative of the health state of the plant was recorded. The reported values were then averaged to obtain the Average Length of Shoots for each ROI. The LWA was calculated for each plant by multiplying the height of its leaf wall (hence, the previously obtained Average Length of Shoots) with the horizontal length of its leaf wall. Then, similarly to what was previously performed with the Average Length of Shoots, the reported values were then averaged to obtain the LWA for each ROI. The position of the ROIs was decided to respect the shape and variability of the vineyards as shown in Figure 3. Once determined, the spatial coordinates of all ROIs were recorded with the differential GNSS system.

To test the agronomic variables’ correspondence to what is visible through Satellite and UAV images, their distribution over NDVI clusters was tested using the ANOVA test. The ANOVA test, in this case, allowed to verify whether the observed distribution of the agronomic variables over the NDVI clusters created from all three datasets had a meaningful statistical difference [30].

3. Results

After detailing the context of the study and the proposed methodology, in this section, the obtained results are presented. In the first section, the results obtained from classical statistical approaches are shown in

Table 6. Coefficient of determination ($R^2$) between NDVI S2 and NDVI UAV_plain.

	F1	F2	F3	F4	F5	F6
Casadei	0.709	/	0.808	0.745	0.778	0.814
Petra	0.920	0.902	/	0.703	0.721	0.792
Rigoli	0.802	0.668	0.890	0.693	0.825	0.642
Tua Rita	/	0.338	0.589	0.693	0.825	0.642

F1–F6 represent the monitoring dates. Values are formatted according to the consistency of the relationship: $R^2$ ≤ 0.6 (normal); 0.6 ≤ $R^2$ < 0.7 (italic); 0.7 ≤ $R^2$ < 0.8 (bold-italic) and $R^2$ ≥ 0.8 (red bold-italic). The empty cells (/) means missing data. Pearson Correlation coefficient is not shown since all correlations were positive with values above 0.6.

and

Table 7. Coefficient of determination ($R^2$) and Pearson Correlation coefficient (r) between NDVI S2 and NDVI UAV_vines.

	F1		F2		F3		F4		F5		F6
	${\mathit{R}}^{\mathbf{2}}$	$\mathit{r}$	${\mathit{R}}^{\mathbf{2}}$	$\mathit{r}$	${\mathit{R}}^{\mathbf{2}}$	$\mathit{r}$	${\mathit{R}}^{\mathbf{2}}$	$\mathit{r}$	${\mathit{R}}^{\mathbf{2}}$	$\mathit{r}$	${\mathit{R}}^{\mathbf{2}}$	$\mathit{r}$
Casadei	/	/	/	/	0.435	0.660	0.655	0.809	0.715	0.845	0.801	0.895
Petra	/	/	0.527	0.726	/	/	0.118	0.343	0.041	0.203	0.287	0.536
Rigoli	/	/	0.625	0.790	0.787	0.887	0.256	0.506	0.469	0.684	0.678	0.823
TuaRita	/	/	0.369	0.607	0.152	0.389	0.396	0.629	0.674	0.821	0.763	0.873

F1–F6 represent the monitoring dates. Values are formatted according to the consistency of the relationship: r, $R^2$ ≤ 0.6 (normal); 0.6 ≤ r, $R^2$ < 0.7 (italic); 0.7 ≤ r, $R^2$ < 0.8 (bold-italic); and r, $R^2$ ≥ 0.8 (red bold-italic). The empty cells (/) mean missing data.

. Afterward, the results from the spatial statistical analysis, such as the Bivariate MI, described in Section 2.5, are presented in

Table 8. Bivaraite MI between lagged NDVI S2 and NDVI UAV_plain.

	F1	F2	F3	F4	F5	F6
Casadei	0.733	/	0.719	0.773	0.797	0.799
Petra	0.853	0.806	/	0.622	0.640	0.701
Rigoli	0.638	0.750	0.869	0.771	0.822	0.714
Tua Rita	/	0.348	0.679	0.616	0.776	0.439

Values are formatted according to the consistency of the relationship: MI score ≤ 0.6 (normal); 0.6 ≤ MI score < 0.7 (italic); 0.7 ≤ MI score < 0.8 (bold-italic) and MI score ≥ 0.8 (red bold-italic). The empty cells (/) means missing data.

and

Table 9. Bivariate MI between lagged NDVI S2 and NDVI UAV_vines.

	F2	F3	F4	F5	F6
Casadei	/	0.546	0.737	0.765	0.836
Petra	0.628	/	0.189	0.057	0.417
Rigoli	0.757	0.840	0.430	0.639	0.747
Tua Rita	0.411	0.342	0.509	0.691	0.765

Values are formatted according to the consistency of the relationship: MI score ≤ 0.6 (normal); 0.6 ≤ MI score < 0.7 (italic); 0.7 ≤ MI score < 0.8 (bold-italic) and MI score ≥ 0.8 (red bold-italic). The empty cells (/) mean missing data.

. At the end, the distribution of both agronomic variables over the NDVI clusters created from the three datasets is evaluated (

Table 10. Average Length of Shoots distribution on NDVI clusters based on NDVI S2.

	Casadei			Petra			Rigoli			Tua Rita
	Low	Medium	High	Low	Medium	High	Low	Medium	High	Low	Medium	High
F2	9.83	10.5	14.83	8.8	13.18	14.27	5.13	3.67	7.14	8	9.9	12.67
F3	47.86	59.57	63	74	80	86	42.5	53.29	62	54.83	53.5	61.25
F4	64.75	70.40	90.33	90.38	98.92	105	75	92.43	93.33	60.29	62.86	70.60
F5	69.86	68.40	71.75	97.50	98.17	98.56	77.50	93.86	96.67	57.89	62.17	65.67
F6	63.83	71.86	92	100.60	93.36	96	91	95.75	93.33	60.38	62	62.62

,

Table 11. Average Length of Shoots distribution on NDVI clusters based on NDVI UAV_plain.

	Casadei			Petra			Rigoli			Tua Rita
	Low	Medium	High	Low	Medium	High	Low	Medium	High	Low	Medium	High
F2	10.43	10.75	14.8	9.17	11.83	14.67	4.05	5.7	7.8	7	8.78	13.11
F3	47.86	57.83	64.40	74.38	84.29	84.09	42.4	50.57	62	52	54.75	64
F4	59.14	76.71	95	93.17	91.7	106.27	75.43	97.80	90	58.60	64.43	67.29
F5	62.4	72	75.5	108	90.89	97.82	79.29	94.4	96.25	58.11	63.67	64.44
F6	59.20	75.86	83.75	95.57	94.36	96.78	89.5	97.25	93.33	57.5	62.88	67

and

Table 12. Average Length of Shoots distribution on NDVI clusters based on NDVI UAV_vines.

	Casadei			Petra			Rigoli			Tua Rita
	Low	Medium	High	Low	Medium	High	Low	Medium	High	Low	Medium	High
F2	8.5	11.17	14.83	8.88	9.33	16.85	3	5	9.36	4.83	9.14	13.27
F3	47.86	57.71	66.25	64.5	80.23	97.14	42	51.75	62.67	45.5	56.25	64.63
F4	61.57	64.80	96.5	83.86	96.10	109.70	73.3	89.8	97.6	58.8	64.33	70
F5	56.14	79.83	82	94.67	97.22	102.56	74	94	95.60	56.67	64.14	65.75
F6	56.80	78.8	90	92.75	95.5	97.89	92.8	94.33	93.33	58.57	62.20	63.25

), before being deepened with the ANOVA test (Tables 16–18).

3.1. Classical Statistical Approach: Ordinary Least Square Regression

In Table 6, the coefficients of determination ($R^2$) obtained using a classical OLS model are reported. The model was implemented using NDVI S2 and NDVI UAV_plain as variables, the VIs computed from Sentinel-2 and full UAV imagery, respectively. It is possible to notice that about 90% of the resulting $R^2$ values were higher than 0.6 with the only low value being at F2 at Tua Rita, possibly due to the mowing of the cover crops having just taken place at the time of the survey. All Pearson Correlation values showed good correlations with values above 0.6.

The same procedure was applied considering NDVI S2 and NDVI UAV_vines; indeed, the NDVI maps derived segmenting pure vine pixels from UAV_plain images (Table 7). In this case, at F1, no values are reported since no foliage to be segmented was available before the budburst stage (24 March 2022). As expected, lower coefficient of determination values were obtained as compared to previous analysis (Table 6), with about 44% of registered $R^2$ values being higher than 0.6. The relationship between NDVI S2 and NDVI UAV_vines was expected to provide lower values, as under such discontinuous tree crop (vineyard), only a small part is represented by the vine canopy, and thus, in the mixed pixels of the NDVI S2, most of the effect is related to the background. Both Table 6 and Table 7 show values in line with what was highlighted by previous studies [10,11,12,14,15].

3.2. Spatial Statistical Analysis: Bivariate Moran’s Index

The spatial statistic approach was used to give a more unbiased insight into the correspondence of two spatially distributed variables by determining the presence or absence of spatial autocorrelation among them. The Bivariate MI analysis was performed comparing the NDVI maps of Sentinel-2 imagery with the NDVI computed on UAV_plain (Table 8) and UAV_vines (Table 9), respectively. In about 90% of the cases (flight dates and sites), the MI score resulted in being higher than 0.6, while the values resulted in being lower when segmenting UAV pure vines pixels (UAV_vines) where only 52% of maps showed $R^2$ values higher than 0.6 (Table 9). While these results could be apparently considered similar to those obtained using the classical OLS model, they are more representative of the inner relationships existing between S2 and UAV maps, as MI compare the values in a nearby matrix and not only cell by cell. As observed for OLS analysis, the relationship between NDVI S2 and NDVI UAV_vines (pure vine pixels) resulted in being coherently lower than the one performed on full UAV multi-spectral images (NDVI UAV_plain). Only two observations such as F4 and F5 at Petra site showed unexpectedly low values, yielding $R^2$ values of about 0.189 and 0.057, respectively. Even though at F4 and F5 cover crops were registered as absent from every vineyard in Table 3, at the Petra site, upon visual analysis there still were a few green patches that sent the segmentation algorithm astray, thus making behaviors such as this and that observed in Table 7 for the same flights possibly imputable to shortages in the segmentation process.

The comparison of NDVI clustered maps allowed to deepen the analysis between NDVI S2, NDVI UAV_plain and NDVI UAV_vines. All produced maps (provided as supplementary material in Appendix A (Figure A1, Figure A2, Figure A3, Figure A4, Figure A5 and Figure A6) showed appreciable results in terms of the tightness of the clusters. This was expected considering that in the results obtained using the Monovariate MI analysis, the NDVI values were deemed to be not randomly clustered. This notion was further strengthened by the use of the Chi-squared test [31]. The results obtained by the Chi-squared test performed between the datasets showed high significance (p-values < 0.001) in all cases, meaning that the categorical variables obtained (clusters) were correctly associated with each other, and, coherently, a spatial correlation between them existed.

3.3. Distribution Assessment of Agronomic Variables over NDVI Clusters

To assess for significant relationships resulting between NDVI clusters and ground-truth agronomic variables, indeed, the Average Length of Shoots and the LWA, different methods are required than those used to compare NDVI maps. To evaluate the accuracy of the agronomic variables’ correspondence with the three NDVI datasets, the values of both variables were further divided into the three NDVI clusters and averaged. For each variable and each dataset, it was possible to observe whether cells from NDVI maps classified as Low were showing lower values than those clustered as Medium and whether the latest were showing lower values than those clustered as High. This relationship was further quantified using the ANOVA, a statistical test used to analyze the difference between the means of more than two groups.

Table 13. ANOVA scores for Average Length of Shoots on NDVI S2.

	Casadei		Petra		Rigoli		Tua Rita
	F Statistic	p  Value	F Statistic	p  Value	F Statistic	p  Value	F Statistic	p  Value
F2	3.359	0.067	2.202	0.132	1.802	0.204	1.758	0.197
F3	3.854	0.045	0.637	0.538	5.982	0.014	0.455	0.640
F4	1.504	0.259	2.185	0.134	3.476	0.062	1.836	0.192
F5	0.023	0.978	0.008	0.992	4.707	0.029	2.664	0.093
F6	2.405	0.129	0.729	0.493	0.249	0.786	0.185	0.833

Values are formatted according to the consistency of the relationship: 0.01 < p-value ≤ 0.05 (bold-italic) and p-value ≤ 0.01 (red bold-italic).

,

Table 14. ANOVA scores for Average Length of Shoots on NDVI UAV_plain.

	Casadei		Petra		Rigoli		Tua Rita
	F Statistic	p  Value	F Statistic	p  Value	F Statistic	p  Value	F Statistic	p  Value
F2	2.169	0.154	3.029	0.067	2.634	0.11	4.427	0.025
F3	4.582	0.028	0.671	0.521	6.293	0.012	1.836	0.184
F4	2.915	0.09	3.849	0.035	5.465	0.019	1.179	0.333
F5	0.42	0.665	2.963	0.071	3.914	0.047	1.895	0.175
F6	2.155	0.156	0.102	0.903	0.739	0.508	4.193	0.029

Values are formatted according to the consistency of the relationship: 0.01 < p-value ≤ 0.05 (bold-italic) and p-value ≤ 0.01 (red bold-italic).

,

Table 15. ANOVA scores for Average Length of Shoots on NDVI UAV_vines.

	Casadei		Petra		Rigoli		Tua Rita
	F Statistic	p  Value	F Statistic	p  Value	F Statistic	p  Value	F Statistic	p  Value
F2	4.42	0.034	18.014	<0.001	25.897	<0.001	12.496	<0.001
F3	5.212	0.019	6.459	0.006	5.323	0.02	5.918	0.009
F4	5.259	0.021	13.433	<0.001	6.649	0.01	7.002	0.007
F5	3.159	0.076	0.632	0.54	8.056	0.005	4.432	0.025
F6	3.172	0.076	0.405	0.672	0.023	0.977	0.901	0.421

Values are formatted according to the consistency of the relationship: 0.01 < p-value ≤ 0.05 (bold-italic) and p-value ≤ 0.01 (red bold-italic).

and Tables 19–21 report the correspondence of the agronomic variables collected in the field with the NDVI clustered maps for the NDVI S2, NDVI UAV_plain and NDVI UAV_vines maps, respectively.

Table 10 reports the NDVI S2 dataset. It is possible to observe that in 80% of the cases, cells classified as Low NDVI well reported Average Length of Shoots values lower than those found in Medium NDVI ones. This percentage remained roughly consistent (85%) when considering NDVI UAV_plain (Table 11), while it is significantly increased to 100% of the cases when considering the NDVI UAV_vines (Table 12). In the cells classified as High NDVI, higher Average Length of Shoots values were registered as compared to those found in Medium NDVI on 95%, 85%, and 95% of the cases when considering NDVI S2 (Table 10), NDVI UAV_plain (Table 11), and 95% when dealing with NDVI UAV_vines (Table 12), respectively.

The ANOVA test increases the depth of the analysis around the distribution of agronomic variables over NDVI clusters. Table 13 shows how only in 15% of the cases, at least two out of three categories have significant differences in their means (p-value < 0.05) when considering NDVI S2, while this percentage increases to 35% in the case of NDVI UAV_plain (Table 14) and to 80% in case of NDVI UAV_vines (Table 15).

The same analysis was repeated for the LWA.

Table 16. LWA distribution on NDVI clusters based on NDVI S2.

	Casadei			Petra			Rigoli			Tua Rita
	Low	Medium	High	Low	Medium	High	Low	Medium	High	Low	Medium	High
F3	32.15	47.76	45.63	72.96	87.83	93.43	23.93	48.57	54.27	49.21	42.65	50.87
F4	59.94	59.03	98.66	92.67	106.72	107.16	58.27	87.74	90.18	49.05	59.22	70.31
F5	69.36	66.82	67.12	95.24	100.87	101.14	54.98	85.03	93.53	50.57	57.39	66.94
F6	63.32	72.23	79.39	96	96.19	96.06	85.58	98.17	96.28	50.19	61.73	66.38

shows how, considering NDVI S2, in 81% of the cases, cells with Low NDVI registered LWA values lower than those found in Medium NDVI. This percentage slightly lowers to 75% when considering NDVI UAV_plain (

Table 17. LWA distribution on NDVI clusters based on NDVI UAV_plain.

	Casadei			Petra			Rigoli			Tua Rita
	Low	Medium	High	Low	Medium	High	Low	Medium	High	Low	Medium	High
F3	32.15	46.51	46.51	74.84	98.18	89.02	23.56	43.54	55.96	46.14	39.94	54.92
F4	46.64	78.38	97.78	97.75	95.28	112.07	59.54	101.49	77.51	51.54	55.32	66.37
F5	66.85	69.01	67.68	103.13	97.34	99.47	59.99	83.94	91.52	48.18	59.85	67.69
F6	62.69	70.15	79.79	101.28	94.7	93.83	85.48	89.08	108.53	47.22	63.91	73.8

), while increasing to 94% when dealing with NDVI UAV_vines (

Table 18. LWA distribution on NDVI clusters based on NDVI UAV_vines.

	Casadei			Petra			Rigoli			Tua Rita
	Low	Medium	High	Low	Medium	High	Low	Medium	High	Low	Medium	High
F3	32.15	42.41	55.0	56.1	86.1	115.61	23.56	44.47	57.62	36.9	48.96	56.05
F4	47.13	63.2	106.2	83.22	101.71	117.25	54.85	84.08	96.96	33.81	62.6	70.61
F5	54.78	75.25	84.36	90.8	100.49	107.85	48.03	84.04	92.27	49.92	58.01	68.31
F6	56.77	76.06	79.2	88.64	97.14	101.62	88.46	85.31	108.53	43.25	58.58	69.22

). Moreover, cells with High NDVI register LWA values higher than those found in Medium NDVI in 81% of the cases when considering NDVI S2 (Table 16), 75% of the cases when considering NDVI UAV_plain (Table 17) and 100% of the cases when dealing with NDVI UAV_vines (Table 18).

For this variable, the ANOVA test shows slightly lower results, with 12.5% of the cases where at least two out of three categories have a significant difference in their means (p-value < 0.05) when considering NDVI S2 (

Table 19. ANOVA scores for LWA on NDVI S2.

	Casadei		Petra		Rigoli		Tua Rita
	F Statistic	p  Value	F Statistic	p  Value	F Statistic	p  Value	F Statistic	p  Value
F3	2.928	0.084	0.472	0.629	15.203	<0.001	0.568	0.575
F4	2.919	0.09	0.897	0.421	2.715	0.103	1.808	0.196
F5	0.015	0.985	0.146	0.865	5.059	0.024	3.045	0.069
F6	0.542	0.594	<0.001	1.0	0.726	0.513	2.266	0.129

Values are formatted according to the consistency of the relationship: 0.01 < p-value ≤ 0.05 (bold-italic) and p-value ≤ 0.01 (red bold-italic).

), while this percentage increases to 31% in the case of NDVI UAV_plain (

Table 20. ANOVA scores for LWA on NDVI UAV_plain.

	Casadei		Petra		Rigoli		Tua Rita
	F Statistic	p  Value	F Statistic	p  Value	F Statistic	p  Value	F Statistic	p  Value
F3	2.891	0.087	0.721	0.497	11.755	0.001	1.323	0.288
F4	5.874	0.015	1.36	0.276	5.117	0.023	0.926	0.417
F5	0.009	0.991	0.125	0.883	3.124	0.078	4.946	0.017
F6	0.63	0.548	0.206	0.815	3.045	0.104	8.997	0.002

Values are formatted according to the consistency of the relationship: 0.01 < p-value ≤ 0.05 (bold-italic) and p-value ≤ 0.01 (red bold-italic).

) and to 69% in the case of NDVI UAV_vines (

Table 21. ANOVA scores for LWA on NDVI UAV_vines.

	Casadei		Petra		Rigoli		Tua Rita
	F Statistic	p  Value	F Statistic	p  Value	F Statistic	p  Value	F Statistic	p  Value
F3	4.921	0.023	5.457	0.011	11.543	0.001	4.235	0.029
F4	19.247	<0.001	4.966	0.016	5.212	0.022	12.036	0.001
F5	2.009	0.174	1.368	0.274	9.516	0.003	3.797	0.039
F6	1.43	0.275	0.608	0.553	3.007	0.106	9.583	0.001

Values are formatted according to the consistency of the relationship: 0.01 < p-value ≤ 0.05 (bold-italic) and p-value ≤ 0.01 (red bold-italic).

).

4. Discussion

The obtained results help to better understand the relationship between S2 and UAV data and how they can be used to power applications able to monitor relevant agronomic variables. NDVI S2 and NDVI UAV_plain were indeed strongly related. In fact, when compared with classical methods such as OLS models, they showed high scores in both $R^2$ and r with the only outlier (at Tua Rita at F2) possibly being justified by the fact that just mowed cover crops could represent a non-negligible source of noise. The obtained results as seen in Table 6 show that in 90% of the cases, the value of r is higher than 0.60, resembling that seen in Khaliq et al., with such cases comprising 65% of the total [13].More generally, the correlation between NDVI S2 and UAV_plain, also explored in terms of the Pearson Correlation by Sozzi et al. and in terms of $R^2$ of an OLS model by Pistonchi et al., Di Gennaro et al., and Nonni et al., shows good and consistent values [10,11,14]. The values of the same coefficient were expectedly lower when comparing NDVI_plain and NDVI_vines. Khaliq et al. register a drop in values, where 75% of the r scores are lower than 0.41. This is in contrast to that seen in this study, where the number of cases found to have r scores higher than 0.6 drops but not excessively, from 90% to 72%. A higher drop is seen in $R^2$ though, going from 90% to 44% of cases, showing results higher than 0.6. The work of Pastonchi et al. also suggests a declining trend in OLS scores when comparing the correlation of NDVI S2 with NDVI_plain and NDVI_vines [14]. This seems to suggest that NDVI S2 should not be expected to be able to approximate NDVI values computed by looking at pure-vine pixels from UAV images. As already indicated though, classical statistical methods are to be taken lightly in a setting such as that considered in the present work. Stronger indications come by looking at more reliable measures such as the multivariate MI as already suggested by Matese et al. [15]. Table 8 and Table 9 indeed show how NDVI S2 and NDVI UAV_plain consistently produced better results than NDVI S2 and NDVI UAV_vines, confirming what was already glimpsed using classical statistical approaches. This complements what was reported by Pastonchi et al., where in some cases bivariate MI comparisons of NDVI S2 and NDVI UAV_plain show better results than comparisons of NDVI S2 and NDVI_vines [14]. This is indeed plausible but remains a minority among cases in terms of visibility when considering a higher number of cases as performed in this study.The NDVI value, though it is not necessarily valuable information by itself, gains value in its spatial relationship with a neighborhood of other values. The comparison of spatially clustered NDVI maps comes to help in this sense. In this case, in fact, the Chi-squared test returns significant results for all produced maps, indicating a relationship between NDVI clusters as created from all three datasets. As it was possible to expect, upon visual assessment, NDVI clusters from NDVI UAV_plain are extremely similar to those generated by grouping NDVI S2 values, while with NDVI UAV_vines, it is actually possible to spot a few relevant differences as seen, for example, at Petra at F4 (Figure A4) and F5 (Figure A5) (actually the two flights with the lowest results with both comparison methods). Even though there is a difference in NDVI S2 and NDVI UAV_vines, as testified by comparisons undergone with both OLS and multivariate MI, the first still managed to convey a relevant portion of the information carried by the second in terms of spatial distribution of NDVI data.

Despite the highest percentage of pixels in an image being represented by soil, NDVI values related to this element have very low variability, as they all physiologically tend to 0, thus enabling the distribution of vines vigor to come out, especially starting from F4 (27 June 2023) when the track of cover crops generally disappeared. Vice versa, at the time of the first three flights, cover crops were still present, even though with a decreasing presence, thus affecting the variation of values in NDVI S2. At the same time, the first flights represented a period of limited expressed vegetative vigor, hence, a period where monitoring a vineyard by the use of S2 images was predictably less feasible. As S2 NDVI clusters are considered to be similar to those producible by NDVI UAV_vines, they can represent a valid tool for monitoring vigor distribution in a vineyard through NDVI on broad scenarios and a valid source of information to support agronomic decisions. A conclusion shared by most other works present in the literature and taken into account by the present study [10,11,14], with the exception of Khaliq et al. that, even though from a smaller data sample, found NDVI S2 capabilities to be insufficient to monitor the spatial variability of vineyards’ vigor [13]. It is thus possible to state that NDVI S2 can be, with higher reliability in periods where cover crops are not present and foliage and vegetation are in advanced stages, a valuable and sufficiently reliable tool for broad applications.

On the other hand, UAV images could be used to complement this information, especially in the first stages of the season when cover crops and weeds might still be present and vine vegetation is at its initial and most delicate moment. For this purpose, it was relevant to analyze the ability of the two platforms to monitor the considered ground-truthed agronomic variables, the Average Length of Shoots and the LWA. Considering the first one, it was possible to observe how its values distribution among NDVI clusters was coherent in the majority of the cases for all three datasets (Table 10, Table 12 and Table 14). More interesting to observe were the results of the ANOVA test, which failed in most cases to denote significant differences between the means of the clusters when considering NDVI S2 (85%) and NDVI UAV_plain (65%) as shown in Table 13 and Table 14. This could be due to the fact that most of the NDVI UAV_plain information is related to soil and partially shadows, flattening the ability to show information about agronomic variables. In NDVI UAV_vines, instead, the selection of pure-vine pixels allowed a better assessment of cluster differences as seen in Table 15. Such conclusions resemble those found in other works, such as those of Khaliq et al. and Di Gennaro et al. [12,13]. The same goes if considering the LWA: its value distribution is mostly well separated within the three clusters for all three datasets (Table 16, Table 17 and Table 18), but their score values for the ANOVA test show a trend similar to that found for the Average Shoots Length. In fact, the ANOVA test scores are significant in determining differences between the means of the clusters when considering NDVI S2 only in 12,5% of cases (Table 19), while increasing to 31% (Table 20) and 69% (Table 21) of the times when considering UAV_plain and UAV_vines. NDVI_vines thus show to have a higher correlation than the other two datasets with agronomic variables. This reinforces the idea that satellite data are not suitable for vine-specific tasks and must be limited to broader monitoring activities, while segmented UAV images can help to gather NDVI values highly correlated with most agronomic variables.

Figure 9 helps to resume the information gathered during the discussion of the results of the statistical analysis and to conjugate the different possible usages of satellite and UAV data. NDVI maps from Sentinel-2 imagery proved to be a consistent and trustworthy tool for large-scale applications, with NDVI clusters such as those provided in the Appendix A (Figure A1, Figure A2, Figure A3, Figure A4, Figure A5 and Figure A6) are to be considered a valuable output for wineries’ technical directors. However, NDVI maps from satellite data suffer from some limitations which come to surface in the case that the monitoring is conducted during the vines’ early phenological stages or in the case that grassing or cover crops are present in the inter-rows of the monitored vineyard. UAV platforms have the ability to cover such shortages: despite having their own flaws, segmentation algorithms allow to separate inter-row pixels from UAV-derived NDVI maps and thus provide clean and reliable information. Moreover, given the existing correlation with quantitative agronomic variables, UAV platforms could provide insights into the early stages of the growing season. More in general, despite their usage coming at a greater cost, UAVs remain the preferred platform for performing vine-specific tasks.

5. Conclusions

Multi-spectral images from aerial platforms may represent an important source of information for the field of PA and PV. In the last decades, different studies have assessed the quality of information provided by satellite platforms with the intent to give agronomists and farmers valuable tools for decision support. Comparisons of S2 and UAV images were reported in several studies by highlighting the limitations and potentialities of both platforms. In this frame, the present study connected the most relevant points touched by previous works and conducted the analysis over a broad and consistent space and time span. The results obtained allow to stress the importance of bypassing the dichotomy between the usage of S2 and UAV and sustaining the necessity of their integration. This study shows evidence that S2 can be a valuable platform for large-scale monitoring, as it passively produces multi-spectral images at a tendentiously irrelevant cost, providing a strong correspondence with UAV imagery. In this case, the creation of NDVI clusters turned out to be a more useful output compared to the values of single cells. An NDVI value does not tell much about what happens in its 10 m/pixel cells, but it better explains if in its neighborhood similar and coherently clustered values are present. It is nonetheless relevant to understand how and how much this information can be disturbed by the presence of cover crops and weeds, especially when they are not uniformly distributed. These flaws can be efficiently bridged by UAV platforms that make the monitoring of grass-covered vineyards possible. Moreover, segmented UAV data show a good relationship with agronomic variables, thus becoming the reference for precision applications, especially in the initial stages of vine growth. Both these tools can take advantage of their integration by representing an invaluable source of information for agronomists and farmers in wineries and in other agricultural businesses.