Using Satellite Thermal-Based Evapotranspiration Time Series for Defining Management Zones and Spatial Association to Local Attributes in a Vineyard

A well-planned irrigation management strategy is crucial for successful wine grape production and is highly dependent on accurate assessments of water stress. Precision irrigation practices may benefit from the quantification of within-field spatial variability and temporal patterns of evapotranspiration (ET). A spatiotemporal modeling framework is proposed to delineate the vineyard into homogeneous areas (i.e., management zones) according to their ET patterns. The dataset for this study relied on ET retrievals from multiple satellite platforms, generating estimates at high spatial (30 m) and temporal (daily) resolutions for a Vitis vinifera Pinot noir vineyard in the Central Valley of California during the growing seasons of 2015-2018. Time-series decomposition was used to deconstruct the time series of each pixel into three components: long-term trend, seasonality, and remainder, which indicates daily fluctuations. For each time-series component, a time-series clustering (TSC) algorithm was applied to partition the time series of all pixels into homogeneous groups and generate TSC maps. The TSC maps were compared for spatial similarities using the V-measure statistic. A random forest (RF) classification algorithm was used for each TSC map against six environmental variables (elevation, slope, northness, lithology, topographic wetness index, and soil type) to check for spatial association between ET-TSC maps and the local characteristics in the vineyard. Finally, the TSC maps were used as independent variables against yield (ton ha-1) using analysis of variance (ANOVA) to assess whether the TSC maps explained yield variability. The trend and seasonality TSC maps had a moderate spatial association (V = 0.49), while the remainder showed dissimilar spatial patterns to seasonality and trend. The RF model showed high error matrix-based prediction accuracy levels ranging between 86% and 90%. For the trend and seasonality models, the most important predictor was soil type, followed by elevation, while the remainder TSC was strongly linked with northness spatial variability. The yield levels corresponding to the two clusters in all TSC were significantly different. These findings enabled spatial quantification of ET time series at different temporal scales that may benefit within-season decision-making regarding the amounts, timing, intervals, and location of irrigation. The proposed framework may be applicable to other cases in both agricultural systems and environmental modeling.


Introduction
Precision irrigation for agricultural activities is a key factor for managing yield and quality, minimizing the use of water resources, and maximizing crop production. Irrigation practices generally relate to the timing, amounts, and spatial distribution of water input [1]. Therefore, both spatial and temporal information about various field characteristics is required to enable precision irrigation applications. Extensive research has been published regarding location and allocation of water input in the field, also involving the optimal timing of irrigation for specific types of crops [2][3][4]. Accurately quantifying the amount, timing, and location of irrigation depends both on the method of irrigation [5,6] and on accurate information regarding meteorology, terrain, plant, and soil conditions [7,8]. Site-specific management is an important factor in precision agriculture. The leading approach for differential management of the field is based on the delineation of the agricultural plot to management zones (MZs), which are homogenous sub-field areas based on within-field spatial variability [9]. There is an ongoing effort in the scientific community to develop and improve accurate and representative methods for defining MZs for specific purposes such as irrigation, fertilization, and pest control [10][11][12]. In vineyards, specifically for wine production, the amount of water provided to the grapevines is highly associated with vegetative growth, yield, and attributes related to quality traits, such as sugar content, acidity, and pigment content [1]. Water stress levels in vines are often highly variable across the vineyard [13,14]. Therefore, differential irrigation treatments are required, based on defined MZs, to comply with dissimilar water demands by the vines according to the target yield quality and quantity specified by the grower [15].
Evapotranspiration (ET) is defined as water loss to the atmosphere by evaporation from the soil surface and transpiration from the canopy cover. ET is a difficult water balance component to measure and estimate since it is affected by multiple atmospheric variables (e.g., wind speed, relative humidity, temperature and solar radiation) [16] as well as plant physiology [17,18]. The ability to quantify ET values in the field is highly valuable for determining the water stress of the crops and planning the amounts and timing of irrigation at different timesteps during the growing season [19,20]. Over the past two decades, many techniques for acquiring information regarding ET levels have been proposed and successfully applied [21][22][23]. However, these approaches lack the ability to provide routine spatial information about sub-field heterogeneity in plant conditions and water use/stress. In response, several satellite remote sensing-based algorithms have been developed and have demonstrated the utility of such an approach when attempting to quantify within-field heterogeneity in crop water use [24][25][26][27]. These algorithms include the Surface Energy Balance Algorithm for Land (SEBAL; [28]), Mapping Evapotranspiration with Internalized Calibration (METRIC; [26]), the Surface Energy Balance System (SEBS; [25]), and the Two-Source Energy Balance model (TSEB; [29]). The commonality between these approaches is that they predominantly rely on thermal infrared (TIR) imagery to approximate land surface temperature (LST) characteristics and then estimate ET as the residual component of the energy balance equation. These approaches have been used extensively and have been validated over many landcover types. A subset of ET mapping algorithms, including the Atmosphere-Land Exchange Inverse (ALEXI) model [27,30] and associated flux disaggregation technique (DisALEXI) [31,32] have been applied over vineyards, showing promising results [33][34][35][36][37].
The paired ALEXI/DisALEXI approach is predicated on the retrieval of LST from TIR imagery, which has been shown to be particularly useful in monitoring ET and plant stress due to its sensitivity to variations in soil moisture availability [38]. The approach calculates ET on both high-spatial/low-temporal (e.g., Landsat) and moderate-spatial/high-temporal (e.g., the Moderate Resolution Imaging Spectroradiometer; MODIS) resolution grids before being merged using the Spatial and Temporal Adaptive Reflectance Fusion Model (STARFM) [39]. Combining data from multiple sensors via data fusion provides daily ET estimates at field scale resolution, allowing for the Remote Sens. 2020, 12, 2436 3 of 23 quantification of both spatial and temporal dynamics in the field. Studies were conducted over two Pinot noir vineyards located near Lodi, CA, as a part of the Atmospheric Profile and Evapotranspiration eXperiment (GRAPEX; [35]). Since the inception of GRAPEX in 2013, the application of ALEXI/DisALEXI over vineyards has expanded to sites in Northern California (Sonoma County) and the southern Central Valley of California (Fresno County) [40].
Spatiotemporal dynamics in agricultural fields have been studied for decades [41][42][43]. The outcomes of these analyses were commonly a series of maps or descriptive statistics highlighting spatial variability in each time step [44]. In some cases, with a minimal number of time steps, this method is sufficient to emphasize the differences. However, in other studies where the timeframe includes many images/maps, an additional analysis was usually conducted to emphasize the temporal variation in space [45,46]. Quantification of ET over space is commonly represented on a daily scale, thus involving hundreds of images. Research in the past years has focused on mapping field characteristics, such as canopy cover and ET, at the field/landscape scale for each time step [37,47]. Alternatively, defining MZs based on within-field spatial variability provides information regarding the spatial heterogeneity in the field [48,49]. However, this approach does not provide any information regarding the temporal variation. Not many studies provided a framework to analyze the effect of spatially varying attributes on the temporal pattern [50]. For irrigation purposes, none of the mentioned approaches contributed an indication regarding the crop water status and the amount of input to be implemented in each MZ. Using daily ET assessments at the vineyard scale not only enables one to determine a spatial separation of the field into MZs based on ET temporal variations, but also to provide a complete status of the typical ET temporal dynamic of each MZ for irrigation purposes.
While time-series clustering (TSC) using remote-sensing data has been proposed in the past [51], very few studies have used this technique, mostly for land-cover monitoring [52]. To the best of our knowledge, no previous efforts have been made to use TSC for delineating MZs in agricultural fields, or to characterize the spatial attributes underlying the TSC spatial patterns.
The main objective of this study was to characterize the spatial and temporal patterns of ET within a vineyard. The specific objectives were: (1) to use time-series decomposition of ET into three components, including trend, seasonality, and remainders; (2) to cluster each of the three ET components into MZs using time-series clustering; (3) To quantify the impact of spatial field attributes on the spatial patterns of ET at different time scales (components); (4) to analyze the spatial effect of ET for each time component on grapevine yield.

