Automatic Coregistration Algorithm to Remove Canopy Shaded Pixels in UAV-Borne Thermal Images to Improve the Estimation of Crop Water Stress Index of a Drip-Irrigated Cabernet Sauvignon Vineyard

Water stress caused by water scarcity has a negative impact on the wine industry. Several strategies have been implemented for optimizing water application in vineyards. In this regard, midday stem water potential (SWP) and thermal infrared (TIR) imaging for crop water stress index (CWSI) have been used to assess plant water stress on a vine-by-vine basis without considering the spatial variability. Unmanned Aerial Vehicle (UAV)-borne TIR images are used to assess the canopy temperature variability within vineyards that can be related to the vine water status. Nevertheless, when aerial TIR images are captured over canopy, internal shadow canopy pixels cannot be detected, leading to mixed information that negatively impacts the relationship between CWSI and SWP. This study proposes a methodology for automatic coregistration of thermal and multispectral images (ranging between 490 and 900 nm) obtained from a UAV to remove shadow canopy pixels using a modified scale invariant feature transformation (SIFT) computer vision algorithm and Kmeans++ clustering. Our results indicate that our proposed methodology improves the relationship between CWSI and SWP when shadow canopy pixels are removed from a drip-irrigated Cabernet Sauvignon vineyard. In particular, the coefficient of determination (R2) increased from 0.64 to 0.77. In addition, values of the root mean square error (RMSE) and standard error (SE) decreased from 0.2 to 0.1 MPa and 0.24 to 0.16 MPa, respectively. Finally, this study shows that the negative effect of shadow canopy pixels was higher in those vines with water stress compared with well-watered vines.


