Vineyard Variability Analysis through UAV-Based Vigour Maps to Assess Climate Change Impacts

Abstract

Climate change is projected to be a key influence on crop yields across the globe. Regarding viticulture, primary climate vectors with a significant impact include temperature, moisture stress, and radiation. Within this context, it is of foremost importance to monitor soils’ moisture levels, as well as to detect pests, diseases, and possible problems with irrigation equipment. Regular monitoring activities will enable timely measures that may trigger field interventions that are used to preserve grapevines’ phytosanitary state, saving both time and money, while assuring a more sustainable activity. This study employs unmanned aerial vehicles (UAVs) to acquire aerial imagery, using RGB, multispectral and thermal infrared sensors in a vineyard located in the Portuguese Douro wine region. Data acquired enabled the multi-temporal characterization of the vineyard development throughout a season through the computation of the normalized difference vegetation index, crop surface models, and the crop water stress index. Moreover, vigour maps were computed in three classes (high, medium, and low) with different approaches: (1) considering the whole vineyard, including inter-row vegetation and bare soil; (2) considering only automatically detected grapevine vegetation; and (3) also considering grapevine vegetation by only applying a normalization process before creating the vigour maps. Results showed that vigour maps considering only grapevine vegetation provided an accurate representation of the vineyard variability. Furthermore, significant spatial associations can be gathered through (i) a multi-temporal analysis of vigour maps, and (ii) by comparing vigour maps with both height and water stress estimation. This type of analysis can assist, in a significant way, the decision-making processes in viticulture.

1. Introduction

About 70% of the available worldwide clean water is used in agriculture [1]. Moreover, by the year 2050, there will have to be an estimated 70% increase in food production [1] to sustain Earth’s population. Therefore, to attain a sustainable agriculture, it is essential to ensure proper water management. Global warming evolution throughout the years means this phenomena is one of the major threats to agricultural production, also with effects on society [2,3,4,5]. Less precipitation, associated with more frequent and longer drought periods [6], ultimately leads to an increase in the use of water in agricultural activity. To improve water usage efficiency, the United Nations (UN) set sustainable development goals with the aim to create an expected increase in efficiency in all sectors by the year 2030. This will ensure sustainable extractions and the implementation of integrated water resources management [7]. It is crucial that the agricultural sector contributes to this effort by developing and implementing controlled irrigation management systems [8,9]. As such, it is necessary to have an efficient analysis of crops’ water status.

The enduring search for resource use optimization, risks reduction, and minimizing environmental impacts led to the emergence of precision agriculture (PA) [10]. To understand both spatial and temporal variabilities of a production unit, PA’s tools and technologies enable the acquisition and processing of large data volumes (e.g., image processing techniques, geo-statistical methods) [10,11]. The precision viticulture (PV) concept derived from PA involves applying different technologies to vineyard management and grape production [12,13]. However, grapevine (Vitis vinifera L.) development is strongly related to spatial heterogeneity, which depends on several factors to determine both its production and quality [14]. Some of the more relevant factors are soil quality and type, vegetation management operations, irrigation systems, nutritional status, pest and disease control, air temperature, and precipitation levels [13,15]. Changes in one of these factors may result in the occurrence of biotic and abiotic problems. Depending on its severity, it may result in a significant decrease in production or quality, and therefore, considerable economic losses [16]. The Douro Demarcated Region (DDR, north-eastern Portugal) spatial variability is high due mainly to the terrain’s topographic profile, climatic variations, and soil characteristics, which causes vineyards to be unique throughout the DDR [12].

In the last few years, due to their flexibility and efficiency in diverse environments, the use of unmanned aerial vehicles (UAVs) emerged in agriculture applications [17]. UAVs can acquire georeferenced data with a high spatial resolution while using different types of sensors (RGB, near infrared, multi and hyper-spectral, thermal infrared (TIR) and LiDAR) [18], which allow for the output of several digital products, such as ortho-rectified mosaics, digital elevation models (DEMs), land surface temperature, and vegetation indices (VIs) [18]. Indeed, their ability to carry different types of sensors make UAVs a suitable solution for agricultural applications. While multispectral sensors acquire data from the electromagnetic spectrum in the near and visible infrared region (400 to 1000 nm), thermal sensors can acquire data in the far infrared zone (5000 to 18,000 nm), where the reflection value of each pixel can be transformed into a temperature value [18]. Among the different VIs, which can be considered as a set of arithmetic operations applied in different bands used to extract different vegetation characteristics [18], the normalized difference vegetation index (NDVI) [19] must be highlighted as it is frequently used in agricultural applications to estimate different crop-related parameters: biomass [20]; canopy structure, leaf area index (LAI), crop management [21]; and mapping vigour zones [22]. Moreover, it was found to correlate well with grape quality properties [23]. As for temperature-based indices, they constitute a quick and practical way to estimate crop water status, therefore indicating the plants’ water content. The crop water stress index (CWSI) [24] is widely used in remote sensing to monitor plants’ water status and consequent irrigation management [25]. TIR-based indices were employed to different crops, such as olives [26], grapevines [27], cotton [28], wheat [29], rice [30], sugar-beet [31], and maize [32]. Remote sensing platforms can also be a helpful tool for a better understanding of spatial variability, which has a significant meaning in vineyard management activities. Actually, UAVs have already been used to, e.g., estimate the leaf area index [33,34], irrigation management and water stress mapping [27,35,36], disease detection and mapping [37,38], and detection of nutritional deficiencies [39].