Study Area
The vineyard under study is located in the Central Valley of California, near the city of Lodi (Figure 1), and is part of the USDA-ARS GRAPEX [35]. This area is characterized by a Mediterranean climate, with warm, dry summers and wet winters, with occasional spring rainfall events. The vineyard was planted in 2009 with Pinot noir vines and installed with drip irrigation. It has an area of 35 ha with an East-West row orientation and spacing of 1.5 m between vines and 3.3 m between rows [36]. Budbreak in the vines typically occurs at mid-March and harvest in late August/early September. The leaves then reach senescence and start to defoliate towards the dormancy period in the winter. Two vernal pools are located within the vineyard ( Figure 1) and were masked out since they do not represent the temporal structure of grapevine ET. In addition, the pixels at the borders of the vineyard were masked out as well, due to the potential impact on the border vines from external sources [53] and mixed-pixel effects. This was done by removing one column from the eastern and western parts of the vineyard, as well as the northern and southern rows. The vineyard is equipped with a meteorological/flux tower located approximately half-way north-south along the eastern edge of the field. Precipitation amounts were measured using a tipping bucket rain gauge (TE525, Texas Electronics, Dallas, Texas).

ET Images
This study applied the ALEXI surface energy balance model [27] and the associated disaggregation algorithm (DisALEXI) [31,32] in conjunction with STARFM [39] to derive daily images of ET over the vineyard of interest. ALEXI relates changes in LST, derived from geostationary satellites, to latent heat flux and surface moisture availability to estimate daily 4 km ET fluxes. These coarse resolution ET fluxes are spatially downscaled (DisALEXI) using higher resolution LST information from MODIS (1 km) and Landsat (30 m). They are then combined using STARFM to produce daily 30 m ET maps.
ALEXI is based on the TSEB model [29], which is unique in that it partitions radiometric temperature, TRAD, into soil and canopy temperatures (Ts and Tc, respectively), given the vegetation cover fraction (f) at radiometer view angle (θ): These component temperatures are then used to solve the energy balance equation for a mixed scene by combining soil and canopy contributions: where Rn is net radiation (W/m 2 ), H is sensible heat flux (W/m 2 ), G is soil heat flux (W/m 2 ), and λE is the latent heat flux (W/m 2 ). As in Equation (1), subscripts "s" and "c" refer to soil and canopy components, respectively. TSEB is applied twice during the morning, roughly 1 hour after sunrise and again 1 hour before local noon, using radiometric temperature acquired from the Geostationary Operational Environmental Satellites (GOES) platform. Energy closure during this time interval is provided by a simple slab model of atmospheric boundary layer development [54] which internally simulates land-atmosphere feedback on near-surface air temperature. As a result of using a timedifferential approach, the flux errors due to absolute sensor calibration, atmospheric corrections, and emissivity corrections are minimized [55].

