Assessing the Feasibility of Using Sentinel-2 Imagery to Quantify the Impact of Heatwaves on Irrigated Vineyards

Heatwaves are common in many viticultural regions of Australia. We evaluated the potential of satellite-based remote sensing to detect the effects of high temperatures on grapevines in a South Australian vineyard over the 2016–2017 and 2017–2018 seasons. The study involved: (i) comparing the normalized difference vegetation index (NDVI) from mediumand high-resolution satellite images; (ii) determining correlations between environmental conditions and vegetation indices (Vis); and (iii) identifying VIs that best indicate heatwave effects. Pearson’s correlation and Bland–Altman testing showed a significant agreement between the NDVI of highand medium-resolution imagery (R = 0.74, estimated difference −0.093). The band and the VI most sensitive to changes in environmental conditions were 705 nm and enhanced vegetation index (EVI), both of which correlated with relative humidity (R = 0.65 and R = 0.62, respectively). Conversely, SWIR (short wave infrared, 1610 nm) exhibited a negative correlation with growing degree days (R = −0.64). The analysis of heat stress showed that green and red edge bands—the chlorophyll absorption ratio index (CARI) and transformed chlorophyll absorption ratio index (TCARI)—were negatively correlated with thermal environmental parameters such as air and soil temperature and growing degree days (GDDs). The red and red edge bands—the soil-adjusted vegetation index (SAVI) and CARI2—were correlated with relative humidity. To the best of our knowledge, this is the first study demonstrating the effectiveness of using medium-resolution imagery for the detection of heat stress on grapevines in irrigated vineyards.


Introduction
The increasing frequency of heatwave events represents a severe threat to viticulture. Prolonged periods of unusual heat might affect yields and the quality of grapevines. The vulnerability of grapevines to heatwaves varies with variety and phenological stage, with flowering and ripening being the most critical phases to heat stress [1,2]. Moreover, the occurrence of heatwaves during key phenological stages affects fruit-set and yield [3][4][5][6]. Reduced yields are probably related to the fact that heat stress affects photosynthetic rate, as demonstrated in Trincadeira in Portugal [7], Semillon in

Study Area
The study was conducted over the 2016-17 and 2017-18 growing seasons at Yalumba Oxford Landing Estate (OLE), a commercial vineyard in the Riverland, South Australia (34 • 06 06.29"S and 139 • 50 39.21"E). The trial area of 213 ha comprised several grapevine (Vitis vinifera L.) cultivars planted at a spacing of 1.8 m between vines and 3 m between rows (1852 vines ha −1 ), with an east-west row-orientation. The interrow was bare (no vegetation) during the growing season. The vines were trained using a quadrilateral cordon system and were mechanically box pruned with a sprawling canopy. Soil type was predominantly loamy sand (5-25 cm) with a sandy loam subsoil.
A preliminary analysis of NDVI values over the observation period revealed a non-productive area (12.07 ha, NDVI < 0.6 on 17 January 2017), a medium-vigor area (77.75 ha, NDVI from 0.6 to 0.8 on 17 January 2017) and a high-vigor area (111.90 ha, NDVI > 0.8 on 17 January 2017) ( Figure 1). By this date (17 January), the canopy reached its maximum size (shoot growth terminated). Therefore, it was considered an appropriate time to assess the variability of vigor within the sample area. The non-productive area, corresponding to headlands and nursery, was excluded from the analysis, which was conducted on the remaining 189.65 ha. The NDVI thresholds were established based on actual pixel values and through the examination of values found through reviewing references. Based on a study conducted in Italy, a minimum NDVI value of 0.6 appeared useful to identify a vegetation wall [52]. King et al. [53] used the same thresholds selected in the present analysis to identify the variability within the vineyard. The region is characterized by a dry, subtropical climate, and is mostly semi-arid. The temperature range can vary significantly both diurnally and seasonally. The region is prone to hot spells during summer, while winters are mild. The hottest month is January. According to the South Australian Regional Office of the Australian Bureau of Meteorology (BOM, [54]), the highest temperature (48.2 • C) was recorded in 2009. The average annual rainfall is 270 mm, mostly in winter.
The 2016-17 and 2017-18 growing seasons were similar in terms of heat accumulation, solar radiation and air temperature, while rainfall was more abundant during the first season (Table 1). According to the definition provided by BOM [54], a heatwave is the persistence of five or more days at or above 35 • C, or three or more days at or above 40 • C. This definition was however not specific to grapevines. Grapevine key physiological processes, specifically maximum photosynthetic rate, are compromised beyond 35 • C [55]. Thus, the temperature threshold for grapevines for the purposes of this study was revised to 35 • C, and heatwaves were redefined as three or more consecutive days at or above this temperature. Based on this temperature threshold, four heatwaves were observed in each season for 2016-2017 and 2017-2018 (Table 2).