UAVs have already proved to be a cost-effective and flexible alternative for remote sensing, within a PA context. They present an improved decision-making process to the farmer and provide greater flexibility, when compared to other remote sensing platforms [13].

As for PV, vineyards have significant areas occupied by elements other than grapevines (e.g., inter-row vegetation, man-made structures, vegetation that usually surrounds the plot, and grapevines’ shadows) [13,40]. These elements can be automatically identified by means of digital image processing methods. Indeed, several methods have been proposed to deal with UAV-based aerial imagery or with the resulting digital products from the photogrammetric processing. For example, grapevine segmentation [41,42], supervised and unsupervised machine learning [43], point clouds derived from photogrammetric processing [44,45], and DEMs [16,33,40]. Regarding VIs, they are one of the most common segmentation techniques applied in a remote sensing [46], mainly to segment a given image into two classes: vegetation or non-vegetation [47]. However, when considering vineyard vegetation, VIs acknowledges all types of vegetation without distinguishing grapevines from non-grapevines (e.g., inter-row vegetation). By using the DEM—or more specifically, the canopy surface model (CSM), which can be obtained by subtracting the digital terrain model (DTM) from the digital surface model (DSM)—quantifying and removing non-grapevine vegetation in a vineyard’s segmentation process can be done as plant height is provided [48].

While different digital outputs can be generated from UAV-based imagery, the amount of data and its complexity can be overwhelming for the common farmer to interpret. Straightforward useful crop-related information is needed. Vigour maps are an example where by using the NDVI, vegetation is classified into different classes according to its characteristics. By applying it to PV, grapevines’ vigour can be defined as the measure of the growth rate during a given time period (e.g., the growing season). This not only enables the classification of vineyard homogeneity zones [49], which is a way to represent the impact of both environmental conditions and soil fertility [50]. There have been some related works done in this area. Khaliq et al. [51] compared satellite imagery with UAV-based multispectral data in four different epochs of the grapevines’ vegetative cycle. Different comparisons were made by considering: (i) the whole vineyard, (ii) only the grapevines’ vegetation, and (iii) only inter-row areas. The authors reported that satellite multispectral imagery presented limitations due to the ground sampling distance (GSD, 10 m) and to the influence of inter-row information. Primicerio et al. [22] evaluated vigour maps produced for the whole vineyard and only encompassing grapevines’ vegetation by applying an automatic segmentation method [41]. Campos et al. [52] used UAV-based vigour maps to create prescription maps for vineyard spraying operations.

Studies supported by imagery acquired in one flight mission alone mainly focused on assessing non-grapevine vegetation removal when considering the whole vineyard, and in creating task-oriented vigour maps [22,52,53,54]. With reference to multi-temporal studies, there are those whose aim is to compare different growing seasons by evaluating biophysical grapevines parameters [54,55,56]. Furthermore, studies utilizing intra-season multi-temporal data, considered the whole vineyard information [57] or vineyard changes were not the main focus [51]. As found in Primicerio et al. [22], vigour maps using only grapevines’ vegetation showed a better representation of the variability within the vineyard. The spatial variability in grapevines’ water status can be assessed thought both multispectral and TIR imagery, where TIR imagery serves as an immediate way to estimate crops’ water status, while multispectral data can show cumulative water deficits [35]. As such, the TIR data has the potential to help understand water stress for near-real-time decision-making support [58]. By integrating TIR and multispectral data, datasets to study grapevines’ response to climate change [59] can be created.

This study aimed to evaluate vineyard vigour maps (NDVI) created using UAV-based multispectral imagery within a multi-temporal context and in different grapevines’ phenological stages. The main goal was to study grapevines’ vegetation dynamics during the growing season up until harvesting. Two approaches were used: (i) considering the whole vineyard area, and (ii) considering only automatically detected grapevines’ vegetation. Spatial assessment between the generated vigour maps, and grapevines’ canopy temperature and height data—obtained from UAV-based TIR and RGB imagery, respectively—were conducted with the objective to correlate vigour maps with potential grapevines’ water stress and canopy height. This allowed for the assessment of non-grapevine features when analyzing vigour maps.

The next section presents the study area and the methods used both for data acquisition and processing. Results are presented in Section 3 and discussed in Section 4. Lastly, the most significant conclusions are shown in Section 5.

2. Materials and Methods

2.1. Study Area and Environmental Context

This study was conducted in a 0.30 ha vineyard located in the University of Trás-os-Montes e Alto Douro campus, Vila Real, Portugal (41°17′13.2′′ N 7°44′08.7′′ W WGS84, altitude: 462 m), in the DDR (Figure 1). The vineyard (cv. Malvasia Fina) is trained in a double Guyot system, where each row has grapevines 1.20 m apart and there is 1.80 m distance in between rows. There is a total of 22 rows with a NE–SW orientation. Furthermore, it is a rainfed vineyard, with fertilization applied using foliar spraying and with phytosanitary management operations taking place throughout the entire season. Inter-row areas are composed of spontaneous vegetation, which is managed using mechanical interventions at least twice per season.

