Introduction

Precision viticulture is a management strategy that uses spatio-temporal data and observations to capture, describe and manage vineyard variability to improve productivity. The aim is to gain insights into how vineyard management can be improved through the adoption of new technologies to increase efficiency, quality, yield and sustainability, while minimising environmental impacts (Balafoutis et al., 2017a, 2017b). This is particularly important in regions where high-quality standards for wine production justify the use of site-specific management practises to increase grape quality and yield. Vineyards are inherently spatially variable and spatial production decisions in response to this variability can be very profitable (e.g. Bramley & Hamilton, 2004).

Canopy characterisation is an important component of site-specific crop management techniques. Remote sensing is widely used in agriculture, especially for monitoring crop growth and estimating crop quality and yield (Mulla, 2013). In viticulture, non-destructive remote sensing methods, especially geo-referenced information on canopy structure and field-level variability, are a highly beneficial alternative to in-situ field sampling, allowing rapid assessment of large areas (De Castro et al., 2018; Giovos et al., 2021; Hall et al., 2002; Johnson, 2003). Maps or images based on a variety of vegetation indices (VI) from remote sensing can be used to monitor the condition of vegetation and quantitatively assess its health and vigour. The spectral response of the vine canopy needed to calculate VIs can be determined directly, accurately and non-destructively either remotely or proximally (at close range) using a variety of sensors and sensor configurations on terrestrial, airborne or satellite platforms (Hall et al., 2002; Sozzi et al., 2020). There are many different VIs that can be calculated based on the available spectral (band) information collected, but the most widely used in viticulture (and crops generally) is normalised difference vegetation index (NDVI) (Rouse et al., 1973; Taylor & Bates, 2021). The NDVI uses a combination of red and near infra-red reflectance, which relates to vegetation vigour and biomass. The NDVI has been extensively researched and associated with various structural and physiological characteristics of grapevines and is often used as data for spatial management decisions in vineyards (Acevedo-Opazo et al., 2008). Recent work has highlighted that the NDVI may not be the most effective VI in viticulture (Matese & Di Gennaro, 2021; Taylor et al., 2021), but it is still a relatively robust VI and remains the default VI used in many research and commercial applications.

Variations in vineyard micro-climate and soil properties lead to significant changes in the size and structure of the vine canopy within vineyard blocks (Gatti et al., 2022; Hall et al., 2011). NDVI data can provide information on vine canopy variation that in turn may affect the critical characteristics of yield and berry composition (Hall et al., 2002). A common approach to relate NDVI to vineyard production is to perform statistical and regression analyses, including descriptive statistics, Pearson correlation and regression models. Pearson’s correlation coefficient has been quantified in several studies to determine the spatial correlation between VIs (usually but not always NDVI) and yield components, suggesting that the zones with the highest vine vigour are also the highest yielding (e.g. Anastasiou et al., 2018; Gatti et al., 2016; Hall et al., 2011; Sun et al., 2017). Others have investigated correlations between NDVI and fruit quality indices, such as berry size, total soluble solids, anthocyanins, grape colour and total phenolic content (Anastasiou et al., 2018; Baluja et al., 2012; Hall et al., 2011; Kasimati et al., 2021a). These studies differ in their relationships between NDVI, yield and quality, often because of differences in the viticulture system. Although higher NDVI values are generally associated with higher yields and lower quality, this is not always the case. The reverse relationship has been observed, with zones of low vigour (NDVI) tending to produce higher quality wine in some instances.

In addition to direct comparisons between VIs, particularly NDVI, and crop production attributes, other research has examined the strength of relationships between different types of NDVI maps in the same vineyards, e.g. Sozzi et al. (2019) and Kasimati et al. (2021b). This has shown that there are differences in the spatial structuring of VIs when the images come from different types of sensors. However, whether comparing NDVI maps from different sensors or comparing NDVI values to vineyard production, a commonality is that these relationships tend to only be reported in a statistical context and mostly via correlation analyses. None of the previous papers have considered how these relationships would be translated into decisions in an operational context, i.e. very little was understood about how differences in NDVI maps might influence spatial production decisions and ultimately affect the management of yield and quality variability in vineyards.

An alternative approach to correlation and regression analyses is to use classification algorithms. Spectral information can be classified, with or without other data such as soil or production information, to ‘group’ these data into potential management units (MUs) or zones (MZs). To this end, clustering techniques have been widely used in precision agriculture and precision viticulture, in particular the k-means algorithm and its variants (Fridgen et al., 2004; Taylor et al., 2007). The k-means algorithm assumes that these different data sources are located in one place. This often requires appropriate interpolation as a pre-processing step, as data collected at different times and/or with different measurement systems are rarely collected at the same location. However, the advantage of this method is that classes (MUs) are generated that have a more homogenous (i.e. VI) response within each MU, and these can be considered as operational management units.

