Adaptive Estimation of Crop Water Stress in Nectarine and Peach Orchards Using High-Resolution Imagery from an Unmanned Aerial Vehicle ( UAV )

The capability to monitor water status from crops on a regular basis can enhance productivity and water use efficiency. In this paper, high-resolution thermal imagery acquired by an unmanned aerial vehicle (UAV) was used to map plant water stress and its spatial variability, including sectors under full irrigation and deficit irrigation over nectarine and peach orchards at 6.12 cm ground sample distance. The study site was classified into sub-regions based on crop properties, such as cultivars and tree training systems. In order to enhance the accuracy of the mapping, edge extraction and filtering were conducted prior to the probability modelling employed to obtain crop-property-specific (‘adaptive’ hereafter) lower and higher temperature references (Twet and Tdry respectively). Direct measurements of stem water potential (SWP, ψstem) and stomatal conductance (gs) were collected concurrently with UAV remote sensing and used to validate the thermal index as crop biophysical parameters. The adaptive crop water stress index (CWSI) presented a better agreement with both ψstem and gs with determination coefficients (R2) of 0.72 and 0.82, respectively, while the conventional CWSI applied by a single set of hot and cold references resulted in biased estimates with R2 of 0.27 and 0.34, respectively. Using a small number of ground-based measurements of SWP, CWSI was converted to a high-resolution SWP map to visualize spatial distribution of the water status at field scale. The results have important implications for the optimal management of irrigation for crops.