During the studied period (May to September 2018), a total of 170 mm of precipitation was registered, along with 590 mm of potential evapotranspiration. Mean values for maximum, mean, and minimum air temperatures were 29 °C, 20 °C, and 13 °C, respectively. Monthly values are presented in Figure 2. Higher air temperature values were observed in July, August, and September, while May and June presented higher precipitation values. In contrast, there was almost no precipitation in August. This environmental data was acquired using a weather station located some 400 m away from the study area.

2.2. UAV-Based Data Acquisition

RGB, multispectral and TIR imagery were acquired using both a DJI Phantom 4 (DJI, Shenzhen, China) and a Sensefly eBee (senseFly SA, Lausanne, Switzerland). The former is a low-cost UAV equipped with an RGB sensor (12.4 MP resolution) attached to a three-axis electronic gimbal. For the purpose of this study, it was modified to support a multispectral sensor: the Parrot SEQUOIA (Parrot SA, Paris, France). This sensor consisted of a four-camera array, which was able to acquire data in the green (550 nm), red (660 nm), red-edge (735 nm), and near infrared (790 nm) parts of the electromagnetic spectrum, with a 1 MP resolution. Moreover, a Sunshine sensor (Parrot SA, Paris, France) was also added to the UAV’s top. It is responsible for acquiring the irradiance conditions during the flight mission in the same spectral bands as the multispectral sensor and to geolocate the acquired imagery.

As for the Sensefly eBee, it is a fixed-wing UAV used to acquire TIR imagery with the thermoMAP (senseFly SA, Lausanne, Switzerland) sensor (between 7500 nm to 13,500 nm, with 640 × 512 pixels and a temperature resolution of 0.1 °C), with automatic in-flight thermal image-based calibration. Ground control points (GCPs), used for aligning the acquired imagery during the photogrammetric processing, were measured using a Global Navigation Satellite System (GNSS) receiver in real-time kinematic (RTK) mode based on the TM06/ETRS89 coordinate system (GCP’s location in Figure 1). While the multi-rotor UAV was used mainly due to its capability to survey areas at lower flight heights, which provides higher spatial resolution [18], the fixed-wing UAV surveyed a larger area, which included the studied area. Furthermore, the TIR sensor only operated as a fixed-wing UAV.

Data acquisition was conducted in five flight campaigns, from 17 May 2018 to 21 September 2018. Each flight campaign corresponded to distinct grapevine phenological stages: flowering (May and June), fruit set (July), veraison (August), and harvest (September). Details are presented in Figure 3. All flight campaigns were conducted near solar noon to minimize sun and shadow influences. Flights for both the RGB and multispectral sensors were done at a 40 m height, with a forward overlap of 80% and 70% side overlap between images. The GSD was approximately 1.8 cm for the RGB and of 4.4 cm for the multispectral imagery. Regarding flights for TIR imagery acquisition, they were carried out at a 75 m flight height, with a 90% forward overlap and 75% side overlap between images, resulting in an approximate 17.5 cm GSD. All flight campaigns utilized RGB and multispectral imagery, while TIR imagery was only acquired from F3 onward (see Figure 3), due to both in-field observations and the environmental context, since rainfall can induce an error in the remotely sensed grapevine water status in the subsequent days [60]. Moreover, a radiometric calibration was performed prior to each flight for the multispectral imagery using a reflectance panel provided by the manufacturer, along with the irradiance data from the sunshine sensor. Irradiance and reflectance data enabled a reliable radiometric workflow for the collection of repeatable reflectance data over different flights, dates, and weather conditions.

2.3. Data Processing and Parameters Extraction

Imagery acquired in each flight campaign was processed using the Pix4Dmapper Pro (Pix4D SA, Lausanne, Switzerland). This software makes use of structure from motion (SfM) algorithms to identify common points in the images. It can create point clouds, and by interpolating them, generate different orthorectified outcomes depending on the sensor used. Imagery from each sensor was processed in different projects. The default processing options for each sensor were applied, but point clouds were generated with a high-point density. Point cloud interpolation was achieved using inverse distance weighting (IDW) and by applying noise filters. The generated digital outcomes were: (i) RGB–orthophoto mosaic, DSM, and DTM; (ii) multispectral–VIs; and (iii) TIR–land surface temperature. By subtracting the DTM from the DSM, the CSM was obtained. From the multispectral imagery, the NDVI [19] was obtained using a normalization between the near-infrared (NIR) and red bands, as given in Equation (1).

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

(1)

The land surface temperature was used to compute the CWSI through the empirical model presented in Equation (2). It was based in the usage of canopy temperature, T_c, and the lower and upper canopy temperature limits (T_dry and T_wet), corresponding, respectively, to well-watered and non-transpiring leaves. These values can be directly obtained in the field or by using UAV-based thermal infrared imagery [61]. CWSI values can vary between 0 (no stress signs) and 1 (high levels of stress). In this study, T_wet and T_dry values were obtained as described in the work of Matese and Di Gennaro [23]: T_wet was obtained by wetting some leaves and immediately measuring their temperature, while T_dry values were obtained by applying petroleum jelly in the leaves and registering their temperatures after some minutes had gone by. Temperature values were measured using a handheld infrared thermometer (Shenzhen Jumaoyuan Science and Technology Co., Ltd., Shenzhen, China), with a ±1.5 °C precision and operating between 8000 nm to 14,000 nm.