In previous work, NDVI data from four sources, including two terrestrial crop reflectance sensors, a UAV-mounted multispectral camera and Sentinel-2 satellite imagery, were analysed using both a correlation and a MU approach to determine how similar these data sources were in a production environment (Kasimati et al., 2021b). The correlation coefficients and MU maps showed that NDVI data from different sensors, especially from satellite platforms, provided different information for decision-making. When it comes to production considerations, not all NDVI maps are equal. However, this previous research only investigated the temporal evolution of operational relationships between NDVI maps in vineyards over time. It did not examine the operational relationships of different NDVI data with grape yield and quality within the same vineyard.

Therefore, building on previous work, this research compared NDVI data collected by different proximal and remote sensing systems over the same vineyard at the same time and determined the extent to which these data sources correlate with each other and with grape yield and quality data at the operational level. The operational, production context was explored by including management units in the analysis, which currently represent the most common type of spatial and variable management in crop production systems. The aim was to determine whether different NDVI data sources at different times in the season, which respond differently to some extent, ultimately lead to different management decisions.

Materials and methods

Study area

The study was conducted in a commercial vineyard located at 420 m above sea level on a steep slope (16%) on a commercial vineyard, the Palivos Estate in Nemea, Greece (37.8032° N, 22.69412° E, WGS84). It is planted with Vitis vinifera L. cv. 'Agiorgitiko' for wine production, and the training system is a vertical Guyot system with cane pruning. The data used for the analysis were obtained in 2019 and 2020 from an area of approximately 2 ha within the estate. The row orientation is northeast-southwest and the row spacing is 2.2 m (Fig. 1).

Data collection—sensors

A HiPer V RTK GPS (Topcon Positioning Systems Inc., Livermore, CA, USA) was used to record the field boundaries and elevation data. Four surveys were performed during each of the two seasons (2019 and 2020) based on identified key phenological stages and the quality of sensor data obtained in each campaign (Kasimati et al., 2021b). The four dates selected were early June, early July, early August and late August, representing floraison, fruit set, berries pea-size and veraison stages (Table 1).

Table 1 Canopy reflectance, proximal and remote sensing data, acquisition dates and the corresponding phenological development stages of grapevine by Eichhorn and Lorenz (1977), referred to as the EL scale
Full size table

For each campaign, vine vigour was assessed using two vehicle-mounted proximal sensors simultaneously using an ad-hoc mounting structure (Fig. 2). These sensors were a CropCircle active proximal canopy sensor (ACS-470, Holland Scientific Inc., Lincoln, NE, USA) and a SpectroSense + GPS passive canopy sensor (Skye Instruments Ltd, Landrindod Wells, UK) were mounted on a tractor. The CropCircle was mounted horizontally to image the side curtain of the canopy, following the method of Drissi et al. (2009) while the SpectroSense + GPS was mounted above the canopy to collect data from a nadir position. Both sensors were adjusted to maintain an appropriate distance (0.3–0.5 m), according to manufacturer specifications, from the vines at each growth stage. They record proximal reflectance measurements from the side and top of the canopy. All proximally recorded data were georeferenced using either a Garmin GPS16X HVS (Garmin, Olathe, Kansas USA) or the integrated built-in global navigation satellite system (GNSS).

Remote assessments of vine vigour were made using Sentinel-2 satellite imagery and a Phantom 4 Pro UAV (SZ DJI Technology Co., Ltd., Shenzhen, Guangdong, China) equipped with a Parrot Sequoia + multispectral camera (Parrot SA, Paris, France) and associated 3-axis georeferencing metadata using the camera’s integrated positioning system. The aerial UAV imagery were acquired around midday with nadir flights at 30 m above the ground on the same days as the proximal measurements. The acquisition interval of the multispectral camera was set at 2 s, and the flight plan overlap and sidelap were 80% and 70%, respectively. The ground sampling distance (GSD) of the image ortho-mosaics was ~ 30 mm. Atmospherically corrected Sentinel-2 satellite 2A images, with a spatial resolution of 10 m pixels and a 5-day revisit time, were downloaded from the official Copernicus Open Access Hub (scihub.copernicus.eu) for the date closest to the date of the proximal and UAV surveys. This usually occurred within 2 days in the mid-late season surveys but was up to 9 days after the ground-based observations earlier in the season due to high cloud coverage effects in the nearest ingestion dates (Table 1).

Two Vantage Pro 2 weather stations (Davis Instruments Corp., Hayward, CA, USA) (Fig. 2, Left) with a rain sensor to detect precipitation, anemometer to measure wind speed and direction, air temperature sensor, humidity sensor and barometer to monitor air pressure were installed inside the vineyard. The weather data were recorded during the entire growing season.

Data collection—production attributes