Meteorological Data
The environmental parameters considered for the analysis of the correlations with spectral bands/VIs (Section 2.5) across the growing season are reported in Table 3. These data were obtained from the automatic weather station of the Natural Resources Management (NRM) weather station network located in Qualco, South Australia (https://www.awsnetwork.com.au/station/2770). Environmental data used for the analysis were based on the same day of the satellite image, the day before the satellite image, and the average of three days before the satellite image. These three observation times intended to consider both the immediate effect of heat stress and the time-lag effects of vegetation responses to weather conditions [56,57]. Differences between two successive time points (referred to henceforth as "variation" and denoted by "∆") were also calculated for all environmental parameters. Based on the availability of images, the second analysis (Section 2.5) instead used the average values from 12 days prior to the satellite image of the meteorological parameters, shown in Table 4.  Growing degree day GDD * RH = ambient relative humidity, VPD = vapor pressure deficit, T = ambient (air) temperature. Growing degree day last 12 days GDD -12 Cumulative sum (∆) of max T above 35 • C last 12 days ** ∆ 35 • C -12 * RH = relative humidity, VPD = vapor pressure deficit, ET 0 = evapotranspiration, T = temperature. ** Assigned the value zero to the air temperature of 35 • C (basal temperature), the ∆ 35 • C was calculated as the difference between the basal and daily maximum temperatures during the 12-day windows.

Data Acquisition and Processing
Thirty-one Sentinel-2 images were acquired from 9 October 2016 to 13 March 2018: 14 during the 2016-2017 growing season and 17 during the 2017-2018 growing season (see Supplementary Materials). The Sentinel-2 [58] multispectral sensor covers 13 spectral bands, ranging from 443 to 2190 nm, with spatial resolutions of 10 m (four visible and near-infrared bands), 20 m (six red-edge/shortwave-infrared bands) and 60 m (three atmospheric correction bands). To date, to date Sentinel-2 provides the best spatial resolution among open-source satellite images. The images were pre-processed using the SEN2COR tool [59], and available in the Sentinel-2 SNAP (Sentinel Application Platform) toolbox [60] to perform the atmospheric correction (level 2A images). This procedure ensured the conversion of top-of-atmosphere reflectances into top-of-canopy ones. After the correction and removal of border pixels, the NDVI values of each date were calculated using R statistical software (Version 3.5.2, RStudio Version 1.0.463) [61]. This procedure allowed the filtration of reflectance values of the digital numbers (DNs), excluding those pixels with an NDVI below a designated threshold (<0.6 on 17 January 2017), as they were considered to be relative to non-productive areas or headlands. These pixels were used to generate a mask to apply to the whole set of images, allowing the exclusion of non-productive areas from the analysis.
The spatial resolution of 10 m did not allow for the discrimination between rows and interrows. When using mixed pixels, the risk of major influence of the interrow in the reflectance measurement is high. To reduce this bias, one WorldView-2 high-resolution image (17 February 2017) was acquired from DigitalGlobe [62] to calculate the NDVI values of the experimental site after removing the interrows. The WorldView-2 image partially matched with the Sentinel-2 image from 16 February 2017 (105 ha).
DigitalGlobe WorldView-2 images have a spatial resolution of 0.46 m (panchromatic imagery) and 1.84 m (multispectral imagery). The WorldView-2 image from 17 February 2017 was compared with the Sentinel-2 image from 16 February 2017. The high-resolution of panchromatic imagery enabled removal of pixels corresponding to interrows to yield NDVI pixels corresponding to grapevine canopies. The pixels corresponding to rows and interrows were selected based on a preliminary visual assessment by establishing a threshold (pixels with DNs from 180 to 290 were considered vines). The threshold was used to create a mask layer using the "Raster Calculator" tool in QGIS Madeira version. To improve the precision of the process, manual correction was applied to the resulting raster layer. A preliminary downsampling procedure of the WorldView-2 image was performed applying the nearest-neighbor algorithm [63] to allow the comparison of the different resolution images. The NDVI values of Sentinel-2 and WorldView-2 were compared performing a Pearson's correlation and Bland-Altman concordance test [64]. Similar studies tested the consistency of datasets with different spatial-resolution by implementing a Pearson's correlation test. However, Pearson's correlation determines the strength of a relation between two variables but not their agreement. Two datasets can have a good correlation but strongly disagree. The choice to perform a Bland-Altman test seemed, therefore, more exhaustive.

Spectral Bands and Vegetation Indices
Spectral bands and VIs were computed for all the medium-spatial-resolution images using a custom script in R statistical software. The Sentinel-2 multispectral instrument covers 13 spectral bands of which eight were used in this analysis to compute VIs ( Table 5). The first analysis (Section 2.5), assessing the relation between environmental conditions and spectral bands, focused on the red edge, NIR (near-infrared) and SWIR bands. For the second analysis (Section 2.5) focusing on heatwaves, the ensemble of eight bands was tested.
A total of 13 VIs were calculated using the available spectral reflectance bands from Sentinel-2, chosen from the ones commonly used to assess water and heat stress. The VIs evaluated were both narrow-and broadband. Owing to the presence of bare soil in the interrow, the ensemble included four indices considering the effect of soil background.
The VIs assessed during the first analysis were NDVI, SAVI and EVI. The list of VIs used in the second analysis, along with the corresponding equations, is presented in Table 6. These VIs were chosen based on previously research showing their strong correlation with plant physiology under stress conditions.  Table 6. Vegetation indices used in the study calculated using a Sentinel-2 multispectral instrument set of bands.

Analysis of Remote Sensing Data
The study consisted of two different analyses. The first one aimed to assess the influence of environmental parameters on the spectral response of grapevines. The objective of the second analysis was to verify the effectiveness of medium-resolution remote sensing to identify and quantify the effects of heatwaves on grapevines.
The first analysis investigated Pearson's correlation coefficients between main environmental parameters (Section 2.2) and specific spectral bands/VIs (Section 2.4) across the growing season. This analysis aimed to provide baseline information on the relationship between key environmental parameters and the reflectance of the bands/VIs. The bands and VIs were chosen within the ones currently used in other crops to characterize heat stress [31,[73][74][75][76]. The examination was carried out on the entire trial area using the Sentinel-2 time series available over the two growing seasons, 2016-2017 and 2017-2018 (31 images). The correlations considered not only the environmental parameters but also their variation (∆) from the values recorded during the date of previous image (approximately 10 days).
The second analysis focused on heatwaves, which were defined as three or more consecutive days with daily maximum temperatures at or above 35 • C [54]. The evaluation was implemented by computing Pearson's correlation coefficients between the environmental parameters and spectral bands/VIs. Environmental parameters were assessed during two 12-day windows, the first window prior to the heatwave and the second following the heatwave and which included the heatwave period. The 12-day window was chosen based on the availability of the satellite images. For this analysis, only the images from pre-and post-heatwave were considered (11 images, Table 7). The correlation was calculated across the two different growing seasons and for the two seasons jointly. This second analysis was carried out at two levels: first on the whole trial area, and second on only the mediumand high-vigor areas (based on NDVI values) separately.

Statistical Analysis
Data analysis was carried out using R statistical software (Version 3.5.2, RStudio Version 1.0.463). Pearson's correlation was used to access the relationship between environmental parameters and spectral bands/VIs. The analysis of the correlations was carried out considering the high-vigor areas and the medium-vigor areas separately. The correlations were performed between aggregated values across the two individual growing seasons separately (five observations per each spectral feature and environmental parameter in 2016-2017 and six observations in 2017-2018), and for the two seasons jointly (11 observations per each spectral feature and environmental parameter). After investigating significant (R ≥ 0.6, P ≤ 0.05) correlations, the analysis also took into consideration significant and non-significant correlations with R ≥ ± 0.4 occurring in each growing season as well as combined for both seasons. Even if under significance level, recurring and consistent correlations might identify some trends worth exploring.
The Bland-Altman test was used to evaluate the agreement of WorldView-2 and Sentinel-2 images. Since the difference of the two measurements was not-normally distributed, the systematic bias and limits of agreements were estimated by computing the median and the first and third quartiles, respectively.

Comparison between Sentinel-2 and WorldView-2 Images
After the removal of the pixels relative to the interrows in the image from WorldView-2 Mission, the NDVI values of the medium-resolution image collected from Sentinel-2 Mission were compared to the ones of WorldView-2 Mission image oversampled at 10 m ( Figure 2). The aim of this analysis was to assess if the medium-resolution data could reproduce what can be observed using high-resolution images.
The comparison was implemented performing first a Pearson's correlation and then a Bland-Altman concordance test. The Pearson's correlation coefficient (R = 0.74, p < 2.22 × 10 −16 ) is shown in Figure 3. With regards to the Bland-Altman test, the estimated difference was −0.093 (95% confidence interval −0.127 --0.065, Figure 4). Differences between the two images tended to be larger for low NDVI values, but the concordance increased substantially above the threshold value of 0.5.   WorldView-2 and Sentinel-2 images. The plot describes the agreement between two quantitative measurements (NDVI). The systematic bias and limits of agreement were calculated by computingthe median and the first and third quartiles of the differences between the two measurements, respectively. The difference of the two paired measurements is plotted against the mean of the two measurements.

Analysis 1: Correlations between Environmental Conditions and Spectral Data
The results of the Pearson's correlation test highlighted significant correlations between the environmental conditions over the two growing seasons and some spectral bands and VIs. Specifically, the analysis highlighted a positive correlation between growing degree day (GDD) and NDVI (R = 0.62, p = 2.30833 × 10 − ), GDD and SAVI (R = 0.62, p = 2.3115 × 10 −4 ) and a negative correlation between GDD and the SWIR spectral band (R = -0.64, p = 1.17113 × 10 −4 ). Positive correlations were also found between (i) the variation (∆) of avg RH -1 and the red edge band 5 (R = 0.65, p = 1.29385 × 10 −4 ), and (ii) the ∆ of min RH 0 and EVI (R = 0.62, p = 3.40695 × 10 −4 ).
The NIR spectral band did not exhibit any significant correlation with the weather parameters.

Analysis 2: Spectral Features of the Heatwaves
After calculating the correlations across the two different growing seasons and for the two seasons jointly, the analysis focused on the recurring correlations (i.e., the correlations that showed up in both seasons). Both significant and trend results were considered. The correlations in the medium-vigor area, shown in Figure 5, tended to be higher, but the results were substantially homogenous across the three specific areas. Heat stress (avg T -12 , MAX T -12, GDD -12 , avg ST -12 , MAX ST -12, ∆ 35 • C -12 ) was better correlated with the spectral bands in the green and the red edge (bands 6 and 7) regions. The VIs that best described the spectral behaviors of the plants under heatwave conditions were TCARI and CARI. When considering avg ST -12 , EVI also showed a negative correlation trend, which was significant when considering the two growing seasons jointly (R = -0.73, p = 0.010).
The spectral bands more sensitive to a decrease in avg RH -12 during the heatwaves were the red edge 5 and the red. The correlation between red edge and avg RH -12 showed a trend in 2016-2017 and 2017-2018 and was significant when considering both seasons jointly (R = -0.63, p = 0.037). The same results were found considering the correlation between the red band and avg RH -12 (R = -0.63, p = 0.039) in both seasons considered jointly. The CARI2 VI exhibited a negative correlation trend with the relative humidity, which was significant in both seasons considered jointly (R = -0.73, p = 0.010), while the SAVI index positively correlated with avg RH -12 (R = 0.66, p = 0.028) in both seasons considered jointly.
The other VIs used in this study failed to produce useful results. Figure 5 shows the correlogram for the medium-vigor area over the two growing seasons considered separately and jointly. The weather parameters and the spectral bands/VIs not showing interesting results are not reported in the correlogram. Additional information is available in Supplementary Materials. The results relative to the high-vigor and the whole areas were consistent with those shown for the medium-vigor area. The correlograms for the high-vigor and the whole area are shown in the Supplementary Materials.

Discussion
This study analyzed eight heatwaves over the 2016-2017 and 2017-2018 growing seasons. The main challenge posed by the use of medium-resolution spectral information when analyzing row crops is the mixed nature of pixels due to the presence of interrows. The interrows likely affect the calculation of VIs and, therefore, crop status evaluation [77]. Borgogno-Mondino [78] utilized a Pearson's correlation test between high-resolution (0.5 m) UAV images and medium spatial resolution (15 m) satellite (Landsat-8) images on a vineyard after selecting pixels representing only canopy in the high-resolution image. The results showed a notable consistency between NDVI and vigor maps of the two datasets. Although the removal of interrow pixels is ideally the goal of imaging, this can become challenging in some contexts. The presence of a mid-row cover crop is not uncommon. Furthermore, interrows can have green covers. In this study, the validation of a Sentinel-2 dataset was implemented using two different tests-a Pearson's correlation and Bland-Altman concordance test-using a high-resolution image to compare NDVI values. Figure 3 shows that the correlation between the two images was statistically significant. The Bland-Altman test provided greater detail compared to simple correlation analysis, with the former showing that when observing low NDVI values the differences between the two images tended to be larger, but the concordance was satisfactory above the threshold of 0.5. (Figure 4). Considering that the Sentinel-2 image contained pixels with lower NDVI (as the interrows were not removed), these results indicate the feasibility of using spectral data from medium-resolution images to obtain useful information on vineyards.
The analysis of Pearson's correlation coefficients between environmental and spectral features (Analysis 1) provided a general framework to identify the spectral bands/VIs more sensitive to the changes in environmental conditions during the heatwaves. The positive correlations of GDD with NDVI and SAVI were not considered valuable information on the heatwaves, and it is likely that the higher values of the VIs at increasing GDDs was caused by a natural intensification of vigor over the growing season [20,[79][80][81].
The results suggested that monitoring RH during the heatwave periods was key to assessing in real time the effects of heatwaves on grapevines. Since RH indirectly influences plant water status [82,83], the positive correlations found for the red edge (band 5) and EVI indicate the importance of monitoring these spectral features during heatwaves. It has already been demonstrated that the red edge band has a good potential for estimating plant water status [84][85][86]. EVI was assessed in several studies to monitor drought stress [87][88][89]. This could be explained by changes in greenness-related features being less dynamic than modifications in vegetation water status, suggesting that indices in the SWIR spectrum region might be more efficient for drought stress detection [90,91]. The complexity of reflectance dynamics, especially for partial vegetation coverage areas, probably requires a combination of several VIs to monitor drought stress in crops [92][93][94]. The results of the current analysis suggest that a combination of EVI and SWIR might be a valuable indicator in identifying grapevine heat stress in bare interrow vineyards.
The specific analysis of the heatwaves (Analysis 2) focused on the correlations between spectral features and environmental parameters. The analysis aimed at monitoring the reflectance in different heat stress situations, considering images before and after the occurrence of individual extreme heat events. The Sentinel-2 mission reached the five-day revisit periodicity in March 2017. Before this time, the temporal resolution was 10 days. In some instances, cloud cover may have adversely affected the revisit time, leading to reduced availability of useful images. To overcome these limitations, the present study considered the average values of the previous 12 days (of the satellite flyover) for all meteorological parameters, therefore guaranteeing that the heatwave was always included in the analysis. The results of the experiment ( Figure 5) suggest that spectral bands and VIs can provide useful information about heat stress in vineyards. In particular, the green band was negatively related to the thermal parameters considered in this study (avg T -12 , MAX T -12, GDD -12 , avg ST -12 , MAX ST -12, ∆ 35 • C -12 ). The Pearson's correlation coefficients (R) were higher than -0.63 in all cases, achieving the highest correlation (R = -0.89, p = 0.044) with maximum soil temperature in the 2016-2017 growing season.
The reflectance of the green band is an indicator of healthy vegetation. Decreases in the green band, associated with chloroplast conformational changes, have been shown to be related to drought stress status [95]. As a result of the capacity of green reflectance to detect drought stress, several VIs using this band have been proposed. The Vis, including the green spectral band utilized for this study, were CARI, CARI2, MCARI, TCARI, NDVI and GNDVI. As shown in Figure 5, TCARI and CARI revealed good correlations with the thermal parameters, suggesting a feasibility of assessing grapevine heat stress using these VIs. Other spectral bands providing promising results were those using the red edge band, especially Band 6, thus confirming the results of Analysis 1.
The results of this study suggest that ∆ 35 • C -12 is a reliable indicator of heat stress in grapevines. Particularly, ∆ 35 • C -12 showed significant negative correlation (R = -0.91, p = 0.012 in 2017-2018; R = -0.81, p = 0.002 in both seasons considered jointly) with the NIR spectral band, which is typically used to evaluate plant water status of leaves and canopy [96][97][98]. With regards to humidity, the analysis highlighted a negative correlation trend between avg RH -12 and reflectance in the red spectral band. The increased reflectance in the red band is an early indicator of water stress [99]. The persistence of the negative correlation between the red band and the avg RH through the different growing seasons in all diverse vigor areas suggests that this spectral band could be used effectively for early detection of heat stress in grapevines.

Conclusions
This study investigated the feasibility of using free, medium-spatial resolution data from Sentinel-2 Mission to assess the effects of heat stress in grapevines.
The comparison between the medium-and high-resolution imagery showed a good agreement between the two datasets. The analysis of the spectral features of vineyards affected by subsequent heatwaves highlighted the spectrum regions and the VIs more suited to heat stress detection. The spectral bands and the VIs exhibiting high correlation with environmental conditions were consistent with the results of other studies on water stress performed using high-resolution imagery. These results suggest therefore that medium-resolution data can provide valuable information on the possible effects of heatwaves on grapevines. The possibility of using low-cost and open-access imagery could be used to support decision-making towards canopy management, irrigation scheduling, calamity alleviation, insurance services and recovery management. Further investigation is needed to confirm the results of this study, expanding the analysis to non-irrigated vineyards. Furthermore, it will be necessary to analyze the possible effects of heatwave timing, in relation to phenology, on the correlations between VIs/spectral bands and environmental parameters.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/11/23/2869/s1, Table S1: List of the satellite images used during the 2016-17 growing season; Table S2: List of the satellite images used during the 2017-18 growing season; Table S3: Pearson correlation coefficients (R and, in brackets, p) between spectral data and weather parameters in the medium-vigor area over 2016-2017 growing season. NR = not recurring correlation; red font = high significance (p < 0.05); blue fill = correlation (R > ±0.4) recurring in high-vigor, medium-vigor and whole areas; Table S4: Pearson correlation coefficients (R and, in brackets, p) between spectral data and weather parameters in the medium-vigor area over 2017-2018 growing season. NR = not recurring correlation; red font = high significance (p < 0.05); blue fill = correlation (R > ±0.4) recurring in high-vigor, medium-vigor and whole areas; Table S5. Pearson correlation coefficients (R and, in brackets, p) between spectral data and weather parameters in the medium-vigor area over the two growing seasons considered jointly. NR = not recurring correlation; red font = high significance (p < 0.05); blue fill = correlation (R > ±0.4) recurring in high-vigor, medium-vigor and whole areas; Figure S1: Correlogram of the input variables for high-vigor area in 2016/17 growing season (a), 2017/2018 growing season (b) and combined for both seasons (c). Only recurrent correlations between spectral features and weather conditions are shown. Positive correlations are displayed in blue and negative correlations in red colour. Colour intensity is proportional to R, while the magnitude of the circles is proportional to P; Figure S2