$$\mathrm{CWSI}=\frac{{\mathrm{T}}_{\mathrm{c}}-{\mathrm{T}}_{\mathrm{wet}}}{{\mathrm{T}}_{\mathrm{dry}}-{\mathrm{T}}_{\mathrm{wet}}}$$

(2)

To remove non-grapevine elements from the acquired imagery, segmentation was performed by using the method proposed in Pádua et al. [62]. Both the CSM and the G% index [63], computed from the orthophoto mosaic, were used as inputs, and through thresholding, considering both vegetation and height thresholds, it identified all vegetation within a given height range. While G% was automatically obtained using Otsu’s method for thresholding, CSM used a defined height range. An accurate grapevine segmentation was obtained, filtering out non-grapevine objects such as soil and inter-row vegetation.

This method has already been used in a multi-temporal analysis of grapevines’ vegetation evolution throughout a season in two vineyard plots in Pádua et al. [64]. Similarly, in this study, the method to segment grapevines’ vegetation [62] was applied to evaluate the multi-temporal vineyard evolution when regarding grapevine area and canopy volume, as well as the inter-row vegetation area. The grapevine canopy volume was computed according to Pádua et al. [64], using the mean height of each cluster of pixels obtained during the segmentation process multiplied by its area, where the sum of the volume of each cluster represents the total vineyard volume.

As such, the orthorectified outputs from each flight campaign were used for different purposes. The grapevines’ vegetation was detected and then CSM, NDVI, and CWSI values from the detected parts were considered, while non-grapevine pixels were discarded. Within the scope of this study, three different approaches were tested to create vigour maps. Figure 4 describes the main steps in each approach. Moreover, vigour classes were set to low, medium, and high. The workflow consisted in loading the orthorectified outcomes, followed by the vineyard segmentation method, depending on the used approach. Then, vigour maps were created by a applying a mean filter to the image, using a 2 × 2 m sliding window. Data could then be normalized according to Equation (3) before the vigour map was created. Again, this last step depended on the approach being used.

The first approach relied on the usage of data from the whole vineyard. The outcome was directly smoothed and divided into three classes, using terciles. As for the second approach, it was similar to the first, but it only considered the grapevines’ vegetation. Lastly, the third approach, similar to the second approach, considered only normalized grapevines’ vegetation. Normalization was done based on the mean value of the 10% higher and lower values of the smoothed grapevines’ vegetation values. Then, three vigour classes are created by dividing the values in the normalized raster according to fixed thresholds: (i) values lower or equal to 0.4 were considered low vigour; (ii) between 0.4 and 0.7 were considered medium vigour; (iii) and values above 0.7 were considered high vigour.

2.4. Vigour Maps versus Spatial Statistics

Vigour maps obtained from each flight campaign were compared with the CSM and the CWSI using statistical techniques that consider geospatial variability. This comparison was done by converting the three vigour classes maps to a 4 × 4 m grid. The grid size was selected by considering the studied vineyard’s characteristics: each grid square was confined to two vine rows. This pipeline was proposed by Matese et al. [56]. Regarding the methods used in this comparison process, they were the local bivariate Moran’s index (MI) and the bivariate local indicators of spatial association (LISA) [65]. Local MI (LMI) is based in the Moran’s index [66], which measures the global data correlation. While a positive correlation represents similar values in the area’s neighbourhood, a negative value represents the opposite, and zero represents a random spatial agreement. Regarding the LMI, a value is provided for each observation through permutation. The local bivariate MI was used in this study to assess the correlation between a defined variable and a different variable in the nearby areas. In turn, LISA measures the local spatial correlation, providing maps of local clusters with a similar behaviour, which is based on MI. This way, spatial clusters and its dispersion can be assessed. Bivariate LISA (BILISA) [65] was used as in Anselin [67] to examine the spatial relationship between the CSM and CWSI and the vigour maps. This comparison was made using GeoDa software [68]. Spatial weights were necessary to perform these analyses: a eight-connectivity approach (3 × 3 matrix) was used to create the weights map and BILISA was executed with 999 random permutations. The computed cluster maps and its significance were used. Cluster maps specify positive and negative spatial associations and are divided into four classes, based on the correlation of the value with its neighbourhood. The obtained associations are: (i) high–high (HH), where high values correlated with high values in the neighbourhood; (ii) low–low (LL), in which low values correlated with low values in the neighbourhood; (iii) high–low (HL); (iv) and low–high (LH). The three classes of vigour maps computed through the different approaches were compared with their correspondent vigour map in the following flight campaign, as well as with the CSM and CWSI three classes maps.

3. Results

This study yielded different digital products through the methods employed, from which it is important to highlight the vineyard status, vigour areas, potential water stress areas, and a multi-temporal vineyard characterization.