ET Images
This study applied the ALEXI surface energy balance model [27] and the associated disaggregation algorithm (DisALEXI) [31,32] in conjunction with STARFM [39] to derive daily images of ET over the vineyard of interest. ALEXI relates changes in LST, derived from geostationary satellites, to latent heat flux and surface moisture availability to estimate daily 4 km ET fluxes. These coarse resolution ET fluxes are spatially downscaled (DisALEXI) using higher resolution LST information from MODIS (1 km) and Landsat (30 m). They are then combined using STARFM to produce daily 30 m ET maps.
ALEXI is based on the TSEB model [29], which is unique in that it partitions radiometric temperature, T RAD , into soil and canopy temperatures (Ts and T c , respectively), given the vegetation cover fraction (f ) at radiometer view angle (θ): These component temperatures are then used to solve the energy balance equation for a mixed scene by combining soil and canopy contributions: where Rn is net radiation (W/m 2 ), H is sensible heat flux (W/m 2 ), G is soil heat flux (W/m 2 ), and λE is the latent heat flux (W/m 2 ). As in Equation (1), subscripts "s" and "c" refer to soil and canopy components, respectively. TSEB is applied twice during the morning, roughly 1 h after sunrise and again 1 h before local noon, using radiometric temperature acquired from the Geostationary Operational Environmental Satellites (GOES) platform. Energy closure during this time interval is provided by a simple slab model of atmospheric boundary layer development [54] which internally simulates land-atmosphere feedback on near-surface air temperature. As a result of using a time-differential approach, the flux errors due to absolute sensor calibration, atmospheric corrections, and emissivity corrections are minimized [55]. ALEXI produces fluxes at 4 km spatial resolution, which is too coarse for field-scale applications. In order to provide field-scale flux estimates, the associated spatial disaggregation technique known as DisALEXI [32] was implemented. DisALEXI utilizes higher resolution TIR imagery (MODIS, Landsat, airborne sensors) to spatially disaggregate the coarser-resolution ALEXI flux estimates. This is done by running TSEB over each ALEXI pixel within a user-defined modeling domain, using TIR information at time of overpass from a satellite of interest (MODIS, Landsat, airborne sensor) as input. The air temperature boundary condition (set at a nominal blending height of 50m) for an ALEXI pixel is iteratively adjusted until the DisALEXI ET fluxes, spatially averaged over the ALEXI pixel, converge. This ensures consistency between the ALEXI and DisALEXI fields.
DisALEXI ET fluxes are upscaled from overpass time (instantaneous) to daily total estimates using the ratio of instantaneous to daily insolation: where λE t2 is the latent heat at t2, Rs t2 is the insolation rate at t2, R s24 is the time-integrated daily insolation rate, λ is the latent heat of vaporization, and ρw is the density of water. The reader is referred to [39,56,57] for detailed descriptions of the ET modeling process. Applications over the current vineyard (Lodi) are described by [33,36,40]. Previous studies report good agreement between ALEXI/DisALEXI-derived daily ET and eddy covariance flux tower measurements, regardless of vineyard location. Specifically, Semmens et al. [33] reported root mean square error (RMSE) of under 1 mm day −1 at the Lodi site for the 2013 growing season. Knipper et al. [36] reported slightly lower RMSE values (0.8 mm day −1 ) at the same site when extending the analysis through 2016. When compared over multiple vineyards characteristic of having different grape varieties, trellis designs, canopy structure, irrigation management and climate conditions, the ALEXI/DisALEXI model produced daily ET estimates with good fidelity (RMSE values of 0.88 mm day −1 on average between sites) [40].
ET images from four consecutive growing seasons, between 2015 and 2018, were used in this study. The growing season of grapevines is typically from early spring until mid-fall, while the vines are covered with transpiring leaves. Therefore, a subset was taken from the four-year dataset, from March 1 to October 31 (days of year (DOY) 60-304) of every year ( Figure 2). The selected images included 71 Landsat scenes and 980 MODIS scenes. Following the STARFM fusion process, there were 980 high resolution (30 m) ET images, with 245 images per growing season. Few images had missing ET values due to weather conditions that did not allow for clear MODIS images, and the missing data for each pixel was interpolated using a temporal moving average with a window of 4 days. This procedure was performed using the "imputeTS" package in R [58].
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 24 ALEXI produces fluxes at 4 km spatial resolution, which is too coarse for field-scale applications. In order to provide field-scale flux estimates, the associated spatial disaggregation technique known as DisALEXI [32] was implemented. DisALEXI utilizes higher resolution TIR imagery (MODIS, Landsat, airborne sensors) to spatially disaggregate the coarser-resolution ALEXI flux estimates. This is done by running TSEB over each ALEXI pixel within a user-defined modeling domain, using TIR information at time of overpass from a satellite of interest (MODIS, Landsat, airborne sensor) as input. The air temperature boundary condition (set at a nominal blending height of 50m) for an ALEXI pixel is iteratively adjusted until the DisALEXI ET fluxes, spatially averaged over the ALEXI pixel, converge. This ensures consistency between the ALEXI and DisALEXI fields.
DisALEXI ET fluxes are upscaled from overpass time (instantaneous) to daily total estimates using the ratio of instantaneous to daily insolation: where λEt2 is the latent heat at t2, Rst2 is the insolation rate at t2, Rs24 is the time-integrated daily insolation rate, λ is the latent heat of vaporization, and is the density of water. The reader is referred to [39,56,57] for detailed descriptions of the ET modeling process. Applications over the current vineyard (Lodi) are described by [33,36,40].
Previous studies report good agreement between ALEXI/DisALEXI-derived daily ET and eddy covariance flux tower measurements, regardless of vineyard location. Specifically, Semmens et al. [33] reported root mean square error (RMSE) of under 1 mm day −1 at the Lodi site for the 2013 growing season. Knipper et al. [36] reported slightly lower RMSE values (0.8 mm day −1 ) at the same site when extending the analysis through 2016. When compared over multiple vineyards characteristic of having different grape varieties, trellis designs, canopy structure, irrigation management and climate conditions, the ALEXI/DisALEXI model produced daily ET estimates with good fidelity (RMSE values of 0.88 mm day −1 on average between sites) [40].
ET images from four consecutive growing seasons, between 2015 and 2018, were used in this study. The growing season of grapevines is typically from early spring until mid-fall, while the vines are covered with transpiring leaves. Therefore, a subset was taken from the four-year dataset, from March 1 to October 31 (days of year (DOY) 60-304) of every year ( Figure 2). The selected images included 71 Landsat scenes and 980 MODIS scenes. Following the STARFM fusion process, there were 980 high resolution (30m) ET images, with 245 images per growing season. Few images had missing ET values due to weather conditions that did not allow for clear MODIS images, and the missing data for each pixel was interpolated using a temporal moving average with a window of 4 days. This procedure was performed using the "imputeTS" package in R [58].

Environmental Attributes
Six different spatial environmental attributes were acquired to check for their influence on the spatial pattern of the ET time series in the vineyard. These covariates included the following ( Figure 3):

•
Digital elevation model (DEM): Represents shifts in elevation (m) at a 30 m spatial resolution. The variability in elevation is known to affect plant conditions [59] and ET across space [60].
Within the study area, the elevation ranges between 37 and 47 m above sea level, with a median of 40.5 m above sea level. • Slope: Derived from the DEM, with identical spatial resolution and coordinate system. It was computed as the highest difference between a certain pixel and its 8 neighboring cells [61] and was applied using the "raster" package in R [62]. The slope is known to have an impact on soil water content [63], thereby affecting the water availability to the plant [64] and its water consumption levels. The slope in the vineyard is considered moderate, ranging between 0 and 3.16 • , with a median of 0.75 • . • Northness: Derived from the DEM with an inherent 30 m spatial resolution. Aspect was computed as the gradient vector at each pixel location according to the direction between neighboring pixels, such that the North is defined as 0 • and south as 180 • . The aspect was derived using the R package "raster" [62]. To avoid using a circular variable in a regression model, aspect was transformed to northness, which is the cosine of aspect [65,66]. Aspect is known to affect the levels of exposure to solar radiation by the plants, thereby influencing plant water conditions, with significant effects between North-and South-facing slopes [67]. The west side of the vineyard is mostly facing West, North, and North-West, while the East side of the vineyard is mostly facing East and South-East.