The quantity and quality of the harvest was sampled manually, with grapes harvested by hand at the end of each growing season (mid-September). A regular 100-cell grid (10 m × 20 m) covering the entire area (Fig. 1) was created using ArcMap v10.3 (ESRI, Redlands, CA, USA) and then laid out in the field by creating landmarks to facilitate field measurements of grape yield and quality. This sampling strategy was agreed with the vineyard manager and the grid size was considered by the vineyard manager to be the basic minimum size for management. Yield was determined by counting the total number of crates filled per cell and multiplying it by the average crate weight of the harvested grapes (Fig. 3). Grape quality characteristics were assessed by randomly selecting fifty berries from each cell in the vineyard. Qualitative analysis of the common vineyard quality indicators, total soluble solids (°Brix) in must, total titratable acidity and pH, was carried out according to Stavrakaki et al. (2018) at the Laboratory of Viticulture, Agricultural University of Athens. Total soluble solids in the must were measured without temperature correction using an ATAGO N1 refractometer (Atago Co. Ltd., Tokyo, Japan) with a measuring range of 0–32 °Brix in 0.28°Brix increments. For the determination of the titratable total acidity, a titration with a 0.1 N NaOH solution was carried out. The most abundant organic acid in grapes, tartaric acid, was used to represent the total titratable acid. The must pH was determined using a pH meter.

Fig. 1
figure 1

Left: the national (inset) and local position of the study vineyard, with the target area within the Palivos Estate indicated by the white rectangle and; Right: close up of the satellite image of the commercial vineyard and the 100-cell grid developed parallel to the trellis lines (Google Earth Pro, 2021)

Full size image
Fig. 2
figure 2

Pictures of vineyard-based data collection with different sensors during the survey. Left: two vehicle-mounted sensors, SpectroSense + GPS and CropCircle, for measuring canopy reflectance; and Right: The Phantom 4 Pro UAV in use within the vineyard to collect proximal imagery. An on-site weather station is also shown in the left image

Full size image
Fig. 3
figure 3

Pictures taken during the collection of manual grape production data. Left: the regular 100-cell grid (10 m × 20 m) plan covering the entire area and empty baskets ready to be filled; and Right: filled baskets in the vineyard awaiting collection and recording

Full size image

Data preparation

The geographic co-ordinates for all the tractor-mounted, proximally acquired canopy reflectance data collected during the season were first converted to projected co-ordinates (UTM Zone 34 N) and pre-processed by cleaning and removing data points outside the field boundaries and removing outliers (values >  ± 2.5σ) (Taylor et al., 2007). Data were interpolated to a common 2.2-m grid (corresponding to the row spacing in the vineyard) using block kriging on a 5 × 5-m block size with a local variogram with the Vesper software (Minasny et al., 2005). Using ArcMap, the interpolated data were scaled up to 10 m × 20 m plots (corresponding to the cells used for manual sampling) using a mean aggregation of the 2.2 m grid points. Ultimately, this produced a time series of NDVI maps with a spatial resolution of 10 m × 20 m, aligned parallel to the vineyard rows.

The multispectral images captured by the UAV were mosaicked using Pix4D software (Pix4D S.A., Prilly, Switzerland) through the ‘Ag Multispectral’ photogrammetric model pipeline. Radiometric calibration was applied to the generated ortho-mosaic using the reference images of a radiometric calibration target (airinov aircalib) acquired after each flight. Finally, the generated NDVI ortho-mosaic was fitted to the boundaries of the vineyard and the data were scaled up to the same 100 management plots (10 m × 20 m) using a mean aggregation approach.

The Sentinel 2A images have a spatial resolution of 10 m and the NDVI was calculated using bands 4 and 8 (red and near-infrared, respectively). A spatial correction was made along the boundaries of the experimental field based on ground control points from the detailed UAV map before scaling up the data to the 10 m × 20 m plots by averaging the NDVI of all pixel centroids within the plots. A schematic and examples of the collection and processing flow for each data source is shown in Fig. 4.

Fig. 4
figure 4

(Reproduced from Kasimati et al., 2021b)

Sequence of data processing and created NDVI maps. The first row shows the sensor type and positioning relative to the vine canopy (orientation and relative distance) during data acquisition. NDVI images derived from either interpolated ground points or captured images (second row) at their original resolution were used to create a series of transformed NDVI maps with a spatial resolution of 10 m × 20 m, aligned parallel to the vineyard rows (third row). This was done at four dates in both years to track the spatio-temporal change in the vine canopy response.

Full size image

The final step of data processing was the creation of a geodatabase at the grid cell level that contained the centroid co-ordinates of each cell (management plot), the mean NDVI values of each sensor for a given cell at a given date, the date associated with each sensor survey and the manually acquired yield and quality values for each cell.

Data analysis

Exploratory analysis of NDVI and crop production attributes