3.1. Multi-Temporal Vineyard Characterization

Figure 5 presents the orthorectified outcomes from the photogrammetric processing. There was a noticeable overall NDVI decline throughout the season (Figure 5a). However, grapevines’ canopy height (Figure 5b) presented a growth from the first to the third flight campaign, while remaining constant from then on. As for the temperature (Figure 5c), a high temporal variability was observed due to both the day temperature and the inter-row vegetation. For example, in the third flight campaign, temperature differences between areas with or without grapevines’ vegetation were smaller, about 1.0 °C, than in the other flight campaigns: approximately 2.2 °C for F4 and 1.4 °C for F5. Moreover, registered land surface temperatures presented the same behaviour as the maximum air temperature (Figure 2) registered in each month. Indeed, they were lower in July (followed by September), and higher in August.

Due to early vegetation development in grapevines by the time the first flight campaign took place, the minimum height to consider as grapevines’ vegetation was 0.2 m. As for the remainder of the flight campaigns, minimum and maximum heights were set to 0.5 and 1.9 m, respectively.

Table 1. Maximum, mean, and minimum values of the normalized difference vegetation index (NDVI), crop surface model (CSM), surface temperature, and crop water stress index (CWSI) when considering the whole vineyard plot and only grapevines’ vegetation in the five flight campaigns.

Type	Outcome	Parameter	F1	F2	F3	F4	F5
Whole area	NDVI	Max	0.88	0.91	0.89	0.78	0.78
		Mean	0.57	0.74	0.68	0.42	0.38
		Min	0.13	0.26	0.27	0.17	0.01
	CSM (m)	Max	1.17	1.48	1.59	1.51	1.53
		Mean	0.06	0.19	0.35	0.22	0.19
		Min	0.00	0.00	0.00	0.00	0.00
	Temp (°C)	Max	–	–	38.74	59.90	45.84
		Mean	–	–	29.89	44.35	37.20
		Min	–	–	27.12	37.26	32.49
	CWSI	Max	–	–	1.00	1.00	1.00
		Mean	–	–	0.60	0.83	0.78
		Min	–	–	0.04	0.23	0.07
Grapevines’ vegetation only	NDVI	Max	0.87	0.89	0.89	0.75	0.78
		Mean	0.70	0.82	0.80	0.62	0.59
		Min	0.41	0.59	0.64	0.37	0.25
	CSM (m)	Max	1.07	1.48	1.59	1.51	1.53
		Mean	0.40	0.89	1.16	1.01	0.99
		Min	0.20	0.47	0.52	0.27	0.20
	Temp (°C)	Max	–	–	31.20	47.81	39.36
		Mean	–	–	28.92	42.17	35.84
		Min	–	–	27.12	37.26	32.49
	CWSI	Max	–	–	0.82	1.00	0.91
		Mean	–	–	0.38	0.68	0.48
		Min	–	–	0.04	0.23	0.07

presents the differences in NDVI, CSM, land surface temperature, and CWSI values when considering the whole vineyard plot and when analyzing only detected grapevines’ vegetation. Generally, mean and minimum NDVI values were higher when considering only grapevines’ vegetation. As for maximum values, some high values were accounted for in areas other than with grapevines’ vegetation. The same tendency was verified in the mean and minimum height values, obtained through the CSM. However, maximum values were practically similar, except for the first flight campaign. An inverse tendency was verified when analyzing the land surface temperature and CWSI, i.e., higher values were found when analyzing the whole vineyard plot.

Extracted vineyard parameters allowed for a multi-temporal analysis of both grapevines’ vegetation area and volume, as well as for other vegetation present in the studied area. Figure 6 contains these results. The first flight campaign presented the lower values for the grapevines’ vegetation area: 82 m², representing 3% of the vineyard plot. The grapevines’ vegetation area increased until the third flight campaign, from which a significant decline was verified in the following flight campaigns. The grapevines’ canopy volume presented the same behaviour. As for inter-row vegetation, a growth happened between the first and the second flight campaigns, from 6% to 20% of the vineyard plot. After the fourth flight campaign, inter-row vegetation area decreased to 26 m² (1% of the vineyard plot), whilst a small increase was verified in the last flight campaign.

3.2. Generated Vigour Maps

Vigour maps were generated as described in Section 2.3 and assessment values are presented in this section. Each map was classified as one of three classes, namely as a low, medium, or high vigour area.

3.2.1. Visual Assessment

Figure 7 presents the vigour maps generated using three approaches. When encompassing the whole vineyard (i.e., considering bare soil and all existing vegetation), as presented in Figure 7a, a perspective of the plot’s homogeneity throughout the season was obtained. Approaches considering only detected vineyard vegetation presented a higher diversity, providing a deeper perspective on the grapevines’ vegetation spatial variability (Figure 7b,c). Still, a tendency for a lower vigour classification in the left part of the studied area was noticeable in all approaches. The same situation was verified in the southern central part of the vineyard plot. This assessment was more pronounced in the first approach but had more detail in both the second and third approaches.