Introduction
Water availability is a critical limiting factor in the agricultural industry; therefore, a wide range of new technologies and strategies have been adopted to optimize the agricultural water consumption [1][2][3][4]. Granier et al. [5] argued that the measurements of physiological parameters can provide better information about the whole-plant-level water use with changing atmospheric water demands. For example, the water potential has been used to characterize the plant water stress and to schedule irrigation in vineyards [6][7][8], as well as for nuts trees [9,10], and olive trees [11,12]. However, water potential is typically measured on a plant-by-plant basis leading to high costs and requiring a considerable time when these measurements are extended to cover a large area [13,14]. This limitation has motivated the development of cost-and time-effective alternatives to evaluate plant water status.
Multispectral imagery to capture images at the leaf and canopy levels has been proposed as an effective tool for agricultural applications [15] to indirectly and remotely assess plant water status. For example, Rapaport et al. [16] reported that estimating the water balance index (WABI-2) using visible (538 nm) and short-wave infrared (1500 nm) spectrum is a good indicator of water stress in grapevines. Rallo, et al. [17] suggested that spectral information between the near infrared (NIR) (750 nm) and short-wave infrared (SWIR) (1550 nm) ranges can improve the prediction of leaf water potential. In addition, Pôças et al. in [18,19] showed that the wavelength information of visible (VIS) and NIR spectra can be used to predict water status. Poblete, et al. [20] suggested that artificial neural networks using information obtained from 500 to 800 nm could be used to predict the stem water potential (SWP) spatial variability in vineyards.
Furthermore, the Crop Water Stress Index (CWSI) derived from the radiometric temperature of a plant canopy measured using thermal infrared (TIR) sensors has been suggested as a reliable tool to assess water stress [21][22][23][24][25] showing good correlations with ground measurements of water potential. However, as in the case of ground-based water potential measurements, when large crop areas are to be assessed, the ground-based TIR measurements can still be time-consuming and impractical. Thus, remotely collected TIR imagery has been suggested as an alternative tool that can provide crop status information over large regions in a non-invasive manner [26][27][28][29]. In particular, unmanned aerial vehicles (UAV) have become a useful remote sensing tool, having significant advantages in terms of cost, versatility, and high spatial resolution [30]. The CWSI studies using UAV-borne sensors have achieved a high correlation with the plant water status measured using ground-based measurements [14,29,[31][32][33].
However, the UAV-borne TIR sensing for plant water stress suffers from the technical issue of the potential degradation of the canopy temperature information by the pixels of a shaded (or shadow) canopy; this is because the surface temperature of sunlit canopy is known to better represent the plant water stress. Existing methods to remove these shaded pixels from remote-sensing images can be divided into two principal steps: shadow detection followed by a de-shadow process [34]. The first step, shadow detection, can be conducted by either thresholding or modeling [35]. The thresholding process is more common as it is less complicated than modeling, because modeling requires prior information of shadows and mathematical conceptualization; consequently, modeling is applied only to specific cases. The thresholding process involves finding the optimal threshold value of a digital number based on histograms to segregate shadow information from other types of information. Previous studies have used different wavelengths to elucidate thresholds for shadow deletion. For example, NIR (757-853 nm) [36], the ratio between blue (450-520 nm) and NIR (760-900 nm) [37], Infrared (10.4-12.5 µm) [38], and indices [39][40][41] have been used to separate undesired information. However, the TIR information obtained by the commonly used thermal imaging devices (based on an uncooled microbolometer) does not provide sufficient sensitivity for subtle temperature variation [15]; therefore, this method often fails to distinguish shadow canopy pixels from shadow soil pixels. Considering the issues with the shadow canopy pixels, an important process in thermal image processing is shadow pixel removal to improve the resampling of the sunlit canopy information [42]. Zarco-Tejada et al. [43] and Suárez et al. [44] highlight the importance of resampling sunlit canopy pixels using hyperspectral and multispectral imagery, respectively, to assess the plant water stress. Using UAV-borne thermal imagery, several studies have proposed different methodologies to achieve shadow removal and avoid the shadow effect in the case of thermal images. For example, Zarco-Tejada et al. [45] suggested that only the center portion of the canopy row be sampled to minimize the inclusion of shadow canopy pixels. Gonzalez-Dugo et al. [46] sampled the central 50% of the crown pixels of the canopy. Santesteban et al. [29] detailed the complexity of avoiding shadow information, especially in thermal imageries, and proposed a Digital Elevation Model (DEM) and Otsu [47] combined methodology to filter shadows using height differences presented in the ground. Despite the proposed methodologies of shadow removal for UAV-borne images and, even if capturing images in overcast conditions can minimize the intensity of shadowing [48], the identification of shadow canopy pixels produced by the canopy over itself as information to be deleted in thermal images is not considered. For a drip-irrigated vineyard, Figure 1 shows an example of a thermal image in which the shadow canopy pixels cannot be identified on comparing with the visible imagery VIS (490 nm). avoiding shadow information, especially in thermal imageries, and proposed a Digital Elevation Model (DEM) and Otsu [47] combined methodology to filter shadows using height differences presented in the ground. Despite the proposed methodologies of shadow removal for UAV-borne images and, even if capturing images in overcast conditions can minimize the intensity of shadowing [48], the identification of shadow canopy pixels produced by the canopy over itself as information to be deleted in thermal images is not considered. For a drip-irrigated vineyard, Figure 1 shows an example of a thermal image in which the shadow canopy pixels cannot be identified on comparing with the visible imagery VIS (490 nm). Considering the effect of shadow on water stress estimation, it is crucial to determine shaded pixels and remove them [49][50][51]. Möller et al. [23] proposed a methodology to detect grapevine crop water status using thermal and visible images collected using truck-mount sensors at 15 m above the ground to sample the sunlit canopy information; they used Ground Control Points (GCP) made of cross-marked aluminum plates to geo-reference, align, and coregister the images from two different sensors. Leinonen and Jones [42] also proposed a methodology to assess water stress in grapevine and broad bean fields using ground-obtained thermal and visible images; their methodology was based on non-automatic (by expert user) selection of GCP to overlay the images and later warp and resample the images to obtain the sunlit canopy information. Finally, Smith et al. [52] proposed a methodology to detect regions of soil moisture deficit from a spinach plantation using thermal and visible images. Bulanon et al. [53] proposed a methodology for fruit detection using thermal and VIS imagery in which four corners of a ground-marked region of interest were used to coregister VIS and TIR images and perform shadow removal. However, in all of these studies, challenges in coregistering optical and TIR images were reported when the images were combined for shadow removal [54] using non-automatic coregistration. Considering this, our study proposes an automatic scheme based on Scale Invariant Feature Transformation (SIFT) computer vision algorithm and an improved matching pairs point selection to remove shaded pixels in a UAV-borne thermal image to improve the estimation of the CWSI for a drip-irrigated Cabernet Sauvignon vineyard grown under Mediterranean climate conditions.

Site Description and Experimental Design
The study site has a predominant typical Mediterranean climate with a summer period from December to March that is usually dry (2.2% of annual rainfall) and hot with an average daily temperature of 21 °C, and spring that is usually wet (16% of annual rainfall). Average annual rainfall in the region is about 500 mm, which falls primarily during April to August. Considering the effect of shadow on water stress estimation, it is crucial to determine shaded pixels and remove them [49][50][51]. Möller et al. [23] proposed a methodology to detect grapevine crop water status using thermal and visible images collected using truck-mount sensors at 15 m above the ground to sample the sunlit canopy information; they used Ground Control Points (GCP) made of cross-marked aluminum plates to geo-reference, align, and coregister the images from two different sensors. Leinonen and Jones [42] also proposed a methodology to assess water stress in grapevine and broad bean fields using ground-obtained thermal and visible images; their methodology was based on non-automatic (by expert user) selection of GCP to overlay the images and later warp and resample the images to obtain the sunlit canopy information. Finally, Smith et al. [52] proposed a methodology to detect regions of soil moisture deficit from a spinach plantation using thermal and visible images. Bulanon et al. [53] proposed a methodology for fruit detection using thermal and VIS imagery in which four corners of a ground-marked region of interest were used to coregister VIS and TIR images and perform shadow removal. However, in all of these studies, challenges in coregistering optical and TIR images were reported when the images were combined for shadow removal [54] using non-automatic coregistration. Considering this, our study proposes an automatic scheme based on Scale Invariant Feature Transformation (SIFT) computer vision algorithm and an improved matching pairs point selection to remove shaded pixels in a UAV-borne thermal image to improve the estimation of the CWSI for a drip-irrigated Cabernet Sauvignon vineyard grown under Mediterranean climate conditions.

Site Description and Experimental Design
The study site has a predominant typical Mediterranean climate with a summer period from December to March that is usually dry (2.2% of annual rainfall) and hot with an average daily temperature of 21 • C, and spring that is usually wet (16% of annual rainfall). Average annual rainfall in the region is about 500 mm, which falls primarily during April to August.
Flight campaigns and climate measurements were carried out in a drip-irrigated Cabernet Sauvignon vineyard located in the Pencahue Valley, Maule Region, Chile (35 • 20 L.S; 71 • 46 L.W). The three-year-old wine grapes were trained on a vertical shoot positioned (VSP) system. The vineyard fractional cover, which represents the dimensionless parameter of ground covered vegetation over uncovered ground [55], was 19%. In addition, the vineyard with east-west oriented rows (at 1 m × 2 m) Sensors 2018, 18, 397 4 of 17 was irrigated daily using 2 L·h −1 drippers spaced at intervals of 1 m. The soil is las doscientas type with a compact arsenic soil texture with high levels of Fe and Mn.
The experimental design consisted of two completely randomized treatments (well-watered and deficit-irrigated vines) with four replications (six vines per replication). The SWP for well-watered vines showed values that ranged between −0.6 and −0.8 MPa and the deficit-irrigated vines showed values that ranged between −0.9 and −1.25 MPa. The SWP was measured at the time of UAV overflight [56] using a pressure chamber (PMS 600, PMS Instrument Company, Corvallis, OR, USA) from the middle vines for each repetition. A total of 32 leaves from the middle zone of the canopy were measured corresponding to two mature and healthy sun-exposed leaves that were previously covered with plastic bags and coated with aluminum foil for at least 1 h before measurements [6].

Cameras and Image Processing Description
A multispectral camera was used to collect VIS-NIR images for shadow identification. The images were obtained from a Micro MCA-6 camera (Tetracam's Micro Camera Array), which has an array of sensors with band-path filters whose center-wavelengths are 490, 550, 680, 720, 800, and 900 nm with a resolution of 1280 (H) × 1024 (V). For thermal infrared imaging, the FLIR TAU2 640 (FLIR Systems, Inc., Wilsonville, OR, USA) was used. This camera consists of an uncooled microbolometer of 640 (H) × 512 (V) with a pixel pitch of 17 µm and spectral band ranging between 7.5 and 13.5 µm. The thermal calibration was conducted using the methodology proposed by Ribeiro-Gomes et al. [57], in which an artificial neural network is used with the sensor temperature and the digital response of the sensor as input and a Wallis filter to improve the photogrammetry process. Further, the multispectral calibration was performed using the methodology proposed by Poblete et al. [20] in which normalization of the reflectance was performed using a "white reference" Spectralon panel (Labsphere Inc., Sutton, NH, USA) and compared comparison was made with that obtained using a spectroradiometer (SVC HR-1024, Spectra Vista Cooperation, Poughkeepsie, NY, USA) to account for any relative spectral response of each band of the camera as proposed by Laliberte et al. [58].
All images from both sensors were processed using a photogrammetric software PhotoScan (Agisoft LLC, Saint Petersburg, Russia) to stitch the images together to increase the Field of View (FOV) while maintaining the intrinsic characteristics of both cameras [59]; the same software parameters proposed by Ribeiro-Gomes et al. [57] for the same type of sensor were used for stitching.
Finally, the meteorological conditions and flight description on the day of SWP are detailed in Table 1.

CWSI Calculation
The calculation of the CWSI was first proposed by Jones [60] and was described as follows: where T canopy represents the canopy temperature obtained using the UAV-borne TIR. T wet represents the temperature of a fully transpiring canopy and T dry represents the temperature of a non-transpiring fully stressed canopy. As proposed by King et al. [21] and Grant et al. [51], these values do not necessarily need to be an absolute canopy temperature limit value, but serve rather as indicator temperature to scale measured canopy temperatures to the environment for calculating relative water stress. The process for obtaining the values of T dry and T wet was based on the methodology proposed by Park et al. [31]. The process involved using an adaptive approximation based on the TIR histograms derived from the images, and then, identifying the T dry and T wet values after the shadow filtering process by considering the highest and the lowest parts of the histograms, respectively.

Scale Invariant Feature Transformation (SIFT) and Random Sample Consensus (RANSAC)
This algorithm was originally proposed by Lowe [61] to extract characteristic features from images in a robust manner, which is independent of variations in scaling, rotation, translations, and illumination. The algorithm workflow was summarized and explained in detail by Ghosh and Kaabouch [59]. Based on the study by Ghosh and Kaabouch [59], the five primary steps involved in this algorithm are discussed briefly in the following lines. The scale-space construction step is based on applying several Gaussian filters to the image to compute the differences between the adjacent resulting images. Then, in the scale-space extrema detection, a selection of the highest and smallest values between each point and the 26 consecutive neighbors is conducted. Further, in the keypoint localization step, low contrast and edge response points are discarded. For the resulting keypoints, the orientation assignments based on the gradient directions are computed. To define the keypoint descriptors, histograms over each keypoint orientation is calculated considering the highest peak and values under 80% as predominant directions of the local gradients.
After these five steps are performed, the nearest neighbor of a keypoint in the first image is identified from the keypoints of the second image. To remove the outliers and filtering the incorrectly matched points, the RANSAC algorithm is applied. The RANSAC algorithm was first proposed by Fischler and Bolles [62] as a resampling technique for estimating the parameters of a model, using data that may be contaminated by outliers [63]. As suggested by Derpanis [64], the RANSAC algorithm can be summarized in five principal steps: (1) randomly selectf the minimum number of points required to determine the model parameters; (2) solve for the parameters of the model; (3) determine the number of points under a tolerance value; (4) if the ratio of points resulting from the previous step over the total number of points exceeds a predefined threshold, estimate the model with a new set of points; and (5) otherwise, repeat Steps 1-4 (with a maximum of n iterations). Because the value selected for n is high to avoid mismatching, the RANSAC algorithm is time consuming [63] and has a high computational complexity when coupled with the SIFT algorithm [65]. In addition, as RANSAC is a non-deterministic algorithm [66,67], it does not guarantee the return of an optimal solution [68], resulting in different results for different runs [69]. Furthermore, when computed with few SIFT-derived keypoints, it can be sensitive to initial conditions [70]. Considering these issues, and because thermal and visible images have different characteristics with, for most cases, different spatial resolutions, their coregistering process is complex and the assumption of global statistical dependence is not completely satisfied [71]. The RANSAC algorithm between both images leads to different pairing points, which affects the overall performance and consistency in results. This statement is consistent with Turner, et al. [72], who, using RANSAC algorithm, concluded that thermal mosaics showed lower accuracy when coregistered with multispectral images, compared with visible mosaics. To address this issue, we propose an alternative filtering of matching points based on statistical parameters between previous matched pairs. Image analysis and processing were performed using the MATLAB 2017a (Mathworks Inc., Matick, MA, USA) based on the methodology proposed by Vedaldi and Fulkerson [73].

Slope Filtering of Matching Points
As discussed above, in our method, a statistical filtering method was applied to filter previously mismatched points and our results were compared with those obtained using the RANSAC method. Our process involved both images (thermal and multispectral) as a continuous image joined by the resulting matching pair point (Figure 2).

Slope Filtering of Matching Points
As discussed above, in our method, a statistical filtering method was applied to filter previously mismatched points and our results were compared with those obtained using the RANSAC method. Our process involved both images (thermal and multispectral) as a continuous image joined by the resulting matching pair point (Figure 2). The slope of previous matched points was calculated using the Euclidean formula as follows: where ( , ) corresponds to the thermal image descriptor 1 and ( , ) corresponds to the multispectral image descriptor 1. Then, the statistical parameters were calculated for each matched feature and the filtering was conducted based on the mode of the slopes. As an example, in Figure 2, the previous matched descriptor 2 should not be considered because it was identified as a correctly matched feature, but the slope of both descriptors is different from the mode of all the slopes.

Shadow Filtering
As proposed by Shahtahmasseb et al. [34], histogram-based thresholding methods are commonly employed for shadow detection. In this study, with the aim of identifying the optimal wavelength for shadow detection, histograms for 112 UAV-borne images obtained from the vineyard were analyzed. The K-means clustering algorithm with k-means++ was used to optimize the thresholding [74]; this process was applied for shadow detection to six multispectral bands (490, 550, 680, 720, 800, and 900 nm) and their relative performance was compared. Five clusters and 200 iterations were selected for the classification. After performing the previous steps, the classified clusters for shadows from the six multispectral bands were used to build a mask that was applied to their RGB composition to evaluate the accuracy of the classification. Using the abovementioned process and the previously described SIFT algorithm, the resulting mask was coregistered with the thermal images to delete canopy shadow pixels.

Statistical Analysis
For assessing the impact of shadow canopy pixels on the linear correlation between CWSI and SWPl the coefficient of determination (R 2 ) was calculated. In addition, the root mean square error (RMSE) and standard error (SE) parameters were calculated for the comparison. The slope of previous matched points was calculated using the Euclidean formula as follows: where (x 1 , y 1 ) corresponds to the thermal image descriptor 1 and (x 1, y 1 ) corresponds to the multispectral image descriptor 1. Then, the statistical parameters were calculated for each matched feature and the filtering was conducted based on the mode of the slopes. As an example, in Figure 2, the previous matched descriptor 2 should not be considered because it was identified as a correctly matched feature, but the slope of both descriptors is different from the mode of all the slopes.

Shadow Filtering
As proposed by Shahtahmasseb et al. [34], histogram-based thresholding methods are commonly employed for shadow detection. In this study, with the aim of identifying the optimal wavelength for shadow detection, histograms for 112 UAV-borne images obtained from the vineyard were analyzed. The K-means clustering algorithm with k-means++ was used to optimize the thresholding [74]; this process was applied for shadow detection to six multispectral bands (490, 550, 680, 720, 800 and 900 nm) and their relative performance was compared. Five clusters and 200 iterations were selected for the classification. After performing the previous steps, the classified clusters for shadows from the six multispectral bands were used to build a mask that was applied to their RGB composition to evaluate the accuracy of the classification. Using the abovementioned process and the previously described SIFT algorithm, the resulting mask was coregistered with the thermal images to delete canopy shadow pixels.

Statistical Analysis
For assessing the impact of shadow canopy pixels on the linear correlation between CWSI and SWPl the coefficient of determination (R 2 ) was calculated. In addition, the root mean square error (RMSE) and standard error (SE) parameters were calculated for the comparison.

SIFT and Comparison between RANSAC and Slope Filtering for Filtering Matched Features
The comparison between the abovementioned filtering processes was conducted using a complete orthomosaic obtained from the vineyard built using 112 images. The RANSAC algorithm outputs and its fluctuation on the filtered points is shown in Figure 3. Figure

SIFT and Comparison between RANSAC and Slope Filtering for Filtering Matched Features
The comparison between the abovementioned filtering processes was conducted using a complete orthomosaic obtained from the vineyard built using 112 images. The RANSAC algorithm outputs and its fluctuation on the filtered points is shown in Figure 3. Figure   In addition, the statistical parameters of the matching features slope are listed in Table 2. The selected statistical filter was based on the mode of the slope and its result of filtering is shown in Figure 4.  In addition, the statistical parameters of the matching features slope are listed in Table 2. The selected statistical filter was based on the mode of the slope and its result of filtering is shown in Figure 4.

Multispectral Band Selection for Shadow Detection
For shadow identification, the histogram distribution was calculated to detect peaks related to shadow information. An example of one image per band and its distribution is shown in Figure 5 for a drip-irrigated vineyard.
As mentioned previously, five clusters were selected and 200 iterations were conducted for the classification of all the images. The K-means++ methodology [74] was used to set the thresholds in which shadows should be identified for the six bands. After that, each generated mask was applied to an RGB image composition to identify which better represented the shadow. Figure 6 shows the filtered images and the five identified clusters per band. As is clear in Figure 6, for each band, the clustering process allowed the identification of different types of information.

Multispectral Band Selection for Shadow Detection
For shadow identification, the histogram distribution was calculated to detect peaks related to shadow information. An example of one image per band and its distribution is shown in Figure 5 for a drip-irrigated vineyard.

Multispectral Band Selection for Shadow Detection
For shadow identification, the histogram distribution was calculated to detect peaks related to shadow information. An example of one image per band and its distribution is shown in Figure 5 for a drip-irrigated vineyard.
As mentioned previously, five clusters were selected and 200 iterations were conducted for the classification of all the images. The K-means++ methodology [74] was used to set the thresholds in which shadows should be identified for the six bands. After that, each generated mask was applied to an RGB image composition to identify which better represented the shadow. Figure 6 shows the filtered images and the five identified clusters per band. As is clear in Figure 6, for each band, the clustering process allowed the identification of different types of information. As mentioned previously, five clusters were selected and 200 iterations were conducted for the classification of all the images. The K-means++ methodology [74] was used to set the thresholds in which shadows should be identified for the six bands. After that, each generated mask was applied to an RGB image composition to identify which better represented the shadow. Figure 6 shows the filtered images and the five identified clusters per band. As is clear in Figure 6, for each band, the clustering process allowed the identification of different types of information. For the 490-nm and 550-nm group of images, cluster 1 (C1) tends to identify both soil and internal shadows, while cluster 2 (C2) tends to classify vegetation information. On the other hand, for the 680-nm image, shadow is misclassified, nevertheless C1 allows directly identifying vegetation information. Finally, the 720 nm, 800 nm, and 900 nm images seem to misclassify shadow, mixing classified information in both cases with grassy soil and bare soil. To validate our method, a mask was built from the C1-680 nm image to select just vine canopy which included internal shaded canopy pixels. The resulting mask was applied over the images and K-means++ algorithm was carried out to classify vegetation and internal shaded canopy pixels ( Figure 7A). For the 490-nm and 550-nm group of images, cluster 1 (C1) tends to identify both soil and internal shadows, while cluster 2 (C2) tends to classify vegetation information. On the other hand, for the 680-nm image, shadow is misclassified, nevertheless C1 allows directly identifying vegetation information. Finally, the 720 nm, 800 nm, and 900 nm images seem to misclassify shadow, mixing classified information in both cases with grassy soil and bare soil. To validate our method, a mask was built from the C1-680 nm image to select just vine canopy which included internal shaded canopy pixels. The resulting mask was applied over the images and K-means++ algorithm was carried out to classify vegetation and internal shaded canopy pixels ( Figure 7A).
To assess and validate the accuracy of shadow identification, confusion matrices were calculated for the randomly selected marked winegrapes, as shown in Figure 7B, for the six bands to assess the percentage of correct shadow classification. The percentages of well classified shadow for 490, 550, 680, 720, 800 and 900 nm were 90%, 68%, 89%, 77%, 66% and 58%, respectively (Table 3). Cohen's kappa coefficient value, which is used to assess the chance-corrected agreement between two classifications [75], for each band was 0.77, 0.56, 0.76, 0.71, 0.54, and 0.41, respectively.
for the 680-nm image, shadow is misclassified, nevertheless C1 allows directly identifying vegetation information. Finally, the 720 nm, 800 nm, and 900 nm images seem to misclassify shadow, mixing classified information in both cases with grassy soil and bare soil. To validate our method, a mask was built from the C1-680 nm image to select just vine canopy which included internal shaded canopy pixels. The resulting mask was applied over the images and K-means++ algorithm was carried out to classify vegetation and internal shaded canopy pixels ( Figure 7A).  Based on this information, the 490-nm image, which showed the highest percentage of accuracy and Cohen's kappa coefficient value, was selected to be coregistered with the thermal image and for thermal shadow deletion.

Effect of Shadow Removal on the Relationship between CWSI and SWP
To assess the impact of shadow removal on the prediction of the SWP using CWSI, UAV-borne TIR images with and without removal of shadow canopy pixels were compared. Figure 8 shows the thermal image after automatic coregistration and shadow canopy removal. The colored regions correspond to the filtered temperature information, while the background image represents the initial vineyard information without filtered canopy shadow pixels. The mean values of canopy temperature for the cases with and without shadow canopy were 28.84 ± 1.8 • C and 29.95 ± 2.05 • C, respectively. In addition, the relationship between CWSI and SWP is shown in Figure 9. The mean values of CWSI for the non-filtered information were 0.45 ± 0.14, while those for the filtered information were 0.52 ± 0.17. Finally, the results indicated that the relationship between the CWSI and SWP improved after using the automatic coregistration algorithm. In particular, the coefficient of determination (R 2 ) increased from 0.64 to 0.77. In addition, the values of RMSE and SE decreased from 0.2 to 0.1 MPa and 0.24 to 0.16 MPa, respectively.

Discussion
The selection of B1 (490 nm) as the better multispectral band for classifying shadow canopy pixels was consistent with the previous study by Ünsalan et al. [76], who used the k-means and blue information derived from the RGB spectrum to segment information avoiding shadow pixels to extract street networks and detect houses. This band selection was also proposed by Sirmacek et al. [77], who used the blue wavelength spectrum to detect shadows for building detection, suggesting that this region was dominantly better even compared with green and red for shadow pixels identification [78,79]. The selection of 490 nm image was also preferred when compared with upper wavelengths, in which blue spectrum showed better results for shadow detection increasing the performance for near infrared and shortwave infrared [80]. This validates the previous assumption that internal canopy shadow cannot be identified by TIR imagery. Considering this, the importance of coregister thermal and visible images for detecting shadow pixels was also highlighted by Leinonen et al. [42] who using ground cameras with a non-automatic methodology concluded that one of the principal steps is to correct overlapping VIS and TIR images to assess vine water status.
In the present study, the SWP values between the stressed and well-watered vines [81] can be easily identified. The relationship between the CWSI and SWP improved when the shadow pixels were removed from the vine canopy using the suggested automatic algorithms. For the vineyard, the fractional cover was 19%, while the percentage of canopy shaded pixels was 43%. This indicates that

Discussion
The selection of B1 (490 nm) as the better multispectral band for classifying shadow canopy pixels was consistent with the previous study by Ünsalan et al. [76], who used the k-means and blue information derived from the RGB spectrum to segment information avoiding shadow pixels to extract street networks and detect houses. This band selection was also proposed by Sirmacek et al. [77], who used the blue wavelength spectrum to detect shadows for building detection, suggesting that this region was dominantly better even compared with green and red for shadow pixels identification [78,79]. The selection of 490 nm image was also preferred when compared with upper wavelengths, in which blue spectrum showed better results for shadow detection increasing the performance for near infrared and shortwave infrared [80]. This validates the previous assumption that internal canopy shadow cannot be identified by TIR imagery. Considering this, the importance of coregister thermal and visible images for detecting shadow pixels was also highlighted by Leinonen et al. [42] who using ground cameras with a non-automatic methodology concluded that one of the principal steps is to correct overlapping VIS and TIR images to assess vine water status.
In the present study, the SWP values between the stressed and well-watered vines [81] can be easily identified. The relationship between the CWSI and SWP improved when the shadow pixels were removed from the vine canopy using the suggested automatic algorithms. For the vineyard, the fractional cover was 19%, while the percentage of canopy shaded pixels was 43%. This indicates that only 8.2% of final vegetation pixels were used to develop the relationship between CWSI and SWP.

Discussion
The selection of B1 (490 nm) as the better multispectral band for classifying shadow canopy pixels was consistent with the previous study by Ünsalan et al. [76], who used the k-means and blue information derived from the RGB spectrum to segment information avoiding shadow pixels to extract street networks and detect houses. This band selection was also proposed by Sirmacek et al. [77], who used the blue wavelength spectrum to detect shadows for building detection, suggesting that this region was dominantly better even compared with green and red for shadow pixels identification [78,79]. The selection of 490 nm image was also preferred when compared with upper wavelengths, in which blue spectrum showed better results for shadow detection increasing the performance for near infrared and shortwave infrared [80]. This validates the previous assumption that internal canopy shadow cannot be identified by TIR imagery. Considering this, the importance of coregister thermal and visible images for detecting shadow pixels was also highlighted by Leinonen et al. [42] who using ground cameras with a non-automatic methodology concluded that one of the principal steps is to correct overlapping VIS and TIR images to assess vine water status.
In the present study, the SWP values between the stressed and well-watered vines [81] can be easily identified. The relationship between the CWSI and SWP improved when the shadow pixels were removed from the vine canopy using the suggested automatic algorithms. For the vineyard, the fractional cover was 19%, while the percentage of canopy shaded pixels was 43%. This indicates that only 8.2% of final vegetation pixels were used to develop the relationship between CWSI and SWP. Although the relationship between CWSI and SWP improved, the impact of the shadow was significant in those vines with more water stress [51]. In contrast, because no reduction of the transpiration rate occurred in well-watered vines [82], the difference between leaf temperature and air temperature was not representative [83]. These results are consistent with those of Van Zyl [49], who suggested that the impact of shadow for SWP relationships in stressed vines was considerably higher compared with sunlit leaves. In addition, Pou et al. [50] suggested that shaded canopy information negatively affects the relationship of the vine water status and CSWI because the leaf temperature decreases. Furthermore, Jones et al. [82] suggested that a greater sensitivity with respect to leaf temperature with water status measurements might be better when sunlit canopy information is considered. The effect of shadow deletion on the relationship between CWSI and SWP for stressed and well-watered vines is shown in Figure 10. Considering the results of Figure 9, in stressed vines, the shadow deletion process significantly improved the CWSI-SWP relationship with values of R 2 , increasing from 0.05 to 0.35. However, no differences were observed for well-watered vines. to leaf temperature with water status measurements might be better when sunlit canopy information is considered. The effect of shadow deletion on the relationship between CWSI and SWP for stressed and well-watered vines is shown in Figure 10. Considering the results of Figure 9, in stressed vines, the shadow deletion process significantly improved the CWSI-SWP relationship with values of R 2 , increasing from 0.05 to 0.35. However, no differences were observed for well-watered vines.

Conclusions
Using a modified SIFT computer vision algorithm and Kmeans++ clustering, we performed automatic coregister UAV-TIR and UAV-VIS imagery to detect canopy shadow pixels information in thermal images. The deletion of the canopy shadow information in TIR images positively affects the relationship between the CWSI and SWP, showing an increment in R 2 from 0.64 to 0.77. In addition, the relationship showed a decrease in RMSE from 0.2 to 0.1 MPa and in SE from 0.24 to 0.16 MPa. As future work, our methodology should be applied for validation in different cultivars, seasons, and field conditions. In addition, the impact of automatic removal of shadow canopy pixels should be assessed for evapotranspiration modeling using UAV-TIR images of vineyards.

Conclusions
Using a modified SIFT computer vision algorithm and Kmeans++ clustering, we performed automatic coregister UAV-TIR and UAV-VIS imagery to detect canopy shadow pixels information in thermal images. The deletion of the canopy shadow information in TIR images positively affects the relationship between the CWSI and SWP, showing an increment in R 2 from 0.64 to 0.77. In addition, the relationship showed a decrease in RMSE from 0.2 to 0.1 MPa and in SE from 0.24 to 0.16 MPa. As future work, our methodology should be applied for validation in different cultivars, seasons, and field conditions. In addition, the impact of automatic removal of shadow canopy pixels should be assessed for evapotranspiration modeling using UAV-TIR images of vineyards.