Descriptive information relating to the NDVI data used in this study has previously been published (Kasimati et al., 2021b). Pearson’s correlation analysis using R software (R Core Team, 2022) was carried out to investigate the relationships between the NDVI data from all four sensors for the four selected measurement dates to examine whether the different combinations of NDVI maps exhibited a homogeneous trend. Correlation r-values > 0.50 were considered to exhibit a moderate to strong relationship.

Descriptive statistics (mean and standard deviation) for the crop production data were generated to quantify the yield and quality attributes. Pearson’s correlation was again used to directly link the mean NDVI response to crop production attributes at the 100 m grid cell level as a preliminary analysis.

Delineating management units (MUs)

Before delineating the MUs, PCA was performed to analyse the relationships between NDVI maps. The PCA was performed using the built-in R software program function prcomp() (R Core Team, 2022), to group the NDVI data into statistical factors and create linear independent variables that remove multi-collinearity and describe the spatio-temporal information collected. The results of the PCA (principal components—PCs) were plotted to show trends and also saved to be used in the delineation of the potential MUs.

Fuzzy k-means classification using the Management Zone Analyst (MZA) software (Fridgen et al., 2004) was used to classify the data into MUs. The MZA software was operated with the following parameter settings; using the Mahalanobis similarity measure for multivariate clustering, fuzziness exponent = 1.3, the maximum number of iterations = 300, and the convergence criterion = 0.0001. The FPI (Fuzziness Performance Index) and NCE (Normalised Classification Entropy) indices were used to determine the optimal number of zones (Fridgen et al., 2004; Tagarakis et al., 2013) by minimising both indices. Initially, clustering was performed with all selected dates and with all sensor data. This procedure was referred to as the ‘All’ approach and was considered as the ‘reference’ MU map. Then, each individual NDVI sensor time series for both years were classified separately to create sensor-specific MU maps. A MU map was also created by classifying the PCA outputs. The results of the ‘All’ classification in each year were compared with the map results of each individual sensor or PCA. Comparisons were made by determining the Degree of Agreement (DoA) (Tagarakis et al., 2013) between MU maps, i.e. the number of cells or ‘pixels’ belonging to the same class in the reference (All) and target maps. This method provides a visualisation of the spatial NDVI correlation in terms of production, but does not provide a direct linear correlation index.

Analysis of variance (ANOVA) was used to test how well the variability in grape yield and quality is explained by the various sensor-specific MU maps. Each harvest variable (yield, TSS, TA and °Brix) was analysed individually against the sensor-specific MU maps using a one-way ANOVA analysis. As the MUs formed the defacto ‘treatments’ for the ANOVA, and were unequal in size, an unbalanced one-way ANOVA was performed and the adjusted R2 value recorded. The threshold of significance (p-value) for the statistical tests was 0.05. The analysis was performed in R with the package MANOVA.RM (Friedrich et al., 2019).

Results

Weather data

The climate of Nemea is described as Mediterranean, mild, with an average (thirty-year) rainfall of 690 mm distributed from October to April. Total precipitation ranged from 1093 mm (for 2019) to 518 mm (for 2020) and the average annual temperature ranged from 17.0 °C in 2019 to 18.1 °C in 2020, with a maximum temperature of 40.3 °C in 2020, while the lowest temperature was − 7.6 °C.

Descriptive analysis of yield and quality parameters

The mean (μ), standard deviation (σ) and range of the yield and quality attributes across the 100 sampling plots is given in Table 2. Mean yields (and variances) were low (50 kg/plot = 2.5 Mg/ha), but typical for premium wine production of this region. Yields in 2019 were higher than in 2020. The Brix values were similar in both years and these values according to the estate winemaker are indicative of “ripe” and ready to harvest grapes. The pH values were considered high by the winemaker, especially in 2019, and are associated with over-ripe grapes. In contrast, the TTA values were considered normal values for Vitis vinifera L. cv. ‘Agiorgitiko’, although much higher TTA values were observed in 2019 than in 2020.

Table 2 The minimum and maximum values along with the mean and standard deviation of all yield and quality parameters
Full size table

Exploratory correlation analysis between NDVI layers

As reported previously (Kasimati et al., 2021b), there were generally good positive correlations (r > 0.50) between the NDVI maps from different sensors and dates in both 2019 and 2020, indicating high similarity between the different sources of information. For 2019, the correlation coefficients between sensors per date were strong (r > 0.70), indicating sensor reliability. The maximum correlation (r = 0.90) for 2019 was observed for the Sentinel-2 data between early and late August (i.e. in the mid-late season with full canopy growth). For 2020, the correlations were strongest (r = 0.96) for the UAV data, between early June and August. All Sentinel 2 images were strongly correlated (0.69 < r < 0.96); however, in 2020, the satellite NDVI had relatively weak correlations with the other three sensing systems (0.27 < r < 0.40). Among the sensors for specific dates (phenological stages), the SpectroSense + GPS was found to correlate best with the other sensors in both years. However, no temporal pattern was observed as the best correlations were noticed in late season in 2019 and early season in 2020. (Full details on these correlation results can be found in the supplementary material, Table S1).