Vineyard areas classified with high, medium, or low vigour were evaluated in all flight campaigns. Their percentages are presented in Figure 8a. As for the first approach, the vineyard plot showed a higher percentage of vegetation in the high vigour class (mean overall percentage of 48%). However, in the first flight campaign, there was a higher area classified in the low vigour class (mean overall percentage of 31%). The medium vigour class presented the lower mean overall percentage (21%). As for the second approach, the overall mean area percentage was similar: 43% in the high vigour class, followed by 33% in the low vigour class and 24% in the medium vigour class. Regarding the third approach, the medium vigour class presented the higher mean overall occupation area (42%), followed by the high vigour class (31%) and the low vigour class (27%).

The vineyard vigour area behaviour may not correspond to the grapevines’ vegetation. As such, Figure 8b shows the grapevines’ canopy volume present in each class throughout all the flight campaigns. This was achieved by intercepting vigour classes with the detected grapevines’ vegetation canopy volume. There were variations when comparing the applied approach and when analyzing the flight campaigns in the same approach: the overall value corresponded to the grapevines’ canopy volume presented in Figure 6. When considering the NDVI values for the whole vineyard to generate a canopy map, the grapevines’ canopy volume presented a higher predominance in the high vigour class. However, when comparing this with the approaches that consider only the grapevines’ vegetation, the grapevines’ canopy volume was significantly lower in the low vigour class for the latter approach. Regarding the approach where only grapevines’ vegetation was considered, a clear distinction among the grapevines’ vegetation volume was clear: the high vigour class had a greater grapevines’ canopy volume, followed by the medium and low vigour classes. As for the third approach (normalized grapevines’ vegetation), in the last two flight campaigns (F4 and F5), there was a higher volume in the medium vigour class, corresponding to the detected vineyard area (Figure 8a).

3.2.2. Spatial Correlations

To undergo a spatial assessment, the three approaches to generate vigour maps were applied to the CSM and CWSI outcomes of each flight campaign, when available. Ergo, maps with height values sorted in classes—low, medium and high height—could be obtained from the CSM. These results are presented in Figure 9. Height maps presented a high homogeneity among all approaches, especially from the third flight campaign onward.

From the CWSI, maps that could potentially point out grapevines’ water stress were obtained. They are presented in Figure 10. Again, three classes were considered to sort out each value on every map: low, medium, and high water stress. Results from considering all vegetation present in the vineyard (Figure 10a) showed a high homogeneity across the plot for all flight campaigns. However, when considering only grapevines’ vegetation (approaches two and three) the behaviour was different (Figure 10b,c).

Maps presented in Figure 9 and Figure 10 were compared with the vigour maps presented in Figure 7 in a 4 × 4 m grid using the LMI to measure their spatial correlation.