•
Topographic wetness index (TWI): Represents soil moisture conditions and quantifies the distribution of accumulated area across downslope pixels. It was calculated by Equation (4) using the "dynatopmodel" package in R [68]: where A stands for the upstream contributing area (m2) for each pixel and β is the local slope in the steepest downhill direction (degrees). TWI is used to assess the soil moisture and its effects on plant conditions [69]. TWI levels ranged between 6 and 17 with a median of 9. • Soil types: The soil type map represents different soil groups within the vineyard, represented as polygons, as classified in the Gridded Soil Survey Geographic (gSSURGO) Database. The soil types layer was spatially overlaid over the raster dataset to retrieve the indexes of the soil classes for each pixel, using R package "sp" [70]. Soil texture and structure are among the main attributes affecting soil moisture levels and soil water extraction patterns [71] due to variations in soil depth, permeability levels, micronutrients, and more [72,73]. In turn, these soil properties affect water availability to the plant and ET levels [74]. There are four soil type classes within the vineyard domain (SSHCP, 2018), including the following: Redding gravelly loam: Associated with the Laguna Formation, which is associated with high terrace landform, and composed of interbedded alluvial gravel, sand, and silt. This soil series developed on older, higher alluvial terrace landforms, is a well-drained, gravelly loam mineral soil that is moderately deep and formed in gravelly valley fill [75]. The water holding capacity in low and permeability is very slow, reaching 50-100 cm [72].
Hadselvill-Pentz complex: Consists of very shallow loamy, moderately well-drained intermound soils. This soil series is associated with the Mehrten Formation, which consists of eroded, volcanic mudflow fans. It has organic matter content of 1%-2% and is very shallow [75], with moderately rapid permeability, reaching 10-50 cm. Water holding capacity is very low [72].
Corning complex: Formed of gravelly loam from the Laguna Formation and gravels on high terraces. They have microrelief composed of a mound-intermound structure on smooth terrace side slopes. These are very deep soils, well-drained with a gravelly fine sandy loam surface Remote Sens. 2020, 12, 2436 7 of 23 layer [75]. Permeability in this soil series is very low. After rainfall events/irrigation, the soil is saturated for a short time [72]. Kimball silt loam: The Kimball series consists of well-drained loams. They were formed on older terraces in alluvium derived from sedimentary rocks [76]. Its subsoil permeability reaches 150+ cm, with a moderate water holding capacity and slow runoff [72]. • Lithology: The lithology map represents the spatial distribution of rock formation within the vineyard, represented as polygons, and the same technique was used to extract the classes representing each pixel as for the soil types layer. The lithology may have an effect on water availability to the plant [77] and its ability to consume water. This covariate included two classes -Sandstone on the northwest part of the plot and alluvium in the southeast part of the vineyard.
Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 24 [75]. Permeability in this soil series is very low. After rainfall events/irrigation, the soil is saturated for a short time [72]. Kimball silt loam: The Kimball series consists of well-drained loams. They were formed on older terraces in alluvium derived from sedimentary rocks [76]. Its subsoil permeability reaches 150+ cm, with a moderate water holding capacity and slow runoff [72]. • Lithology: The lithology map represents the spatial distribution of rock formation within the vineyard, represented as polygons, and the same technique was used to extract the classes representing each pixel as for the soil types layer. The lithology may have an effect on water availability to the plant [77] and its ability to consume water. This covariate included two classes -Sandstone on the northwest part of the plot and alluvium in the southeast part of the vineyard.