Introduction
According to projected climatic changes, agricultural drought periods over Australia can increase up to 20% by 2030 [1].Climate change effects on soil water balance and plant transpiration can impact crop productivity and quality significantly.Both quality and yield of crops can be enhanced by adequate and timely irrigation based on real-time monitoring of water status without increasing the level of agricultural water used.Such improvements in water use efficiency include changing irrigation frequency to match the crop water requirement that maximizes yield and quality.As an efficient indicator of crop water status, the crop water stress index (CWSI) has been introduced based on the difference between foliage and air temperature [2,3].Jones [4] proposes a reformulated CWSI by the difference between canopy temperature (T canopy ) to the reference temperature or threshold temperature at full transpiration (T wet ), which is normalized by the temperature difference between the non-transpiring leaf temperature (T dry ) and T wet as the upper and lower reference temperatures, respectively.Consequently, estimations of CWSI require the measurement of the canopy temperature targeting the dominant foliage, which conventionally relies on manual or continuous point measurements.However, the ground-based handheld thermography approach is labor-intensive, costly, and impractical when monitoring large areas; furthermore, agricultural fields commonly feature spatially heterogeneous biophysical conditions and spatial variation exists in water stress levels between crops/trees and even within a plant [5].
Since numerous data from various remote sensing platforms have become available, several studies of crop water stress have been conducted using remote sensing images as a major input, aiming to simplify the calculation of CWSI, to replace ground-based measurements, and to make it applicable to larger areas [6][7][8].Furthermore, research on water stress detection with unmanned aerial vehicles (UAV) has been recently conducted, providing a platform to monitor water status at field scale with higher spatial resolution [9][10][11].In addition, research on combining remote sensing data and ground-based data in vineyards has been also carried out to predict water status at a temporal and a spatial basis [12,13].
Many remote-sensing CWSI approaches have focused on an image-based analysis from near-ground or proximal platforms, assessing plant water status at the crop level.These approaches use both visible and thermal infrared images to obtain the level of crop water stress as a quantitative index, where the image co-registration method is applied to delineate crop canopy pixels and to obtain crop temperature separately from soil temperature [14].The canopy-related temperature is derived from the thermal image, in which non-canopy temperature such as soil is masked out by the co-registered visible image based on color segmentation.Despite the advantages in obtaining canopy pixels, the co-registration is often time-consuming.In order to eliminate the co-registration process and consequently to facilitate the estimation of stress index, a histogram approach using only thermal images has been proposed.In the approach, an individual histogram of each thermal image is analyzed to distinguish between canopy-related pixels and non-canopy pixels by empirical or statistical methods [15,16].However, a fundamental problem associated with the histogram-based clustering of canopy pixels and the subsequent calculation of canopy temperatures is the presence of mixed pixels that represent partial coverages of both canopy and soil.Due to the nature of raster imagery, mixed pixels of canopy and non-canopy background such as soil exist along their boundaries.The mixed pixels of thermal imagery can cause significant bias in canopy temperature during midday when they are included in the clusters of canopy pixels.Thus, the mixed pixels should be excluded to reduce errors in the histogram-based analysis when determining canopy-related temperatures.
In addition to the delineation of canopy pixels, the determination of T wet and T dry for CWSI has a critical influence on the accuracy of CWSI values as an indicator of crop water stress.The determination of T wet can be made using several methods, such as: (1) the non-water stress baseline (NWSB) by a linear regression function between T canopy -air temperature (T air ) and vapor pressure deficit (VPD) [3]; (2) the wet artificial reference surface (WARS), acting as crops in abundant water status regarded as the temperature of a fully transpiring crop [17]; and (3) a statistical method based on a histogram analysis of canopy temperature [16,18].T dry can be derived from: (1) ancillary meteorological variables based on isothermal radiation [14,19]; (2) a dry reference leaf coated with petroleum jelly (Vaseline) to prevent transpiration [19]; (3) an empirical method by adding a fixed threshold (such as 5-7 • C) to T air [16,[20][21][22]; and (4) a statistical method based on histogram analysis of the canopy temperature.Rud et al. [16] compared three methods to calculate T wet -from WARS, a method based on energy balance, and a statistical histogram approach-and showed that all CWSI values have similar correlations with actual reference data.Since CWSI based on statistical analysis does not require either any meteorological data for the energy balance method or equipment for WARS, it reduces the time related to field work and the CWSI calculation complexity.The statistical histogram approach, however, would provide a valid range of CWSI values only when representative wet and dry references are available in the scene.For example, if all the crops are under either appropriate water supply or water stressed, the statistical approach alone can yield a biased result in the CWSI calculation.Although the method proposed here is not applicable to an area of similar levels of water stress (uniform temperatures), it still remains a practical method to estimate CWSI across an area where canopies representing different water stress levels exist.In addition, the limitations associated with the representative reference temperature could be alleviated using a minimum number of ground truth data (i.e., stem water potential or stomatal conductance).Since the biased results are shifted and scaled from true results, the general water stress model can be adjusted by a simple relationship between CWSI and the ground truth data.
Most estimations of the upper and lower reference temperatures employ single values of T wet and T dry as a reference.However, the reference value of T wet and T dry can change with field circumstances.In the case of temperate tree crops, orchards may contain several crop cultivars and canopy training systems even under the same irrigation regime.Therefore, applying a single set of T wet and T dry could result in inaccurate and misleading CWSI which do not reflect the actual water status, particularly in the orchards with a combination of multi cultivars and different canopy structures.
Thus, this study proposes a new method to estimate CWSI.The main goal is to develop an automated method to interpret very high resolution data to obtain water stress indices at field scale by using UAV-borne sensing at a near real-time window for early detection of water stress in crops.This method can be used to guide optimal irrigation management and to regulate crop quality and yield.The specific objectives are: (1) to develop an adaptive methodology to obtain T wet and T dry based on feature extraction and probability modeling.The proposed method does not require any artificial references or meteorological data to determine the reference baselines, since it is based on classification of sub-regions which have the same property in an entire field, employing automated thresholding of T wet and T dry for each sub-region and; (2) to provide a practical method which can show a realistic representation of the water stress map of plants at high-resolution, utilizing a small number of ground truth biophysical measurements.

Site Description
The site is a nectarine/peach orchard (0.97 ha, 138 m × 70 m), located near Tatura, Victoria, Australia (36 • 26'08"S, 145 • 16'13"E, 114 m AMSL) and administrated by the Department of Economic Development, Jobs, Transport, and Resources (DEDJTR).The region has a temperate climate with average annual rainfall of approximately 480 mm.Annual average reference crop evapotranspiration is approximately 1190 mm [20].A nectarine (Prunus persica (L.) Batsch cv.Autumn Bright) and a peach (Prunus persica (L.) Batsch cv.August Flame) were planted in the winter of 2013.The experiment was carried out on 22 February 2015.The nectarine and peach were at full leaf-up and late maturity phenological stage (i.e., stage 3 of fruit growth).Both cultivars are spaced at an inter-row distance of 4.5 m and 1.0 m between trees with rows orientated in a north-south direction on a fine, sandy, loamy, Shepparton soil.Trees are trained to two canopy systems: vertical leader and Tatura trellis (Figure 1).For the Tatura trellis trees, the canopy dimension of each leader was approximately 1.85 m high and 0.5 m wide.For vertical leader trees, canopy dimensions were approximately 2.0 m high and 0.8 m wide.

Experiment Design
The study site contained eight experimental plots (Figure 2), where control and deficit treatments were applied as an irrigation regime (Table 1).A drip irrigation system was used along the tree lines spaced at every 0.5 m with an emitter flow rate of 1.6 Lh −1 .Irrigation amount was based on a crop water use (crop evapotranspiration, ETc) model [23] in four control plots and the remainder of orchards.Four deficit plots were imposed by withholding irrigation through valves for five days prior to the UAV sensing and ground-based measurements.Each plot contained at least 20 trees.

Experiment Design
The study site contained eight experimental plots (Figure 2), where control and deficit treatments were applied as an irrigation regime (Table 1).A drip irrigation system was used along the tree lines spaced at every 0.5 m with an emitter flow rate of 1.6 L h −1 .Irrigation amount was based on a crop water use (crop evapotranspiration, ET c ) model [23] in four control plots and the remainder of orchards.Four deficit plots were imposed by withholding irrigation through valves for five days prior to the UAV sensing and ground-based measurements.Each plot contained at least 20 trees.

Experiment Design
The study site contained eight experimental plots (Figure 2), where control and deficit treatments were applied as an irrigation regime (Table 1).A drip irrigation system was used along the tree lines spaced at every 0.5 m with an emitter flow rate of 1.6 Lh −1 .Irrigation amount was based on a crop water use (crop evapotranspiration, ETc) model [23] in four control plots and the remainder of orchards.Four deficit plots were imposed by withholding irrigation through valves for five days prior to the UAV sensing and ground-based measurements.Each plot contained at least 20 trees.

Aerial Images
A multi-rotor unmanned aircraft (S900 manufactured by DJI) was used to carry payload with a thermal infrared (TIR) camera (A65, FLIR Systems, Inc., Wilsonville, OR, USA) to measure canopy temperature.The TIR camera was integrated with an on-board computer and GPS for geo-tagging.The TIR camera captures temperature with a spectral range of 7.5-13 µm, a thermal sensitivity of <0.05 • C at +30 • C, a spatial resolution of 640 × 512 pixels, a focal length of 25 mm, and a field of view of 25 • (H) × 20 • (V).Images were acquired across the entire study site at nadir view angle at midday with a cloudless sky to minimize shadows cast by the crop canopy, and to capture the period of peak ET c and likely symptoms of water stress.The UAV was flying at a speed of 2 ms −1 at an altitude of 90 m above ground level (AGL) with the footprint of 39.2 m × 31.3 m and a ground sample distance (GSD) of 6.12 cm, considering sufficient image overlaps for photogrammetric processing.

Ground Targets
Regular surface patterns of orchards can cause an ill-conditioned image matching between consecutive images due to either lack of key points (interesting points) or mismatched points.Thus, two types of ground target, ground control point (GCP), and ground artificial feature (GAF) were designed for a reliable photogrammetric bundle adjustment with a large number of consecutive images.In particular, the GCP target was made of an aluminum sheet with a size of 0.6 m × 0.6 m.Since aluminum material has a low emissivity, the target appears distinctively as a cold object in thermal imagery.In total, 24 GCPs were distributed in the site (Figure 2) and surveyed by a differential GPS (DGPS) with higher than 3-cm positional accuracy.The GAF targets were made of aluminum-coated cardboard with various shapes (triangle, trapezoid, and rectangle) and were randomly distributed across the site at known locations and configurations in order to enhance image matching quality in consecutive TIR images.

Leaf Temperature and Gas Exchange Measurements
Leaf temperature for the calibration of TIR images was measured on five sunlit leaves and five shaded leaves per plot (n = 10 leaves × 8 plots = 80 observations), using a handheld infrared thermometer (TN410LCE, ZyTemp, Radiant Innovation Inc., HsinChu, Taiwan) within a 40-min time window during the UAV image acquisition.Stomatal conductance (g s , mol m −2 s −1 ) and the transpiration rate (E, mmol m −2 s −1 ) were measured on three mature and fully expanded sunlit leaves, corresponding to the fifth leave from the tip of a main shoot, per plot (n = 3 leaves × 8 plots = 24 observations), using a portable photosynthesis system (LI-6400, LI-COR Inc., Lincoln, NE, USA) during midday (solar noon) at the same time of aerial image acquisition.

Stem Water Potential Measurements
Stem water potential (ψ stem , MPa) was measured on two fully expanded shaded leaves per plot (n = 2 leaves × 8 plots = 16 observations) during midday using a Scholander pressure chamber (Model 3000, Soil Moisture Equipment Co., Santa Barbara, CA, USA).The leaves were selected from branches near the main trunk and were covered with an aluminum foil bag for a minimum of one hour of equilibration time prior to the measurement of ψ stem .

Aerial TIR Image Processing
All TIR images were captured in a 14-bit raw format of signal-based values emitted from the objects.The raw images were converted to 16-bit temperature-based images by using a customized code written in Matlab R2014b (Mathworks Inc., Matick, MA, USA).A one-point calibration method [24] was used to retrieve the adjusted surface temperature using an averaged canopy temperature of a sampled tree, which was measured with the handheld infrared thermometer concurrently with the UAV operation.All consecutive images were processed via aerial image triangulation with the geo-tagged flight log and the GCPs by using a photogrammetric software (PhotoScan, Agisoft LLC, Saint Petersburg, Russia).Digital elevation model (DEM) was generated based on the point cloud, which is a set of matching points between overlapping images.Then, a georeferenced orthomosaic was built using the DEM as the surface parameter in the software.

Feature Extraction
Mapping CWSI requires extraction of canopy pixels to sample canopy temperature values excluding soil and other non-leaf materials.Inclusion of the non-canopy temperature can result in significant errors in the interpretation of CWSI.However, the extracted canopy pixels may contain "edge pixels" that feature mixed canopy-soil temperature values.The edge pixels typically show higher temperature than pure canopy pixels during midday periods.Tree rows were oriented north-south in the study site, and the majority of edges ran in corresponding direction in the orthomosaic image (Figure 3a).An edge detection method, combined Sobel and Canny, was applied to exclude ambiguous mixed pixels in the canopy-soil boundary in Matlab R2014b (Figure 3b).Sobel and Canny are designed to detect the gradient changes maximally to edges in vertical and horizontal directions [25].Other common existing edge methods such as Prewitt and Roberts can be also applicable; however, the combined Sobel and Canny outperformed the other options in detecting tree boundaries along the rows in the orthomosaic image.The detected edges were then dilated for more conservative exclusion of the mixed pixels, which were distributed up to six pixels across the edges, along the direction of thermal gradient (Figure 3c).Finally, the pixels of edges were masked out from the original orthomosaic image (Figure 3d).

Adaptive Crop Water Stress Index (CWSI)
As mentioned above, orchards may have a number of cultivars and tree training systems that have an impact on thermal responses.Applying a single set of T wet and T dry over the entire area may result in non-representative water stress indices due to an inaccurate normalized span of CWSI.For instance, the range of CWSI can be shifted, shrunk, or extended resulting in an inconsistent indicator of the actual crop water status.In this study, an adaptive scheme that utilizes multiple sets of reference T wet and T dry for heterogeneous field conditions is proposed.The proposed method is based on dividing the entire field into sub-regions which can be classified with the same properties (Table 2).

CWSI Calculation
The CWSI algorithm applied in the study was suggested by Jones [4], which can be represented as follows: (1 where Tcanopy is canopy temperature from the aerial TIR image, Twet the temperature of a fully transpiring leaf or lower reference, and Tdry is the temperature of a non-transpiring leaf, also considered as upper reference.

Adaptive and
A specific statistical analysis is proposed to estimate adaptive thresholds of Twet and Tdry based on sub-regions (Figure 4).Firstly, a temperature histogram is generated from TIR image subset of each sub-region.This study assumes that Twet can be taken from the coldest part of the histogram from the TIR image [16], and Tdry is also the temperature of a non-transpiring leaf, which can be derived from the highest part of the histogram.

CWSI Calculation
The CWSI algorithm applied in the study was suggested by Jones [4], which can be represented as follows: where T canopy is canopy temperature from the aerial TIR image, T wet the temperature of a fully transpiring leaf or lower reference, and T dry is the temperature of a non-transpiring leaf, also considered as upper reference.

Adaptive T wet and T dry
A specific statistical analysis is proposed to estimate adaptive thresholds of T wet and T dry based on sub-regions (Figure 4).Firstly, a temperature histogram is generated from TIR image subset of each sub-region.This study assumes that T wet can be taken from the coldest part of the histogram from the TIR image [16], and T dry is also the temperature of a non-transpiring leaf, which can be derived from the highest part of the histogram.
The histograms feature distinctive bimodal and normal density distributions, representing vegetated pixels and soil pixels.The temperature differences of the two features are obvious, corresponding to a range of 30-35 • C at midday.The Gaussian mixture modelling (GMM) is fit to the temperature distribution to objectively cluster canopy/soil pixels and to estimate representative T wet and T dry values for canopy.The GMM assigns each observation to a cluster by maximizing the posterior probability, and is composed of k multivariate normal density components, where k is the number of clusters (e.g., k = 2 for bimodal distribution).The higher-temperature component of Gaussian distribution, representing non-canopy pixels, is eliminated to exclude soil background effects in the Gaussian distribution.Then, an automated threshold of T wet and T dry for each sub-region is determined by the critical values of 99% confidence interval limits from the lower-temperature component of Gaussian distribution.The histograms feature distinctive bimodal and normal density distributions, representing vegetated pixels and soil pixels.The temperature differences of the two features are obvious, corresponding to a range of 30-35 °C at midday.The Gaussian mixture modelling (GMM) is fit to the temperature distribution to objectively cluster canopy/soil pixels and to estimate representative Twet and Tdry values for canopy.The GMM assigns each observation to a cluster by maximizing the posterior probability, and is composed of k multivariate normal density components, where k is the number of clusters (e.g., k = 2 for bimodal distribution).The higher-temperature component of Gaussian distribution, representing non-canopy pixels, is eliminated to exclude soil background effects in the Gaussian distribution.Then, an automated threshold of Twet and Tdry for each sub-region is determined by the critical values of 99% confidence interval limits from the lower-temperature component of Gaussian distribution.

Canopy Temperatures from Sub-Regions
The entire study site, except for the four water deficit plots (Figure 2), was submitted to the same irrigation level based on daily water needs calculated using meteorological data.Canopy temperature from each sub-region, however, was characterized by two major features: cultivar type and canopy structure determined by the tree training system (Figure 5a). Figure 5b-f show the distributions of surface temperature (mainly canopy and soil) from the entire region and sub-regions.The histograms featured distinctive bimodal distribution, representing vegetation (23-38 °C) and soil background (38-57 °C), due to apparent temperature differences up to 34 °C.The temperature distribution of soil represents a mixture of dry grass and bare soil, whereas that of vegetation shows actively growing and transpiring canopies.
Due to the bimodal distribution, the GMM technique with two univariate components was fitted to the histograms, featuring different mean and standard deviation for each sub-region.The outputs of GMM fit to canopy temperatures are presented in Table 3.Both canopy temperatures of N_TT and P_TT were distributed with relatively high standard deviations (SDs) of 2.5 °C and 2.9 °C, respectively, whereas the SDs of N_VL and P_VL were in the range of 1.9 °C and 1.5 °C, respectively.The canopy structure of TT is in Y-shape and wider than VL.Consequently, the structure of TT causes the canopy temperatures to spread wider due to greater differences in the leaf energy balance from different vertical positions and lower leaf area density within the canopies.Another interesting feature was the lower mean value of P_VL (27.7 °C) compared to the mean values of all the other

Canopy Temperatures from Sub-Regions
The entire study site, except for the four water deficit plots (Figure 2), was submitted to the same irrigation level based on daily water needs calculated using meteorological data.Canopy temperature from each sub-region, however, was characterized by two major features: cultivar type and canopy structure determined by the tree training system (Figure 5a).Due to the bimodal distribution, the GMM technique with two univariate components was fitted to the histograms, featuring different mean and standard deviation for each sub-region.The outputs of GMM fit to canopy temperatures are presented in Table 3.Both canopy temperatures of N_TT and P_TT were distributed with relatively high standard deviations (SDs) of 2.5 • C and 2.9 • C, respectively, whereas the SDs of N_VL and P_VL were in the range of 1.9 • C and 1.5 • C, respectively.The canopy structure of TT is in Y-shape and wider than VL.Consequently, the structure of TT causes the canopy temperatures to spread wider due to greater differences in the leaf energy balance from different vertical positions and lower leaf area density within the canopies.Another interesting feature was the lower mean value of P_VL (27.

Mapping Adaptive CWSI
Figure 6a shows the estimated CWSI map based on adaptive of T wet and T dry , showing the water status variability over the orchard.In Figure 6b, details of water deficit and control treatment groups in Nectarine VL are presented.The CWSI values estimated in all four water stressed groups were significantly higher (CWSI = 0.79) than CWSI in the control groups (CWSI = approx.0.37).Figure 6c shows the magnified pixel values of CWSI (6.12 cm GSD).

Mapping Adaptive CWSI
Figure 6a shows the estimated CWSI map based on adaptive of Twet and Tdry, showing the water status variability over the orchard.In Figure 6b, details of water deficit and control treatment groups in Nectarine VL are presented.The CWSI values estimated in all four water stressed groups were significantly higher (CWSI = 0.79) than CWSI in the control groups (CWSI = approx.0.37).Figure 6c shows the magnified pixel values of CWSI (6.12 cm GSD).

Validation of Adaptive CWSI
Midday ψstem and gs measurements were used to provide the relationships with the adaptive CWSI.In addition, a single Twet and Tdry, as used in conventional methods, were derived from the

Validation of Adaptive CWSI
Midday ψ stem and g s measurements were used to provide the relationships with the adaptive CWSI.In addition, a single T wet and T dry , as used in conventional methods, were derived from the VPD formula [23,24] in order to compare them with the proposed adaptive method.The GAF targets were placed next to the eight sampled trees, where ψ stem and g s were measured, in order to identify the tree locations in the orthomosaic image and the CWSI map.CWSI values of the sampled trees were extracted, aggregated, and averaged from the CWSI map.The ψ stem measurements in the sampled trees showed a good negative relationship with the adaptive CWSI (Figure 7a).Specifically, the relationship was significant with determination coefficients (R 2 ) of 0.72.Similarly, a strong relationship was found between g s measurements and the adaptive CWSI with R 2 of 0.82 (Figure 8a).The CWSI estimated using the single reference temperature (baseline), however, exhibited a weaker correlation between ψ stem and g s (Figures 7b and 8b).In particular, the T8 tree showed the largest bias among the sampled trees in the single reference CWSI method.The CWSI value sampled from the T8 was calculated as 0.29, which represents a no-water-stress condition, although T8 was in fact in the water deficit treatment group (of P_VL), and the ground ψ stem measurement of T8 indicated a water deficit condition (−2.91 MPa).Overall, the average canopy temperature of P_VL was lower than the values in the other treatment groups.When examined separately, the mean temperature of T8 (28.5 • C) was higher than the mean value of T7 (26.0 • C) from the control plot in the same crop canopy system (P_VL).However, the mean canopy temperature of T8 was in a similar range with of control-group trees in the other sub-regions.The results imply that the reference temperature values for CWSI need to be adjusted with plant cultivars and/or the tree training system in order to make CWSI a general predictor of the plant water stress (i.e., SWP) across different plants and management types.In other words, the proposed adaptive thresholds of T wet and T dry would be required for the horticultural fields that contain multiple plant cultivars and tree training systems.
identify the tree locations in the orthomosaic image and the CWSI map.CWSI values of the sampled trees were extracted, aggregated, and averaged from the CWSI map.The ψstem measurements in the sampled trees showed a good negative relationship with the adaptive CWSI (Figure 7a).Specifically, the relationship was significant with determination coefficients (R 2 ) of 0.72.Similarly, a strong relationship was found between gs measurements and the adaptive CWSI with R 2 of 0.82 (Figure 8a).The CWSI estimated using the single reference temperature (baseline), however, exhibited a weaker correlation between ψstem and gs (Figure 7b and Figure 8b).In particular, the T8 tree showed the largest bias among the sampled trees in the single reference CWSI method.The CWSI value sampled from the T8 was calculated as 0.29, which represents a no-water-stress condition, although T8 was in fact in the water deficit treatment group (of P_VL), and the ground ψstem measurement of T8 indicated a water deficit condition (−2.91 MPa).Overall, the average canopy temperature of P_VL was lower than the values in the other treatment groups.When examined separately, the mean temperature of T8 (28.5 °C) was higher than the mean value of T7 (26.0 °C) from the control plot in the same crop canopy system (P_VL).However, the mean canopy temperature of T8 was in a similar range with those of control-group trees in the other sub-regions.The results imply that the reference temperature values for CWSI need to be adjusted with plant cultivars and/or the tree training system in order to make CWSI a general predictor of the plant water stress (i.e., SWP) across different plants and management types.In other words, the proposed adaptive thresholds of Twet and Tdry would be required for the horticultural fields that contain multiple plant cultivars and tree training systems.sampled trees were extracted, aggregated, and averaged from the CWSI map.The ψstem measurements in the sampled trees showed a good negative relationship with the adaptive CWSI (Figure 7a).Specifically, the relationship was significant with determination coefficients (R 2 ) of 0.72.Similarly, a strong relationship was found between gs measurements and the adaptive CWSI with R 2 of 0.82 (Figure 8a).The CWSI estimated using the single reference temperature (baseline), however, exhibited a weaker correlation between ψstem and gs (Figure 7b and Figure 8b).In particular, the T8 tree showed the largest bias among the sampled trees in the single reference CWSI method.The CWSI value sampled from the T8 was calculated as 0.29, which represents a no-water-stress condition, although T8 was in fact in the water deficit treatment group (of P_VL), and the ground ψstem measurement of T8 indicated a water deficit condition (−2.91 MPa).Overall, the average canopy temperature of P_VL was lower than the values in the other treatment groups.When examined separately, the mean temperature of T8 (28.5 °C) was higher than the mean value of T7 (26.0 °C) from the control plot in the same crop canopy system (P_VL).However, the mean canopy temperature of T8 was in a similar range with those of control-group trees in the other sub-regions.The results imply that the reference temperature values for CWSI need to be adjusted with plant cultivars and/or the tree training system in order to make CWSI a general predictor of the plant water stress (i.e., SWP) across different plants and management types.In other words, the proposed adaptive thresholds of Twet and Tdry would be required for the horticultural fields that contain multiple plant cultivars and tree training systems.

Discussion
TIR imagery can be used to extend a small number of the ground measurements to a spatially distributed measure of crop water status at the field scale effectively.However, delineating plant canopy pixels that are not influenced by the soil background and sampling reference temperature values automatically are not trivial tasks.In this work, an edge detection and elimination approach is employed to exclude ambiguous mixture pixels, and the Gaussian mixture model (GMM) is applied to automatically sample the reference canopy temperature.The lower-temperature GMM parameters are then used to estimate the reference (cold and hot) temperature values at each sub-region.The adaptive thresholds proposed are obtained at 0.5 and 99.5 percentiles of temperature distributions from TIR imagery to specify the T wet and T dry .According to the GMM results, the variability of canopy temperatures between the different combinations of cultivar and the tree architecture was evident.Characterizing the whole study site with a single set of reference temperature values leads to biased and non-representative CWSI results in case the field is composed of multiple combinations of cultivar and tree architecture.The adaptive CWSI proposed can automatically estimate the proxy plant water stress status consistently with the SWP, standardizing temperature distributions into an integrated span of CWSI.
The proposed method is carried out based on the assumption that there a wide range of water stress levels exists in the field, resulting in representative canopy temperature values for stressed and non-stressed plants.This assumption can be justified by assuming that, even under the same irrigation regime, temperature variations caused by different water availabilities for individual crops could occur by external factors such as soil properties, terrain elevation, and uneven irrigation management.However, deviation from the assumed condition can limit the application of the method used in this work.
The proposed method of CWSI poses a few limitations when it is to be applied to all possible conditions of water status in crop fields.A hypothetical example is when the distribution of the field canopy temperature features a narrow range, indicating that the actual crop water status of all crops is at a similar level.In such a case, the estimated CWSI would still arbitrarily allocate a range of water stress levels within the narrow range of the actual crop water status.This is in fact a fundamental limitation of the proposed CWSI method.However, the drifting CWSI range can be anchored in advance with the help of auxiliary data, such as the air temperature.T dry can be predicted by T air + threshold temperature (e.g., 5 • C) in the empirical method [16,[20][21][22].If the estimated T dry in GMM distribution is lower than T dry predicted by T air , it would indicate non-severe crop water stress.On the contrary, if the range of canopy temperature distribution is very narrow and close to T dry predicted by T air , it would mean that most crops are under a water deficit condition.
In addition, remotely estimated CWSI by itself is not a direct measure of actual crop water status, since its mean level and dynamic range can vary with plant cultivars and environmental factors.Nevertheless, the estimated CWSI can be converted to water stress measures, such as ψ stem and g s , when a small number of the measurements are taken concurrently with the TIR imagery.Moreover, SWP (ψ stem ) has been widely accepted as an indicator of crop water status [26].Thus, many crop growers and physiological researchers have determined water status of crops by ψ stem using the pressure chamber method, particularly for grapevines and stone fruits [27,28].However, even though the SWP-based water stress estimation is a reliable method, it poses limitations due to the resources (labour, equipment) required, and hence limits spatial coverage when applied at field scale.However, the SWP-based water stress mapping can be extended to cover large areas when it is combined with the adaptive CWSI map via the CWSI vs. ψ stem relationship.The SWP map in Figure 9 is created by utilizing the CWSI vs. ψ stem model shown in Figure 7a.The map shows the spatially explicit and detailed map of the SWP as a quantitative indicator for crop water stress.
This study was carried out with one flight campaign in the late summer before leaf senescence.A single flight campaign is not sufficient to explain water stress variations at different crop phenological stages.However, canopy cover of stonefruit crops in terms of fractional radiation interception or leaf area index (LAI) does not change much after full leaf-up (late spring) until commencement of leaf senescence (late autumn).The rapid change in cover occurs early in spring, when evaporative demand is low, and the soil water profile is typically full after winter rainfall in Victoria, Australia.The likelihood of crop water stress during the spring period is therefore low, with growers having the option to irrigate if dry soil conditions prevail.Nevertheless, further experiments are necessary to explore crop water stress at different phenological stages.This study was carried out with one flight campaign in the late summer before leaf senescence.A single flight campaign is not sufficient to explain water stress variations at different crop phenological stages.However, canopy cover of stonefruit crops in terms of fractional radiation interception or leaf area index (LAI) does not change much after full leaf-up (late spring) until commencement of leaf senescence (late autumn).The rapid change in cover occurs early in spring, when evaporative demand is low, and the soil water profile is typically full after winter rainfall in Victoria, Australia.The likelihood of crop water stress during the spring period is therefore low, with growers having the option to irrigate if dry soil conditions prevail.Nevertheless, further experiments are necessary to explore crop water stress at different phenological stages.

Conclusions
This study proposed for the first time an adaptive CWSI calculation method using two main approaches: (1) pure canopy extraction by edge detection and statistical analysis of the distribution of surface temperatures; and (2) the adaptive and systematic determination of lower (wet) and upper (dry) references, obtaining 99% CI in the distribution of canopy temperatures in each sub-region.A strong linear relationship between the adaptive CWSI and the ψstem (and gs) was obtained in this work.The present approach can potentially provide a feasible method of assessing the real water stress of plants with high spatial resolution to assess effectively its variability at the plant and field scale for precision irrigation purposes.As a future work, further research on various crop fields and different crop phenological stages will be investigated to make the present method applicable to general cases.

Conclusions
This study proposed for the first time an adaptive CWSI calculation method using two main approaches: (1) pure canopy extraction by edge detection and statistical analysis of the distribution of surface temperatures; and (2) the adaptive and systematic determination of lower (wet) and upper (dry) references, obtaining 99% CI in the distribution of canopy temperatures in each sub-region.A strong linear relationship between the adaptive CWSI and the ψ stem (and g s ) was obtained in this work.The present approach can potentially provide a feasible method of assessing the real water stress of plants with high spatial resolution to assess effectively its variability at the plant and field scale for precision irrigation purposes.As a future work, further research on various crop fields and different crop phenological stages will be investigated to make the present method applicable to general cases.

Figure 1 .
Figure 1.Canopy systems of the nectarine/peach trees in the study site: (a) vertical leader; (b) Tatura Trellis.

Figure 2 .
Figure 2. (a) Location of the study site, Stonefruit Field Laboratory orchard (Tatura, Victoria.Australia); (b) water deficit plots are presented in yellow, and control plots in blue.T1-T8 represents each irrigation treatment and red dots represent the biophysical sampled tree from each plot.Ground control point (GCP) targets in a black square with a white cross.

Figure 1 .
Figure 1.Canopy systems of the nectarine/peach trees in the study site: (a) vertical leader; (b) Tatura Trellis.

Figure 1 .
Figure 1.Canopy systems of the nectarine/peach trees in the study site: (a) vertical leader; (b) Tatura Trellis.

Figure 2 .
Figure 2. (a) Location of the study site, Stonefruit Field Laboratory orchard (Tatura, Victoria.Australia); (b) water deficit plots are presented in yellow, and control plots in blue.T1-T8 represents each irrigation treatment and red dots represent the biophysical sampled tree from each plot.Ground control point (GCP) targets in a black square with a white cross.

Figure 2 .
Figure 2. (a) Location of the study site, Stonefruit Field Laboratory orchard (Tatura, Victoria.Australia); (b) water deficit plots are presented in yellow, and control plots in blue.T1-T8 represents each irrigation treatment and red dots represent the biophysical sampled tree from each plot.Ground control point (GCP) targets in a black square with a white cross.

15 Figure 4 .
Figure 4. Flow chart of the decision of adaptive reference temperatures based on sub-regions (SR).

Figure 4 .
Figure 4. Flow chart of the decision of adaptive reference temperatures based on sub-regions (SR).
Figure 5b-f show the distributions of surface temperature (mainly canopy and soil) from the entire region and sub-regions.The histograms featured distinctive bimodal distribution, representing vegetation (23-38 • C) and soil background (38-57 • C), due to apparent temperature differences up to 34 • C. The temperature distribution of soil represents a mixture of dry grass and bare soil, whereas that of vegetation shows actively growing and transpiring canopies.
7 • C) compared to the mean values of all the other sub-regions (approximately 30 • C), with a narrower spread in the temperature distribution.The different mean and SD values for the histograms indicate that the distribution of canopy temperatures can vary between different cultivar and tree architecture, even under similar irrigation levels.Thus, a concept of adaptive reference baselines determined by region-specific probability modelling is proposed in this study to standardize the different temperature features for CWSI based on the cultivars and the tree architecture types.

Figure 5 .
Figure 5. (a) The orthomosaic image of surface temperature (°C) derived from thermal imagery using Unmanned Aerial Vehicle (UAV) and classified sub-regions based on cultivars and canopy architectures; (b) mixture modelling in the entire site; (c) in the Nectarine trained to Vertical Leader (VL); (d) in the Nectarine trained to Tatura Trellis (TT); (e) in the Peach trained to TT; and (f) in the Peach trained to VL.

Figure 5 .
Figure 5. (a) The orthomosaic image of surface temperature ( • C) derived from thermal imagery using Unmanned Aerial Vehicle (UAV) and classified sub-regions based on cultivars and canopy architectures; (b) mixture modelling in the entire site; (c) in the Nectarine trained to Vertical Leader (VL); (d) in the Nectarine trained to Tatura Trellis (TT); (e) in the Peach trained to TT; and (f) in the Peach trained to VL.

Figure 6 .
Figure 6.(a) Adaptive Crop Water Stress Index (CWSI) map derived from UAV remote sensing using thermal infrared images.Four dotted rectangles in red represent the area of deficit irrigation treatments and four dotted rectangles in blue the full irrigation of control treatments; (b) example of CWSI map depicting water deficit and control treatment groups for Nectarine VL (Vertical Leader); (c) examples of the pixel-level resolution of CWSI for T1 (in deficit) and T2 (in control) sampled trees.

Figure 6 .
Figure 6.(a) Adaptive Crop Water Stress Index (CWSI) map derived from UAV remote sensing using thermal infrared images.Four dotted rectangles in red represent the area of deficit irrigation treatments and four dotted rectangles in blue the full irrigation of control treatments; (b) example of CWSI map depicting water deficit and control treatment groups for Nectarine VL (Vertical Leader); (c) examples of the pixel-level resolution of CWSI for T1 (in deficit) and T2 (in control) sampled trees.

Figure 7 .
Figure 7. (a) Relationship between stem water potential (ψ stem ) and adaptive CWSI; (b) relationship between ψ stem and CWSI of single reference baseline.

Figure 8 .
Figure 8.(a) Relationships between stomatal conductance (gs) and adaptive CWSI; (b) relationship between gs and CWSI of single reference baseline.Figure 8. (a) Relationships between stomatal conductance (g s ) and adaptive CWSI; (b) relationship between g s and CWSI of single reference baseline.

Acknowledgments:
This research was supported by the Innovation Seed Fund for Horticulture Development grant from the University of Melbourne and Department of Economic Development, Jobs, Transport, and

Table 1 .
Descriptions of sampled trees based on cultivar, training system, and irrigation treatment.

Table 2 .
Classified sub-regions within the studied field based on crop cultivar and tree training system.

Table 2 .
Classified sub-regions within the studied field based on crop cultivar and tree training system.

Table 3 .
Thermal values of canopy temperature using probability modeling.

Table 3 .
Thermal values of canopy temperature using probability modeling.