Exploratory correlation analysis between NDVI layers and grape yield and quality

Summary information relating to the exploratory correlation analysis between sensor-specific NDVI layers and grape yield and quality are shown in Tables 3 and 4 respectively. The upper part of Table 3 shows the best intra-sensor correlations between NDVI data and yield, attaining higher values at the middle of the growing season. The maximum correlation (r = 0.91) for 2019 was observed in early August for UAV data, followed by early and late August for the Sentinel-2 and UAV data, respectively (r = 0.90). In 2020, correlations were stronger in early August for the UAV data (r = 0.81), followed by early July for the SpectroSense + GPS data (r = 0.79). All Sentinel-2 NDVI variables had weak correlations with yield in 2020 (0.32 < r < 0.35). For the given dates during the growing season, the UAV NDVI data was found to correlate best with yield in early and late August.

Table 3 NDVI–Yield Correlation Co-efficient: The first two higher Pearson’s correlation coefficients between NDVI data from all four proximal and remote sensors and yield for a given sensor (upper part) and for given time periods during the season (bottom part)
Full size table
Table 4 NDVI–Total Soluble Solids Correlation Co-efficient: Selected best performed Pearson’s correlation coefficients between NDVI data from all four proximal and remote sensors and total soluble solids for a given sensor (upper part) and for given time periods during the season (bottom part)
Full size table

Table 4 shows the equivalent correlations for TSS instead of yield. The correlation matrices for 2019 and 2020 generally showed a low negative correlation with TSS, with slightly lower values for 2020 than for 2019. The best correlation (|r|= 0.74) was observed for SpectroSense + GPS data at the end of August 2019, followed by CropCircle NDVI in late August (|r|= 0.69). For 2020, correlations were strongest for SpectroSense + GPS data (|r|= 0.70) in early July and late August. All Sentinel-2 NDVI variables had relatively weaker correlations with TSS in both years (0.28 <|r|< 0.58) when correlated with total soluble solids. Across both years, the higher resolution NDVI data (SpectroSense + GPS, CropCircle and UAV) performed consistently well, with r > 0.5 on at least one date.

The pH and TTA data had low correlation values (|r|< 0.30) with the NDVI data, regardless of the sensing platform, at all crop stages. This indicated that pH and TTA was not determined by vine vigour in this vineyard system.

Application of principal component analysis (PCA) to the NDVI data layers

PCA was applied to all the NDVI maps that constituted each year’s dataset (four sensors at four dates). For the 2019 dataset, the first two principal components (PC1, PC2) described 85.7% of the variability in all NDVI maps. For the 2020 dataset, PC1 and PC2 described 85.2% of the variance in the dataset. The results of these three PCAs are summarised in Table 5 and plots of the first two PCs are shown in Fig. 5.

Table 5 Cumulative proportions and standard deviations associated with principal components (PCs) from the analysis of all selected NDVI layers from all sensors in each year (2019 or 2020)
Full size table
Fig. 5
figure 5

PCA cluster plots visualising the values of PC1 vs. PC2 and the relative statistical ellipses for the three classes (95% confidence interval). The closer the major axes of the statistical ellipse of each class are to the parallel, the higher the reliability of the classification. (Class 1/Red: low NDVI = 0.48–0.56; Class 2/Green: medium NDVI = 0.53–0.61; Class 3/Blue: high NDVI = 0.56–0.63 as determined by the MZA software analysis (see next section))

Full size image

Delineating management units from combinations of NDVI layers and PCs

In the fuzzy k-means classification, the FPI and NCE indices were most often minimised when k = 3, so this value was adopted for all classifications to ensure direct comparisons were possible. The fuzzy k-means classification and the MU maps as well as the total DoA analysis between the different sensors are shown in Fig. 6 and Table 7. The classes were labelled according to the mean NDVI response, with the red cluster always representing the lowest mean NDVI value, the green cluster representing intermediate mean NDVI values and the yellow cluster representing the highest mean NDVI value. Table 6 presents the NDVI ranges of each sensor-specific classification for the 3-class reference classification (low, medium and high NDVI) for both years.

Fig. 6
figure 6

(Reproduced from Kasimati et al., 2021b)

Management units (3 classes) resulting from the fuzzy k-means classification of different NDVI values obtained from different sensors in 2019 and in 2020. The analysis was carried out on 10 m × 20 m plots, which were considered as minimum management plots. Each plot was assigned to a class based on its mean NDVI response (Red: low NDVI; Green: medium NDVI; Yellow: high NDVI)