Yield
A yield map was also used for analyzing the association between ET spatial pattern and yield levels [21]. The yield map (Figure 4) was provided by E&J Gallo Winery using Advanced Technology Viticulture (ATV, Joslin, Australia) yield monitoring systems that were applied in the vineyard during the harvest in August 22, 2016. Post-harvest analysis converts mass flow units into tons per hectare, eliminates outliers (yield data greater than three standard deviations from the mean), and normalizes across harvesters. The yield map was interpolated to a 3-meter grid using Vesper Spatial Prediction Software for Precision Agriculture [81]. The Northern-most grapevine rows are masked out due to spurious detected artifacts. The DEM [78], soil types [79], and lithology [80] covariates were acquired from the Geospatial Data Gateway site (https://gdg.sc.egov.usda.gov/) of the United States Department of Agriculture (USDA).

Yield
A yield map was also used for analyzing the association between ET spatial pattern and yield levels [21]. The yield map (Figure 4) was provided by E&J Gallo Winery using Advanced Technology Viticulture (ATV, Joslin, Australia) yield monitoring systems that were applied in the vineyard during the harvest in August 22, 2016. Post-harvest analysis converts mass flow units into tons per hectare, eliminates outliers (yield data greater than three standard deviations from the mean), and normalizes across harvesters. The yield map was interpolated to a 3-m grid using Vesper Spatial Prediction Software for Precision Agriculture [81]. The Northern-most grapevine rows are masked out due to spurious detected artifacts.

Time-Series Decomposition and Clustering of Each Component
The ET time series of each pixel was decomposed into three different time components. In order to comply with specific objective 2, a TSC procedure was applied to each of the components. In general, TSC is a process during which time series data is partitioned into homogeneous groups based on the similarities of the time series values and temporally varying patterns [82]. The data composing a time series is dynamic, its values change as a function of time and clustering of time series may reveal patterns and changes [82]. TSC relied on four elements: Time series representations: A time series is a set of observations, each recorded at a specific time [83]. A continuous time series consists of observations recorded continuously over a time interval. A time series may be displayed as a linear combination of components. These components provide information about patterns of change over time of latent processes composing the data [84] and include the trend, seasonality, and remainder. The trend represents low-frequency variations of the time series, including non-stationary, long-term inter-annual changes. The seasonal component accounts for the seasonal frequency in terms of cycles, in our case, cycles per growing season. The remainder component is the remaining variation (anomaly) in the data beyond the variation accounted by the seasonality and trend components [85]. It may be regarded as the local processes which cause variability among cycles, or shortest-term events [86], while in this present study, the shortest time unit is the daily scale. Decomposition of the time series allows for analyzing different underlying patterns in the data and has been previously used for quantifying multi-seasonal ET patterns [87]. The three components are affected by the inherent characteristics of the time series, such as length and frequency. However, as the suggested framework deals with comparative analysis among the time series of the pixels, it may work on any time period, as long as it consists of two cycles (growing seasons) or more. In other words, the meteorological conditions and effects will be similar for all pixels, regardless of the timeframe that is being used, and what matters is how the plants in each pixel respond to these conditions compared to other pixels. In this study, time-series decomposition was applied in order to analyze each timescale separately; namely, the multi-seasonal timescale (trend), the seasonal timescale (seasonality), and the daily timescale (remainder). The time series in the present project included 980 images. Therefore, each pixel had 980 daily values, and the decomposition was performed at the pixel level. The decomposition method that was used is the

Time-Series Decomposition and Clustering of Each Component
The ET time series of each pixel was decomposed into three different time components. In order to comply with specific objective 2, a TSC procedure was applied to each of the components. In general, TSC is a process during which time series data is partitioned into homogeneous groups based on the similarities of the time series values and temporally varying patterns [82]. The data composing a time series is dynamic, its values change as a function of time and clustering of time series may reveal patterns and changes [82]. TSC relied on four elements: Time series representations: A time series is a set of observations, each recorded at a specific time [83]. A continuous time series consists of observations recorded continuously over a time interval. A time series may be displayed as a linear combination of components. These components provide information about patterns of change over time of latent processes composing the data [84] and include the trend, seasonality, and remainder. The trend represents low-frequency variations of the time series, including non-stationary, long-term inter-annual changes. The seasonal component accounts for the seasonal frequency in terms of cycles, in our case, cycles per growing season. The remainder component is the remaining variation (anomaly) in the data beyond the variation accounted by the seasonality and trend components [85]. It may be regarded as the local processes which cause variability among cycles, or shortest-term events [86], while in this present study, the shortest time unit is the daily scale. Decomposition of the time series allows for analyzing different underlying patterns in the data and has been previously used for quantifying multi-seasonal ET patterns [87]. The three components are affected by the inherent characteristics of the time series, such as length and frequency. However, as the suggested framework deals with comparative analysis among the time series of the pixels, it may work on any time period, as long as it consists of two cycles (growing seasons) or more. In other words, the meteorological conditions and effects will be similar for all pixels, regardless of the timeframe that is being used, and what matters is how the plants in each pixel respond to these conditions compared to other pixels. In this study, time-series decomposition was applied in order to analyze each timescale separately; namely, the multi-seasonal timescale (trend), the seasonal timescale (seasonality), and the daily timescale (remainder). The time series in the present project included 980 images. Therefore, each pixel had 980 daily values, and the decomposition was performed at the Remote Sens. 2020, 12, 2436 9 of 23 pixel level. The decomposition method that was used is the Seasonal and Trend decomposition using Loess (STL) [85] with seasonal window (the smoothing parameter for the seasonal filter) set to periodic (mean seasonal value), so that it was identical across all years. The trend window (the smoothing parameter of the trend) was set to the default 1.5 x period/1-(1.5/seasonal window). Figure A1, in the Appendix A, illustrates a time-series decomposition applied to a time series of one pixel. This process was conducted for all 334 pixels across the vineyard using the R package "stats" [88]. For the trend and remainder, the entire timeframe was used for the TSC procedure, while the seasonality component was reduced to one cycle using the "TSrepr" package in R [89].
Distance matrix: A scalar measure of dissimilarity between all pairs of observations (each observation is a time series of a specific pixel), or similarity matching, in this case each observation consists of daily values during four growing seasons. This whole-sequence matching considers the entire lengths of the time series during the distance calculation [82]. The dynamic time warping (DTW) algorithm is a technique applied to generate the distance matrix and was used to compute the similarity among time series (i.e., pixels). This dynamic programming technique is based on a non-linear warping function that eliminates timing differences between two temporal patterns by warping the time axis of one series so that the similar maximum values would coincide with each other. The time-normalized distance is then calculated as the minimized residual distance between the two series [90].
Clustering algorithm: Uses the distance matrix to group the data based on the similarities between each pair of pixels' time series. In this study, we used the unsupervised fuzzy c-means (FCM) algorithm as the method of clustering. A fuzzy set has fuzzy boundaries where each element receives a degree of membership to each set. The data were partitioned into c fuzzy groups, and the FCM finds a cluster center (centroid) in each group such that the cost function of dissimilarity is minimized [91]. The fuzzification parameter, m, determines the degree of fuzziness in the clusters and, in this case, was set to 2, as conventionally suggested by most studies [92]. The FCM algorithm was applied to the data with 2 ≤ c ≤ 6, and the final number of groups was assigned based on the clustering result with the highest separation between within-cluster similarities and between-clusters dissimilarities. c = 6 was defined as the maximum number of clusters since MZ > 6 would not be practical for irrigation purposes. This partitioning evaluation was performed using the silhouette index, which measures the degree of similarity of each observation to its own cluster compared to other clusters [93]. It ranges from −1 and 1, while higher values indicate a stronger match to the assigned cluster. The average silhouette width (ASW) for the observations was computed using the "cluster" R package [94].
Time-series clustering maps: The outcomes, or output, of the TSC were maps composed of different clusters, while each pixel was assigned to a specific cluster. For each of the three components, the TSC map represented the spatial pattern of temporal similarities among different pixels. The trend TSC map illustrated the long-term, inter-seasonal pattern of change for each group of pixels. The seasonality TSC map represented the seasonal similarity among each group of pixels, and the remainder TSC map displayed the daily within-field spatial pattern according to the dissimilarities between the different clusters. The TSC based on DTW was applied using "dtwclust" package in R [95]. Corresponding plots of the time series for each component were also provided to show the dynamic patterns in time for each cluster. An additional difference plot among the values of the different TSC was presented to show the magnitude of dissimilarity among the clusters using the R packages "reshape2" [96], "data.table" [97], "RStoolbox [98], and "ggplot2" [99].

Precipitation and Irrigation Data Analysis
The irrigation regime and precipitation amount (mm) were examined in an attempt to explain the ET trend pattern. The spring precipitation values from February until June were therefore extracted from the meteorological station located in the vineyard. The accumulated precipitation was calculated, as well as the number of days with rainfall events. In addition, to complement the total water input to the vineyard, the amounts of irrigation per growing season were calculated. Finally, a sum of irrigation and precipitation amounts was provided for each growing season.

Assessment of Spatial Similarity among TSC Maps
To assess the degree of spatial pattern similarity among the three TSC maps, the V-Measure of spatial association between regionalizations [100] was applied. A regionalization is a segmentation of the map into a set of single-connected units. The maps may have areas of disjointed parts. The V-measure is an entropy-based measure that defines the distance from a perfect cluster-matching of two thematic maps, which relies on the homogeneity and completeness measures of two discrete maps. A clustering outcome satisfies homogeneity if all cluster domains contain only data points that are members of a single class. Completeness occurs when all the data pixels that are members of a given group are elements of the same cluster. Homogeneity and completeness are symmetrical measures that compose the V-Measure when averaging them using a harmonic mean [100,101]. Both homogeneity and completeness have values between 0 and 1, resulting in the same scale for the V-Measure. This assessment was applied for all pairs of TSC maps, resulting in three different assessments of spatial associations, using the "sabre" R package [100].

Influence of Field Attributes on ET Spatial Clustering Patterns
To comply with specific objective (3) and assess the rank of impact of different environmental variables within the vineyard on ET levels on different time scales, a random forest (RF) classifier was used. RF is based on a combination of tree classifiers, each generated using a random vector that is sampled independently from the input vector. Each tree casts a unit vote for the most popular class to classify an input vector [102]. Based on studies conducted for the past two decades, RF yields reliable classifications using predictions from an ensemble of decision trees and may be used for ranking the variables with the highest ability to partition the target classes. It is a non-parametric supervised model which creates trees by drawing a subset of training samples through replacement [103]. In this study, the environmental variables detailed in Section 2.2 ( Figure 3) were evaluated for their relative importance to each ET-TSC map. The data were randomly partitioned with no iterations into training and testing using a ratio of 70% and 30%, respectively. The RF algorithm was applied on the training set, after which a prediction over the testing set was performed for evaluation of the model's performance. The results were recorded, including the training error rate and error matrix, and the testing prediction accuracy in % and Kappa coefficient. This procedure was applied using the "RandomForest" package in R [104]. The representation of the results was enabled using the R packages "dplyr" [105], cowplot [106], and "ggplot2" [99].

Effect of ET on Yield
To estimate the extent to which the variability of the yield is explained by the two MZs of each component (trend, seasonality, and remainder) based on the TSC maps and the remainder, we used analysis of variance (ANOVA). ANOVA may be used after meeting a set of assumptions, one of which specifies that the observations should be drawn randomly and independently of each other. In attributes that are spatially autocorrelated, information redundancy is implied, and the duplication of information should be reduced to a minimum [107]. To draw a set of random observations that does not consist of spatial dependence, a minimum distance between observations was defined. Finding the minimum distance was enabled by computing the spatial structure variation with distance using an empirical semivariogram for the yield (Figure 4) [108]. An empirical semivariogram plots half the squared difference (i.e., semivariance) between any two observations against their actual distance, averaged for predefined bins (distance classes). The semivariogram consists of parameters defining its structure, including the range, specifying the maximum distance at which pairs of observations display spatial dependence. The model semivariogram assumes that there is no spatial autocorrelation for distances larger than the range [109]. The model semivariogram was computed using the R package "automap" [110] and the range was used as a minimum distance threshold to extract a random sample of the yield attribute, using the R package "spatialEco" [111]. The selected observations were then used for applying one-way ANOVA with the clusters used as the categorical variable and the yield as the dependent variable. If c > 2, a post-hoc test would be carried out to determine the differences among all pairs of clusters. It should be noted that for c = 2 an unpaired t-test may also be performed. The ANOVA assumptions of normality and homogeneity of variance across groups were tested using the Shapiro-Wilk Normality Test and the Levene's Test, respectively. The procedure was applied using the R packages "raster" [62], "stats" [88], and "car" [112].

Time-Series Clustering
The highest ASW values for all decomposed components were achieved when applying FCM using c = 2, resulting in ASW of 0.53, 0.28, and 0.06 for the trend, seasonality, and remainder, respectively. Consequently, the TSC maps for each of the components were partitioned into two clusters ( Figure 5). The clusters present an overall spatial uniformity with very little fragmentation for all cases. The areal coverage ratio between the clusters for the TSC maps was similar for the trend and seasonality components, with 70% and 30% of the area assigned to clusters 1 and 2, respectively, for the trend component, and 60% and 40% of the area assigned to clusters 1 and 2, respectively, for the seasonality component. The remainder TSC map showed a ratio of 51.5 and 48.5% areal coverage for clusters 1 and 2, respectively.
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 24 carried out to determine the differences among all pairs of clusters. It should be noted that for c = 2 an unpaired t-test may also be performed. The ANOVA assumptions of normality and homogeneity of variance across groups were tested using the Shapiro-Wilk Normality Test and the Levene's Test, respectively. The procedure was applied using the R packages "raster" [62], "stats" [88], and "car" [112].

Time-Series Clustering
The highest ASW values for all decomposed components were achieved when applying FCM using c = 2, resulting in ASW of 0.53, 0.28, and 0.06 for the trend, seasonality, and remainder, respectively. Consequently, the TSC maps for each of the components were partitioned into two clusters ( Figure 5). The clusters present an overall spatial uniformity with very little fragmentation for all cases. The areal coverage ratio between the clusters for the TSC maps was similar for the trend and seasonality components, with 70% and 30% of the area assigned to clusters 1 and 2, respectively, for the trend component, and 60% and 40% of the area assigned to clusters 1 and 2, respectively, for the seasonality component. The remainder TSC map showed a ratio of 51.5 and 48.5% areal coverage for clusters 1 and 2, respectively.  Figure 6 shows corresponding plots to the clusters illustrated in Figure 5, which demonstrate the characteristics and dynamics of the time series of each cluster for each component. Additionally, a third plot shows the difference between the centroids of both clusters to emphasize the dissimilarity between them. For the trend TSC (Figure 6a), the difference plot shows that the highest variation occurred during the 2015 growing season, where Cluster 1 had higher ET values, on average, and the difference ranged between 0.5 mm day −1 in the beginning of the growing season and 0.26 mm day −1 at the end of the season. During 2016 and 2017, the difference was constant around 0.27 mm day −1 , and in 2018 the difference was the lowest, with Cluster 1 having higher values of less than 0.2 mm day −1 than Cluster 2. For the seasonality time series (Figure 6b), the differences throughout the season vary greatly, while the beginning and the end of the growing season were characterized by very low differences. During July and August, Cluster 1 had ET levels of 0.2 mm day −1 higher than Cluster 2. Lastly, the remainder component (Figure 6c) showed high variability with variations of up to 2 mm day −1 . Clusters 1 and 2 showed no differences.  Figure 6 shows corresponding plots to the clusters illustrated in Figure 5, which demonstrate the characteristics and dynamics of the time series of each cluster for each component. Additionally, a third plot shows the difference between the centroids of both clusters to emphasize the dissimilarity between them. For the trend TSC (Figure 6a), the difference plot shows that the highest variation occurred during the 2015 growing season, where Cluster 1 had higher ET values, on average, and the difference ranged between 0.5 mm day −1 in the beginning of the growing season and 0.26 mm day −1 at the end of the season. During 2016 and 2017, the difference was constant around 0.27 mm day −1 , and in 2018 the difference was the lowest, with Cluster 1 having higher values of less than 0.2 mm day −1 than Cluster 2. For the seasonality time series (Figure 6b), the differences throughout the season vary greatly, while the beginning and the end of the growing season were characterized by very low differences. During July and August, Cluster 1 had ET levels of 0.2 mm day −1 higher than Cluster 2. Lastly, the remainder component ( Figure 6c) showed high variability with variations of up to 2 mm day −1 . Clusters 1 and 2 showed no differences.

Irrigation and Precipitation
The amounts of irrigation and precipitation are detailed in Table 1. During 2015 and 2016, the total amount of precipitation (mm) from February to June was much lower than during 2017 and 2018 (about 102, 114, 300, and 227 mm in 2015, 2016, 2017, and 2018, respectively). Consequently, the irrigation inputs to the vineyard were lower during the periods that experienced higher precipitation levels (i.e., 2017 & 2018). The 2015 and 2016 spring rainfall events were very few, while 2015 experienced 14 precipitation events during the spring, which was the lowest number of rainfall days recorded during the research period.

Assessment of Spatial Similarity among TSC Maps
The V-Measure of spatial association between regionalizations was applied to all pairs of TSC maps. The TSC maps of seasonality and trend showed higher similarity in their spatial patterns than each of these TSC maps with the remainder TSC map (Table 2).

Irrigation and Precipitation
The amounts of irrigation and precipitation are detailed in Table 1. During 2015 and 2016, the total amount of precipitation (mm) from February to June was much lower than during 2017 and 2018 (about 102, 114, 300, and 227 mm in 2015, 2016, 2017, and 2018, respectively). Consequently, the irrigation inputs to the vineyard were lower during the periods that experienced higher precipitation levels (i.e., 2017 & 2018). The 2015 and 2016 spring rainfall events were very few, while 2015 experienced 14 precipitation events during the spring, which was the lowest number of rainfall days recorded during the research period.

Assessment of Spatial Similarity among TSC Maps
The V-Measure of spatial association between regionalizations was applied to all pairs of TSC maps. The TSC maps of seasonality and trend showed higher similarity in their spatial patterns than each of these TSC maps with the remainder TSC map (Table 2).

Influence of Field Attributes on ET Spatial Clustering Patterns
The RF algorithm was applied three times and included the dataset of the six environmental variables (Figure 3), and each one of the three TSC maps. The application of RF enabled the generation of a classification model where the cluster map for each component served as the dependent variable and the six covariates received a rank of their relative importance to the clustering maps. Overall, the performance levels for the three RF models were high (Table 3), with overall prediction accuracy >86% based on the error-matrix for the test set. The model for the trend TSC map had the lowest error rate for the training model, as well as the highest prediction accuracy for the testing model, followed by the seasonality and remainder. The most influencing covariate according to mean decrease in accuracy was soil type for the trend and seasonality, followed by elevation (Figure 7). The northness and slope were ranked third and fourth in their degree of importance. The remainder classification, however, showed strong effects of northness, followed by elevation and soil type.

Influence of Field Attributes on ET Spatial Clustering Patterns
The RF algorithm was applied three times and included the dataset of the six environmental variables (Figure 3), and each one of the three TSC maps. The application of RF enabled the generation of a classification model where the cluster map for each component served as the dependent variable and the six covariates received a rank of their relative importance to the clustering maps. Overall, the performance levels for the three RF models were high (Table 3), with overall prediction accuracy > 86% based on the error-matrix for the test set. The model for the trend TSC map had the lowest error rate for the training model, as well as the highest prediction accuracy for the testing model, followed by the seasonality and remainder.
The most influencing covariate according to mean decrease in accuracy was soil type for the trend and seasonality, followed by elevation (Figure 7). The northness and slope were ranked third and fourth in their degree of importance. The remainder classification, however, showed strong effects of northness, followed by elevation and soil type.

Effect of ET on Yield
The Gaussian model semivariogram was constructed using a range of 55.5 m that was used as the minimum distance threshold for the random sampling of yield and the TSC maps. A total of 65 samples were extracted. The mean yield value for the entire vineyard was 9.94 ton ha −1 , while the sample mean was 9.92 ton ha −1 , with a non-significant difference (an unpaired t-test result of t = −0.04, p = 0.97) that ensured that the sample reliably represented the yield for the entire vineyard. ANOVA was carried out after meeting its assumptions, and the findings (Table 4) show that the clusters for the three TSC maps were able to differentiate the levels of yield in the vineyard with statistically significant differences between cluster 1 and 2 for all components. The differences between the mean yield values for the two clusters were the highest for the seasonality TSC map, as well as the difference in variability of yield levels between the two clusters. For Cluster 1, the SD was 2.51, 2.11, and 2.4 ton ha −1 for the trend, seasonality, and remainders, respectively, and the corresponding values for Cluster 2 were 2.86, 2.76, and 2.58 ton ha −1 .

Discussion
This study focused on quantifying the spatial and temporal patterns of ET in a Pinot noir commercial vineyard in California. The framework follows a holistic approach that considered multiple effects of environmental variables, meteorological factors, and different temporal scales to understand the spatial and temporal patterns of ET in the vineyard and partition their characteristics at multiple time scales according to time-series similarities. Each pixel in the dataset included 980 values, representing daily ET levels during four consecutive growing seasons. To analyze the similarities among the time series of the vineyard pixels, a TSC method was applied to each of three time-series components: long-term trend, seasonality, and remainder, which was regarded as the daily variation of the time series. Time-series analysis techniques are widely applied to ET data, due to their enabling long-term records and spatial coverage. Many studies use various techniques of time-series decomposition of ET for forecasting purposes [87,113]. Some used ET for time-series clustering. [114] found that a fuzzy framework for clustering of ET time-series in space was effective for delineating homogeneous ET o regions. [115] used hierarchical time-series clustering of ET based on images acquired from multiple satellites across an entire agricultural river basin. Combining these two approaches resulted in TSC for each time-series component, and enabled us to further explore the spatial characteristics of each time scale.
Decomposition into three time scales enabled us to distinguish between latent temporal dynamics. The processes that define each component vary with space, time, and underlying effects of meteorological conditions and the spatial variability of different attributes in the vineyard. Two clusters were delineated for each component, resulting in TSC maps ( Figure 5) that may also be viewed as the partitioning of the vineyard into MZs for precision irrigation purposes. The remainder clustering was very poor (ASW = 0.11), which indicates strong similarities among the two clusters. The trend and seasonality TSC maps were found to have a stronger similarity of their spatial patterns than the remainder TSC ( Table 2). The supporting plots of ET values for each cluster (Figure 6) may assist in defining the water stress levels of each MZ, as well as the similarities within each MZ. The trend and seasonality plots (Figure 6a,b) uncover the ET differences between the clusters, and the findings show that Cluster 1 had consistently higher ET values throughout the growing season and at a multi-seasonal time scale. During 2017 the ET levels were the highest across the entire vineyard, possibly due to high precipitation amounts. The difference between the cluster centroids enabled investigation of the effect of the meteorological conditions on each cluster's temporal variability during different time scales. The large difference between the two MZs in the trend component (Figure 6a) may be attributed to seasonal variations in spring precipitation events and irrigation amounts (Table 1). These variations might explain the high trend differences among the MZs during 2015. Although the vineyards had received a similar amount of water in terms of rainfall and irrigation inputs, lower differences between MZs were found when spring precipitation levels and a number of rainfall days were higher (2016,2017). During 2018, the total water input to the vineyard was the lowest, resulting in overall lower ET levels in both MZs, as well as a low difference (~0.2 mmday −1 ). A possible explanation could be that higher levels of water availability in the vineyard during the early stages of the season cause a developed and more efficient system of hydraulic conductivity in the grapevines, which in turn become severely stressed towards the end of the growing season [116]. The performance of the RF models for the three components was satisfying (Table 3) and the results of these models (Figure 7) defined the soil type as the most influential static variable affecting ET clustering for the trend and seasonality components. Spring rainfall events highly influence soil moisture content, while the distribution of water input throughout more rainfall events helps to maintain higher soil moisture levels [117] and correspondingly increase water availability to the plants [118]. Additionally, areas with higher elevation levels in the vineyard were characterized by lower ET levels at long-term time scales and during the growing season. Lower, more homogeneous levels of altitudes are related to higher silt and clay content [119], thus allowing for water to be restored in the soil, prevent leakage downhill through both runoff and underground flow, provide better conditions for higher contents of soil moisture [120], and increase the availability of water to the grapevines.
The plots showing the differences in ET levels of the seasonality component ( Figure 6b) enabled exploration of the seasonal dynamics of variation in space. During the beginning of the growing season when the canopy area was small, the differences in ET levels were negligible. In early May, the discrepancies increased, reaching a gap of 0.2 mm day −1 during July. During this phenological stage (bunch closure until veraison), the canopy cover stabilized. After veraison, ET levels were constantly high until the plant suffered from leaf senescence and physiological deterioration, reducing the ability of the stomata to mediate photosynthetic activity [121]. During this phenological stage, the differences were more enhanced due to variation in soil types of each MZ. The vines composing Cluster 1 are located on top of the Kimball silt loam soil group. These are deep soils with moderate water holding capacity. The vines in Cluster 2 were planted on top of shallower soils with low water holding capacity that do not sustain water availability to the plant as efficiently as the Kimball silt loam. Therefore, during June and July when temperature levels are highest, the differences between the two clusters were more pronounced. These findings correspond to previous studies showing that decreasing levels of soil saturation are characterized by higher spatial variability [122] and that higher spatial variability in ET levels was found during the summer season [123].
All components were subjected to impact by the northness, which was the strongest contributor for the remainder. North-and northwest-facing slopes were associated with higher ET levels. These findings correspond with previous studies that demonstrated that south-facing slopes are subjected to stronger radiation [124], therefore lowering the soil moisture content and affecting water availability to the plant and plant density [59,125]. Denser vegetation then transpires more, leading to higher ET values.
Overall, the remainder levels of each pixel were unique and did not enable an accurate differentiation between pixels. The remainder differences were negligible and spanned around zero (Figure 6c). Decomposition of crop ET (i.e., ET c ) in a Mediterranean vineyard also found that fluctuations in ET values are commonly stronger in early spring ( Figure A1) due to highly variable meteorological conditions [87], affecting the ET levels on a daily scale. The extent of influence of the various environmental factors, both meteorological and spatially varying attributes, is site-specific. Although the relations between different attributes and plant conditions are well-documented, the relative importance of these attributes is unique in each agricultural plot.
ET levels during the season may be associated with yield levels, and the analysis of yield against MZ for each component may support this claim. The ANOVA results show that all TSC maps were able to significantly distinguish between yield values spatially. The difference between the two MZs for each component ranged between 1.95 and 2.73 ton ha −1 .
Overall, the daily scale was not optimal for defining MZs and quantifying their underlying processes and averaged temporal dynamics. Practical use of this framework for irrigation modeling would require relating to ET processes at coarser time scales, such as the entire growing season, or during a multi-seasonal timeframe (long-term trend) for optimal differentiation of the vineyard. Site-specific management based on the spatial variability of water stress has been studied for the past years [126][127][128]. Irrigation protocols may benefit from the proposed technique to quantify water stress, both temporally and spatially. Site-specific management for irrigation purposes may be applied while considering the spatial variability of ET time series, along with a supporting quantification of ET levels in each MZ at multiple time scales. It is worth noting that effective irrigation management in the field requires spatially homogeneous clusters and will therefore benefit from undergoing a smoothness procedure prior to application that will generate non-fragmented clusters in space.
Finally, a similar analysis was conducted using the ratio between ET assessment and reference ET (ETo), which related ET to a living reference crop of grass [129]. This ET-ETo ratio (fRET) represents the crop ET while masking out the meteorological conditions and effects and is referred to as the moisture stress factor [36]. The spatial dynamics of fRET were similar to the ET analysis and therefore were not illustrated in this manuscript; however, fRET should also be examined for irrigation purposes while considering time-space-based irrigation application.

Conclusions
The research presented in this manuscript was aimed at defining spatial and temporal patterns of ET within a vineyard, relying on a data fusion framework based on images from multiple satellite sensors. The findings show that the long-term trend and seasonal patterns of the time series provided a similar spatial distribution of pixel similarity. The area of the vineyard planted on top of a silt loam soil had strong within-group similarities and better conditions to consume more water. Overall, the grapevines had a strong response to the altitude and soil characteristics in the vineyard, which was reflected in the temporal dynamics of their ET levels. This spatial heterogeneity also affected yield patterns across the vineyard, while grapevines with higher ET levels were also generally associated with higher levels of yield. The suggested framework may assist in defining site-specific irrigation protocols that consist of delineating the vineyard into zones with spatial similarities. For each zone, the temporal characteristics of ET may also be available, providing a complete illustration of the vineyard's water stress conditions, both spatially and temporally, at different time scales.
With new, high-resolution satellite imagery options, such as Sentinel-2 and Venus, this analysis may be applied at a higher spatial resolution, increasing the detail level of spatial heterogeneity. Moreover, with a higher number of pixels available for a specific vineyard, the statistical time-series analysis would be stronger.
The application of the framework presented in the present paper is not limited to vineyards. It may be used for other agricultural practices, based on other data sources and different variables that quantify plant conditions or soil attributes with regular time intervals. The scope of the framework could be extended to other fields of research, such as environmental phenomena, epidemiology, spatially dependent economic attributes (e.g., consumption dynamics over space), and more.  Gallo Winery made possible the acquisition and processing of the vineyard field data and the financial support for the GRAPEX project by a grant from the NASA Applied Sciences-Water Resources Program (Grant Award NNH17AE39I) supported the development of the satellite-based ET product used in this analysis. Additionally, the research was partially funded by the European Union Horizon 2020 Research and Innovation Programme under grant agreement no. 871128 "European long-term ecosystem, critical zone and socio-ecological systems research infrastructure PLUS".