Table 2. Quantitative comparison using the local Moran’s index of the normalized difference vegetation index (NDVI) vigour classes in the three different approaches considered to the crop surface model (CSM) and crop water stress index (CWSI) classes with a p-value < 0.001, for each flight campaign (F#).

Vigour map	Approach 1		Approach 2		Approach 3
F#	CSM	CWSI	CSM	CWSI	CSM	CWSI
1	0.32	–	0.39	–	0.35	–
2	0.53	–	0.50	–	0.50	–
3	0.41	0.44	0.37	0.43	0.36	0.41
4	0.65	0.40	0.67	0.63	0.70	0.66
5	0.59	0.39	0.66	0.59	0.67	0.57

presents these results. Considering all the vineyards’ vegetation (first approach), stronger correlations were observed for the CSM. In turn, the other two approaches presented a more balanced trend for the CSM and CWSI. Stronger correlation values were found among vigour maps using data from the fourth flight campaign with the third approach (LMI = 0.70 for the CSM and LMI = 0.66 for the CWSI). Lower correlation values were observed in the height maps when considering all the vineyard’s vegetation with data from the first flight campaign. The same was verified in the fourth flight campaign for the CWSI.

The local spatial autocorrelation enabled the creation of clusters maps using BILISA to evaluate HH, LL, LH, and HL patterns between vigour maps of the different flight campaigns and between vigour maps and their correspondent height and water stress maps.

BILISA cluster map for the three evaluated vigour map approaches and its association with height maps is presented in Figure 11. As for the first approach (Figure 11a), there was a clear spatial correlation with a higher significance in the left and right sides of the vineyard plot, corresponding, respectively, to LL and HH associations. However, a smaller number of significant LH and HL clusters were detected. Regarding the other two approaches (Figure 11b,c) that considered only the grapevines’ vegetation, similar spatial patters were found for HH and LL. Furthermore, a significant HL cluster could be found in the southwestern part of the vineyard plot in the fourth and fifth flight campaigns. Significant LH clusters were found in the southeastern part of the vineyard in the second, third, and fourth flight campaigns for the second approach (Figure 11b).

Figure 12 presents the BILISA cluster maps generated from the spatial associations among vigour maps and water stress maps (Figure 10). Significant associations were found when using the first approach, with a representative HH cluster present in the northeastern region of the vineyard plot and a LL cluster in the vineyard’s left side. When considering only the grapevines’ vegetation, a similar behaviour was observed in the third flight campaign. In the remaining flight campaigns, a significant LL cluster existed in the left part of the vineyard, but a lower significance was found for HH in the northeastern part. A high significance among the values was detected in the southern region, which presented HH and LH associations.

Considering the BILISA clusters maps from the vigour maps for each evaluated approach when comparing consecutive flight campaigns (Figure 13), similar patterns were observed in all approaches and significant LH clusters were identified when comparing the first and second flight campaigns considering only the grapevines’ vegetation (Figure 13b,c).

4. Discussion

In this section the most meaningful results achieved in this study are discussed: (i) the multi-temporal analysis of the studied vineyard plot; (ii) the generated vigour maps; and (iii) spatial correlations between vigour maps, grapevines’ height, and potential water stress.

4.1. Multi-Temporal Analysis

The vineyard multi-temporal dynamics can be better understood using the orthorectified results obtained via photogrammetric processing of the UAV-based imagery (Figure 5) though their visual inspection in a geographical information system (GIS) [69].

Orthophoto mosaics can be used to detect missing grapevines and to manage vineyard in-field operations [64]. Vegetation indices (e.g., NDVI) can provide an overall assessment of vegetation vigour and potentially detect phytosanitary problems, such as flavescence dorée [37] and esca [70]. Leaf canopy temperature maps and CWSI can suppress the need to manually measure leaf water potential in the field [35]—a time-consuming approach, usually not performed in the whole vineyard—as well as be used for irrigation management [71].

In this study, an overall NDVI decline was noticeable from the third flight campaign onward (Figure 5a F4, F5). This was related to the grapevines’ vegetative cycle and to the decline of inter-row vegetation. Regarding height values obtained from each flight campaign’s CSM (Figure 5b), a clear distinction existed between grapevine and non-grapevine vegetation (e.g., soil and inter-row vegetation), except for in the first flight campaign (Figure 5b F1). Land surface temperature (Figure 5c) was clear-cut between flight campaigns. In fact, in the fourth and fifth flight campaigns (Figure 5c F4 and F5), there were some signs of the grapevines’ water stress.

Removing non-grapevine elements from the vineyard imagery provided a different perspective on the results, as confirmed in Table 1. Indeed, this enabled the production of estimate parameters such as the overall inter-row vegetation and the grapevines’ area and volume (Figure 6). The estimated grapevines’ vegetation area in the first flight campaign was 81 m² (3% of the vineyard plot) and in the second flight campaign, a 181 m² growth took place (262 m², 9% of the vineyard plot). As for the third flight campaign, there was a growth of 255 m² to 518 m² (18% of the vineyard plot). In the following flight campaigns, the grapevines’ vegetation area was reduced by 199 m² (−38%) to 319 m². Regarding the grapevines’ canopy volume, it was modified by +634%, +160%, −41%, and −12% in between each successive flight campaign, respectively. Concerning the inter-row vegetation area, it presented a behaviour consistent with the available precipitation data (Figure 2). Indeed, it had a 214% growth in between the first two flight campaigns (from 181 m² to 569 m²), representing 20% of the vineyard plot area. A steep decline was noticeable in the third and fourth flight campaigns (a decline of 95% to 26 m²), followed by a growth in the last flight campaign (88 m²). As such, the vineyard inter-row vegetation was a good indicator of soil water status. The same tendency had already been verified in Pádua et al. [64].

By comparing the mean, maximum, and minimum values observed in the different outcomes, either when considering the whole vineyard or only grapevines’ vegetation (Table 1), there was a clear difference among the flight campaigns. Mean NDVI values were superior in all flight campaigns when considering only the grapevines’ vegetation. The same tendency was verified in the CSM. This can be explained by the presence of a significant amount of lower values in the non-grapevine vegetation areas. However, the maximum NDVI values in the first, second, and fourth flight campaigns were registered in non-grapevine vegetation areas. Inter-row vegetation can account for this. Regarding maximum CSM values, they were similar in all flight campaigns, except for the first one, where the maximum height was detected in a non-grapevine area (probably a vineyard post). Minimum CSM and NDVI values were lower in non-grapevine areas. As for temperature-based outcomes (land surface temperature and CWSI), the opposite behaviour was found for the maximum values: they were located in non-grapevine vegetation areas. Mean temperature and CWSI values were lower in the grapevines’ vegetation areas, as it was expected due to the existence of bare soil areas in the vineyard. Minimum temperature and CWSI values were similar in both approaches since they were found in the grapevines’ vegetation areas. These results showed the importance of grapevines’ vegetation segmentation when analyzing a whole vineyard plot. The grapevines’ vegetation segmentation could improve the results in studies where this operation was not automatically performed, which is beneficial for removing non-grapevine elements from the analysis. Such an automatic procedure could help in the evaluation of vegetation indices [21], to detect flavescence dorée and grapevine trunk diseases [72], and to estimate grapevines’ biophysical and geometrical parameters [73].

4.2. Vigour Maps

Vigour maps generated when considering the whole vineyard provided an overall perspective (Figure 7a) about the studied area. Indeed, influences from bare soil and especially inter-row vegetation were clearly noticeable. Generally, the medium vigour class had the smaller area (Figure 8a) and the high vigour class encompassed the majority of the grapevines’ canopy volume (Figure 8b). The latter was, on average, 150% higher than the other vigour classes. A high homogeneity was verified for the last two flight campaigns. The same happened from the second to the last flight campaigns, when computing height maps from the CSM and all campaigns with CWSI. The whole vineyard was considered in both. Positive correlation values were found for the LMI (Table 2). Moreover, the verified homogeneity resulted in meaningful HH and LL areas when comparing vigour maps with CSM and CWSI in the same flight campaign.

Different results were obtained in the other two approaches, where only the grapevines’ vegetation was considered to create vigour maps. The higher incidence of missing grapevine plants in the left area of the vineyard remained almost the same throughout all flight campaigns. This was not noticeable in the first two flight campaigns’ vigour maps when considering the whole vineyard, probably due to an effect caused by inter-row vegetation. Other studies reported similar trends using vigour maps produced from the UAV-based NDVI [22,51] when excluding inter-row vegetation. Moreover, Vanegas et al. [74] found positive correlations when comparing vigour maps created from UAV-based data and a vineyard expert assessment.

As for vineyard area, when considering only grapevine vegetation, it presented a more balanced behaviour. The third approach, normalized grapevines’ vegetation, showed a considerable area of medium vigour class, particularly in the last three flight campaigns due to the fixed cut-off values to create vigour classes. Both approaches, grapevines’ vegetation and normalized grapevines’ vegetation, presented insignificant grapevines’ canopy volume values in the lower classes. Moreover, when considering normalized grapevines’ vegetation, canopy volume values were predominant in the medium vigour class, in agreement with the vineyard’s overall vegetative growth and decline (growth from first to the third flight campaigns and decline onward). Similar relations between vigour and the grapevines’ canopy volume were reported in other studies [73,75]. A higher heterogeneity was verified when observing both the CSM and CWSI maps generated with the approach that considered the normalized grapevines’ vegetation. In fact, when analyzing the CWSI maps from the last two flight campaigns (Figure 10), a period of the grapevines’ water stress was observed. However, this period was not clearly distinguishable in a visual map inspection based on data from the first approach (when the whole vineyard was considered). These correlations were observed in the BILISA cluster maps (Figure 12b,c), where areas with a high vigour showed a HL relationship with the CWSI maps, and significant agreements could be observed in the third flight campaign. A similar trend was reported in Matese and Di Gennaro [23]. Significant spatial associations were found in all approaches—whole vineyard, grapevines’ vegetation, and normalized grapevines’ vegetation—when analyzing the height class maps (Figure 11). Although lesser associations were found in the first flight campaign, this can be explained with the grapevines’ growth cycle. In this case, significant HL areas were found in the approaches considering only the grapevine vegetation. Similarly, Matese et al. [75] observed that some areas with a higher vigour were linked to areas with higher heights.

This study analyzed a vineyard’s behaviour throughout a season with a multi-temporal approach based on multispectral data acquired using a UAV. Furthermore, correlations between the different digital outcomes were found. This presents a potential tool for multi-temporal vineyard assessment and can serve as a base to provide prescription maps, similar to Campos et al. [52], since they can be correlated with agronomical variables (e.g., yield, berry weight, and total soluble solids), as shown in Matese et al. [56]. Indeed, patterns detected when comparing vigour maps from consecutive flight campaigns (Figure 13) highlighted differences in the multi-temporal data, which helps to understand local and spatial grapevines’ vegetative development dynamics throughout the season. However, filtered data considering only values representing grapevines’ vegetation, therefore representing the plants’ physiological status, was proven to be more reliable when comparing the evaluated approaches (Table 2); that is to say, it had a higher overall correlation. As such, it stands to be an excellent tool for decision support systems within vineyard management processes.

5. Conclusions

Climate change can heighten key environmental vectors that negatively impact vineyards. Grapevines can be weakened by both water stress and exposure to higher temperatures, which will increase their vulnerability to phytosanitary issues. UAVs equipped with different sensors can be used to regularly monitor grapevines, documenting changes in the vegetation or signs of diseases/infestation, as well as any stress caused by environmental constraints.

In this context, the need to evaluate current vineyard behaviour is crucial to proceed toward PV. Vigour maps can help to provide relevant insights, helping farmers and/or winemakers to understand their vineyards status and enabling timely actions to tackle problematic areas or observing response to treatments. Furthermore, the methods employed in this study to filter out non-grapevine vegetation presented a better vineyard representation, which can be used to assess a vineyard’s variability, but also to help in managing field-operations, such as those to inspect grapevines or to improve grapevines’ physiological status.

The use of methods to compare spatial correlations allowed us to obtain a spatial distribution of significant clusters among the different approaches evaluated for creating vigour maps. The importance of using different UAV-based outcomes to estimate biophysical and geometrical parameters shows the suitability of UAVs as a remote sensing platform for vineyard multi-temporal monitoring operations. This study allowed us to conclude that the need for UAV-based data can be tracked according to a vineyard’s phenology. Moreover, TIR data should be acquired in periods of higher temperatures to assess areas potentially affected by water stress. Nevertheless, the analysis presented in this study should be assessed in other vineyard types, such as those with irrigation systems, with a lower rate of missing grapevines, and in other wine producing regions with different grapevine training parameters.