Full size image
Table 6 The NDVI ranges of each sensor-specific classification for the 3-class reference classification (low, medium and high NDVI) for both years
Full size table
Table 7 The Degree of Agreement (DoA) table shows the number of plots (10 m × 20 m) that agreed between the 3-class reference classification (All) (low, medium and high NDVI) and the corresponding 3-class output of each sensor-specific classification for 2019 and 2020 across all 100 plots
Full size table

The resulting MU maps (Fig. 6) showed the highest DoA between the different approaches for the high NDVI class, followed by the low and then the medium NDVI classes, i.e. regardless of the sensing platform chosen, the NDVI information was able to consistently represent high vigour zones but was less consistent on the zoning of intermediate or medium NDVI values. The SpectroSense + GPS data showed the highest DoA with the reference map (‘All’) in 2019 (85%), while it was the Sentinel-2 data in 2020 (71%). This latter result is due to the large deviation between the Sentinel-2 maps and the other sensors, such that the All (and PCA) result has a skew within it to account for the satellite data deviation, and subsequently fits the satellite data better. The MZA-derived ‘All’ clustering was almost identical to the PCA-based MU map in both years, which was to be expected as PCs were derived from the ‘All’ dataset.

ANOVA

Each harvest quality variable (TSS, TA and pH) as well as yield was analysed individually from sensor-specific MU maps using a one-way ANOVA analysis. The pseudo-design for the ANOVA was unbalanced because there was a difference in the number of plots assigned to each class based on their average NDVI response. The adjusted R2 for the fits of the ANOVA are shown in Table 8. These values are indicative of the amount of variability in yield and quality attributes that was explained by the MUs created from the reference (All), PCA and sensor-specific NDVI data.

Table 8 Adjusted R2 from ANOVA for each grape production attribute: grape quantity (yield) and quality (total soluble solids in °Brix, pH and total titratable acidity)
Full size table

As for the correlation analysis, the NDVI-derived MUs had good fits with yield and TSS, but poorly explained the variance observed in TA and pH. Yield for both years was better explained by UAV-generated MUs, with an adjusted R2 of 0.76 and 0.59 for 2019 and 2020, respectively. The TSS data were better explained by PCA-derived MUs in 2019 (adj. R2 = 0.49), followed by SpectroSense + GPS and All MUs. The TSS ANOVA results were poor in 2020 with only the UAV-derived NDVI MUs showing any relationship with TSS (adj. R2 = 0.30).

The mean production responses for the MUs are provided in Supplementary Data (Table S2). The low and medium NDVI MUs were characterised by similar yields, that were lower than the high NDVI MU. The medium and high NDVI MUs were characterised by similar but lower °Brix values than the low NDVI MU. There was little differentiation in pH and TA between the MUs, which is consistent with the observed correlation results with NDVI. Therefore, the three MUs can be relatively considered as (i) Low NDVI MU with low yield and high TSS; (ii) Medium NDVI MU with low yield and low TSS; and (iii) High NDVI MU with high yield and low TSS.

Discussion

Comparison among the NDVI data from all four sensors, i.e. the vine-only NDVI values extracted from two proximal sensors (CropCircle and SpectroSense + GPS) and two remote sensors (UAV and Sentinel-2 images) showed increasing correlations as the season progressed and the canopy size expanded. This conforms with previous studies (Hall et al., 2011; Taylor et al., 2013), and reflects the fact that as the relative surface area of the vine canopy increases (and the background signal decreases), there is less difference between sensors that only measure the canopy (proximal sensors) and those that capture the entirety of the vineyard (vines, weeds and soil) and have ‘mixed pixel’ effects. Although there were differences in the absolute recorded NDVI values between sensor systems, mainly due to whether the system sensed the canopy-only or the canopy and the background, a stability in the inter-annual NDVI response was observed across all sensors.

The strongest correlations were observed between NDVI maps derived from similar platforms (either proximal sensor pairs or remote sensing pairs), with correlation coefficients being the highest for the proximal sensor NDVI (full details provide in Table S1). The weaker correlation coefficients between the UAV and Sentinel layers, both assessed using a ‘mixed pixel’ approach, suggest lower reliability between these platforms. However, previous analysis (Kasimati et al., 2021b) resulted in some data from the proximally sensed surveys being removed before analysis due to artefacts in the survey data (one SpectroSense + GPS and one CropCircle dataset for 2019 and two datasets collected with CropCircle for 2020). Erroneous data were possibly related to either equipment malfunction resulting in poor sensing during the survey, illustrating a major disadvantage of proximal sensing because the terrestrial sensors needed to be reattached each time and were more likely to be affected by operator error. Proximal sensing, especially side-mounted (Fig. 2) sensor set-ups are also more likely to be adversely affected by any canopy interventions during the season, such as hedging or defoliation to manage the ripening of the grapes. Canopy interventions will alter the characteristics of the canopy and consequently the measurements taken.

The difference between the strong correlations of the Sentinel NDVI layers with the other sensors in 2019 and the weaker correlations of these satellite layers with all other sensors in 2020, was a concerning outcome from the analysis. It also indicated a divergence between the satellite platform and the high-proximity terrestrial and UAV observations, even when these higher- spatial resolution data were upscaled and the correlations were performed at a similar scale to the satellite imagery. The reason for the lower performance of the satellite imagery in 2020 is unknown and there was no clear indication of system failure. The temporal series of satellite imagery in 2020 were strongly correlated indicating that it was not a ‘one-off’ effect. Khaliq et al. (2019) have previously discussed that satellite imagery resolution cannot be directly used to reliably describe vineyard variability and several research papers (e.g., Darra et al., 2021) have pointed out that the use of medium spatial resolution satellites for vineyard feature assessment is less accurate due to large heterogeneity and narrow line spacing. The results here demonstrate that discrepancies may arise with satellite data acquisition and this might be the reason that few viticulture studies have been reported using such satellite systems (Anastasiou et al., 2019; Erena et al., 2016).

Excluding the 2020 satellite NDVI data, the exploratory correlation analysis for 2019 and 2020 between different NDVI layers and grape yield and TSS showed that the canopy reflectance data had absolute correlations of |r|> 0.60 with both attributes. Both attributes tended to have an increasing strength of correlation with NDVI (canopy vigour) as the season progressed, but a positive correlation with yield and a negative correlation with TSS. An inverse relationship between yield and TSS is not uncommon in viticulture, and reflects situations where crop load (the ratio of yield:leaf area) is unbalanced (Taylor et al., 2019) and the partitioning of photosynthate within the vines varies. Vines with higher yield are often ‘unbalanced’ and have less photosynthate available to convert into carbohydrates (sugars) resulting in a lower TSS.

In line with the weaker inter-season NDVI correlations, the Sentinel-2 NDVI showed much weaker correlations with both yield and TSS in 2020, and a divergence from the 2019 results, which were similar to the other NDVI sensors. The NDVI-yield correlations were strongest around veraison (early August data; Table 3), while the late season NDVI data had the best correlation with TSS (Table 4). Similar results that showed that NDVI at late developmental stages had good correlations with crop yield and TSS attributes were also found by other researchers in Greek vineyard systems (Anastasiou et al., 2018; Fountas et al., 2014; Tagarakis et al., 2013). Mid to late-season NDVI-yield relationships have also been found by Garcia-Estevez et al. (2017) in Spain (veraison NDVI) and Sun et al. (2017) (pre-harvest NDVI) in California.

Among all the sensors, the strongest correlation with yield was observed for UAV-derived NDVI data (2019 and 2020) and with TSS for NDVI data collected with Sentinel-2 (2019) and SpectroSense + GPS, (2020) in the middle/late season with full canopy growth, for both years. While TSS was correlated with the NDVI data (and yield), this was not the case for the other two important grape quality attributes, TA and pH, which had low correlation values (− 0.30 < r < 0.30) with the NDVI data at all crop stages. The development of acidity in these vines appears to be independent of the vigour:yield ratio in these systems.

The performance of each sensor varied and was influenced by the parameters of data collection, such as proximity to the vines and the specific technical characteristics of the equipment used. The UAV and SpectroSense + GPS appeared to perform better in the exploratory analyses and both systems were characterised by a similar mode of operation, i.e. the scan orientation was from above the vine canopy and at high resolution. Although the UAV system is classified as a remote sensor, it provides very high spatial resolution. The value of this over-head high resolution imaging is in agreement with previous studies in wine grape vineyards that have within-season canopy management (Dobrowski et al., 2008; Hall et al., 2011); however, it differs from Taylor et al. (2013), who observed saturation effects with nadir, above canopy proximal sensing systems in high-wire, low-intervention, sprawl juice grape vineyard systems in NY.

Although interesting and informative from a statistical and research perspective, correlations between NDVI and crop production attributes have limited value from an operational and vineyard management perspective. In contrast, deriving management units (MUs) is useful as it presents an alternative way to manage the vineyard. After converting the NDVI layers into MUs, i.e. expressed on an effective production scale, some clear patterns were observed. First, the PCA approach was very similar to the ‘All’ reference approach. The use of PCA, or some other data reduction method, could lend itself to automating this form of the process. The MUs created for low and high NDVI classes had greater consistency and showed very high similarity for the 2 years, with both the remote and proximal sensors showing very similar patterns of variation. Therefore, regardless of the sensing platform used, areas of high and low vigour could be identified for management. In contrast, the MUs of the intermediate NDVI classes were more variable between sensors. This raises an interesting question from a production perspective: Does the farmer want to better manage only the high and low vigour areas or all vigour areas? If the former, then the high-resolution maps in that vineyard can be considered similar, if not equivalent. If management of the entire vineyard is required, then they are less equal.

The DoA of the Sentinel-2 NDVI data with the reference and the other sensors is of concern, as it failed to detect the low NDVI area at the western edge in 2019 and produced a much more random pattern in the low and mid NDVI class in 2020 (Fig. 6). In contrast, it detected the reference high vigour area very well. These results were not due to poor quality imagery at any particular date, as all temporal Sentinel-2 NDVI values were highly correlated. Rather, it suggested that the satellite NDVI may not be ‘the same’ in the production context. Higher spatial resolution satellite imagery, such as imagery provided by PlanetScope (Planet Labs PBC, San Francisco, USA), could be an alternative.

The individual ANOVA results (Table 8) generally supported the notion that MUs derived from UAV data explained the greatest yield variation for both 2019 and 2020, followed by ‘All’ and PCA-derived MUs. This means that it made sense to use the UAV data in zoning for yield management. This was the case when considering both direct correlations and with zoning. For the variation in TSS, PCA, followed by SpectroSense + GPS and ‘All’, explained the largest part of variability for 2019 and UAV for 2020. Thus, in 2019, zoning using SpectroSense + GPS outperformed UAV, and in 2020, the inverse was observed. However, the direct Pearson’s correlation showed that SpectroSense + GPS outperformed UAV. Both sensors explained ~ 10% more variability in TSS than the other two systems (CropCircle and Sentinel-2), indicating that high-resolution NDVI data from above the canopy (nadir positioning) in these vineyards is preferable to side-curtain sensing. From the aforementioned results, it is unclear if zoning with Sentinel-2 data, particularly with the 2020 imagery, is sensible in these systems. Finally, as for the correlation analysis, the NDVI-derived MUs did not adequately explain the other grape quality attributes, TA and pH for any sensing system.

The results of the ANOVA indicated that it was not always appropriate to derive common MUs for both quality and quantity attributes. The amount of variability in yield and total soluble solids that was explained by the MUs created from different approaches did differ, especially between years. Yield for both years was better explained by UAV-generated MUs, explaining 76% for 2019 and 59% for 2020 (Table 8), but the best zoning was achieved from the SpectroSense + GPS data in 2019 (Table 7). The approach used here of assuming a reference layer from using all the data (‘All’) was affected in 2020 by the large divergence in the Sentinel-2 NDVI data from the other sensing platforms. Consequently, the fuzzy k-means classification generated MUs that tried to account for this divergence and, ultimately, generated MUs skewed towards explaining the different spatial patterns in NDVI in the satellite imagery. An unexpected result for this was that the DoA analysis in 2020 generated a stronger result for the Sentinel-2 system (Table 7). As noted earlier, the reason for the divergence in differing spatial patterns between the satellite platform and the other platforms in 2020 was not obvious. It was not seen in 2019. When the Sentinel-2 data was removed and the ‘All’ reference MUs regenerated with the UAV, CropCircle and SpectroSense + GPS data only, the DoA results did not improve (data not shown), although the SpectroSense + GPS was slightly better than the UAV (the inverse of the Table 7 2020 results). The variance of the TSS data explained by the MUs was much lower than the level of yield variance explained. Again, the 2020 NDVI-MUs explained less of the TSS variance than the 2019 NDVI-MUs (Table 8), with only the UAV data providing reasonable ANOVA fits in both years (adj R2 > 0.3), but the UAV was not the best fit in 2019. The seasonal climatic data (not shown) or experiences of the vineyard manager did not show any clear reasoning for the poorer fits (correlations and MU analyses) in 2020 compared to 2019. Further machine-learning approaches are targeted for these data which may help to understand the seasonal difference.

Conclusions

The overall aim of this research was to perform synchronous multi-platform multi-temporal NDVI analysis with vineyard production parameters (grape yield and quality) at the operational level. The high-resolution overhead systems (SpectroSense + GPS and UAV) appeared to provide the best results for managing yield and berry soluble total solids quality in the middle and end of the season when canopy growth is full. In particular, UAV-based data seems to be more relevant for this system, while at the same time this platform provides the benefit of fast data collection. The operational productive context was examined by deriving management units from the NDVI time series, with UAV-derived MUs on average explaining the largest variation in yield for both 2019 and 2020. Approaches that mixed multi-sensor information ('All' and PCA-derived MUs) also performed well. Satellite data, on the other hand, was found to lack consistency, with significant differences between the two years, suggesting that the satellite NDVI may not be 'the same' in the production context. This phenomenon of differences between satellite and UAV data is unclear and cannot be directly explained within the scope of this study. Further research needs to be conducted to repeat the analysis with satellite imagery for additional seasons. The various MU maps showed that NDVI data collected from different sensors did provide different information for decision making, which would lead to different spatial management in the vineyard. Growers in these types of systems should prioritise high-resolution overhead mid-late season NDVI mapping for differential vineyard management.