Random Forest Algorithm Improves Detection of Physiological Activity Embedded within Reﬂectance Spectra Using Stomatal Conductance as a Test Case

: Plants transpire water through their tissues in order to move nutrients and water to the cells. Transpiration includes various mechanisms, primarily stomata movement, which controls the rate of CO 2 and water vapor exchange between the tissues and the atmosphere. Assessment of stomatal conductance is available for gas exchange techniques at leaf level, yet these techniques are not scalable to the whole plant let alone a large vegetation area. Hyperspectral reﬂectance spectroscopy, which acquires hundreds of bands in a single scan, may capture a glimpse of the crop’s physiological activity and therefore meet the scalability challenge. In this study, classic chemometric analyses are used alongside advanced statistical learning algorithms in order to identify stomatal conductance cues in hyperspectral measurements of cotton plants experiencing a gradient of irrigation. Random forest of regression trees identiﬁed 23 wavelengths related to both structural properties of the plant as well as water content. Partial least squares regression succeeded in relating these wavelengths to stomatal conductance, but only partially (R 2 < 0.2). An artiﬁcial neural network algorithm reported an R 2 = 0.54 with an 89% error-free performance on the same data subset. This study discusses implementation of machine learning methodologies as a benchmark for deeper analysis of spectral information, such as required when searching for plant physiology-related attenuations embedded within reﬂectance spectra. There were 18 biological repeats per irrigation treatment times 4 technical repeats in each year for a total of 432 pots. Three cotton varieties— G. hirsutum , G. barbadense , and G. akalpi —were used in the study (they presented no di ﬀ erence in stomatal conductance within each irrigation treatment). Two irrigation experiments were performed: (a) irrigation volume gradient—four daily irrigation rates were applied at 4, 3, 2, and 1 L / pot, corresponding to midday water potentials of − 1.4, − 1.8, − 2.0, and − 2.5 MPa, respectively; water potential (Scholander bomb, MRC, Israel) measured at noon reﬂected the plant’s need for water; each repeat (irrigation treatment * variety) consisted of 18 pots, and each irrigation treatment * variety was replicated 6 times in random blocks; (b) irrigation shut-o ﬀ experiment—all the quads received optimal irrigation volume, which was determined as 3 L per day ( − 1.8 MPa at noon); at t = 0, the irrigation was stopped for 24 h and resumed in an irrigation gradient with the same four volumes treatments as described above per day for a week. The 2019 experiment included a repeat on the volumetric gradient design part performed in 2018 on G. akalpi variety only. Measurements in all the experiments were conducted in the same phenological period. acquired spectra were stitched together. Overall, there were 1222 wavelengths with 1 nm and 6 nm full width at height maximum (FWHM) spectral resolution for STS and FLAME spectrometers, respectively. The unit was of a with a that situated the spectrometer a 2.5 m above the The ground was situated above 94% NH, The spectrometers’ stitched spectra were spectrally calibrated according to the atmospheric water absorption band at 970 nm. Small deviations on the spectral range between STS and FLAME spectrometers corrected to Rascher


Introduction
Transpiration is a physiological process within the plant that allows for water movement through tissues, starting at the root and exiting through specialized pores within the leaves, termed stomata [1]. Complex physiological and regulatory processes drive transpiration, the primary of which is stomatal conductance [2]. The stomata control the exchange of water vapor and CO 2 between the leaf and the atmosphere by promptly changing the stomata aperture size according to the plant needs. Stomatal conductance reduces transpiration when plants experience water stress, therefore, it controls the water status of the plant and its tolerance to drought. Yet, this step reduces photosynthetic rates as well. Therefore, detecting this process on a global scale is a fundamental requirement to developing models for ecological resilience of crops during climate change, global ecology, and, in practice, irrigation management in agriculture [3]. Measuring gas exchange at leaf level through current techniques is not scalable to large vegetation areas due to the technical difficulties of sampling (HL-2000-LL, OceanInsight, FL, USA) according to the manufacturer's instructions. STS and FLAME spectrometers obtained an overlapping region in the spectral range of 936-1120 nm, and on this basis, each two acquired spectra were stitched together. Overall, there were 1222 wavelengths with 1 nm and 6 nm full width at height maximum (FWHM) spectral resolution for STS and FLAME spectrometers, respectively. The air unit was mounted on top of a trolley with a boom that situated the spectrometer pupil at a consistent 2.5 m above the sampled canopy. The ground unit was situated above a 94% reflectance white plate (Permaflect © ; LabSphere, NH, USA). The spectrometers' stitched spectra were spectrally calibrated according to the atmospheric water absorption band at 970 nm. Small deviations on the spectral range between STS and FLAME spectrometers were corrected according to Rascher et al. [16], and the radiometric deviations between the two units were corrected as well (when both were pointing at the reference plate). Thus, the water absorption band at 970 nm obtained the same magnitude and physical properties in both STS and FLAME spectrometers. Reflectance signatures were calculated and corrected for meteorological conditions according to Gordon and Wang's study [17].

Stomatal Conductance Datasets
Various irrigation regimes on cotton plants were the basis for this study, aiming to maintain a wide range of the plants' water status and stomatal conductance values. The irrigation experiments were conducted in 2018 and 2019 ( Figure 1). The cotton plants were grown in 3.9 L pots filled with mixed ground (80:20 Peat soil:Clay, Kekkila BVB, Sweden). Four seeds were sown in each pot in a square geometry and, after emergence, the two weakest plants within each pot were removed. Every four pots within a quad were closely situated so that, early in the season, the canopies of adjacent pots were reducing the ground view. The pots were drip-irrigated with a fertigation solution, thus the fertilizer concentration was relative to the irrigation volume (Sheffer + 3 micro, ICL, Tel Aviv-Yafo, Israel). There were 18 biological repeats per irrigation treatment times 4 technical repeats in each year for a total of 432 pots. Three cotton varieties-G. hirsutum, G. barbadense, and G. akalpi-were used in the study (they presented no difference in stomatal conductance within each irrigation treatment). Two irrigation experiments were performed: (a) irrigation volume gradient-four daily irrigation rates were applied at 4, 3, 2, and 1 L/pot, corresponding to midday water potentials of −1.4, −1.8, −2.0, and −2.5 MPa, respectively; water potential (Scholander bomb, MRC, Israel) measured at noon reflected the plant's need for water; each repeat (irrigation treatment * variety) consisted of 18 pots, and each irrigation treatment * variety was replicated 6 times in random blocks; (b) irrigation shut-off experiment-all the quads received optimal irrigation volume, which was determined as 3 L per day (−1.8 MPa at noon); at t = 0, the irrigation was stopped for 24 h and resumed in an irrigation gradient with the same four volumes treatments as described above per day for a week. The 2019 experiment included a repeat on the volumetric gradient design part performed in 2018 on G. akalpi variety only. Measurements in all the experiments were conducted in the same phenological period.
Remote Sens. 2020, 12, x FOR PEER REVIEW 3 of 19 acquired spectra were stitched together. Overall, there were 1222 wavelengths with 1 nm and 6 nm full width at height maximum (FWHM) spectral resolution for STS and FLAME spectrometers, respectively. The air unit was mounted on top of a trolley with a boom that situated the spectrometer pupil at a consistent 2.5 m above the sampled canopy. The ground unit was situated above a 94% reflectance white plate (Permaflect © ; LabSphere, NH, USA). The spectrometers' stitched spectra were spectrally calibrated according to the atmospheric water absorption band at 970 nm. Small deviations on the spectral range between STS and FLAME spectrometers were corrected according to Rascher et al. [16], and the radiometric deviations between the two units were corrected as well (when both were pointing at the reference plate). Thus, the water absorption band at 970 nm obtained the same magnitude and physical properties in both STS and FLAME spectrometers. Reflectance signatures were calculated and corrected for meteorological conditions according to Gordon and Wang's study [17].

Stomatal Conductance Datasets
Various irrigation regimes on cotton plants were the basis for this study, aiming to maintain a wide range of the plants' water status and stomatal conductance values. The irrigation experiments were conducted in 2018 and 2019 ( Figure 1). The cotton plants were grown in 3.9 L pots filled with mixed ground (80:20 Peat soil:Clay, Kekkila BVB, Sweden). Four seeds were sown in each pot in a square geometry and, after emergence, the two weakest plants within each pot were removed. Every four pots within a quad were closely situated so that, early in the season, the canopies of adjacent pots were reducing the ground view. The pots were drip-irrigated with a fertigation solution, thus the fertilizer concentration was relative to the irrigation volume (Sheffer + 3 micro, ICL, Israel). There were 18 biological repeats per irrigation treatment times 4 technical repeats in each year for a total of 432 pots. Three cotton varieties--G. hirsutum, G. barbadense, and G. akalpi-were used in the study (they presented no difference in stomatal conductance within each irrigation treatment). Two irrigation experiments were performed: (a) irrigation volume gradient-four daily irrigation rates were applied at 4, 3, 2, and 1 L/pot, corresponding to midday water potentials of −1.4, −1.8, −2.0, and −2.5 MPa, respectively; water potential (Scholander bomb, MRC, Israel) measured at noon reflected the plant's need for water; each repeat (irrigation treatment * variety) consisted of 18 pots, and each irrigation treatment * variety was replicated 6 times in random blocks; (b) irrigation shut-off experiment-all the quads received optimal irrigation volume, which was determined as 3 L per day (−1.8 MPa at noon); at t = 0, the irrigation was stopped for 24 h and resumed in an irrigation gradient with the same four volumes treatments as described above per day for a week. The 2019 experiment included a repeat on the volumetric gradient design part performed in 2018 on G. akalpi variety only. Measurements in all the experiments were conducted in the same phenological period.  Hyperspectral acquisition and abaxial stomatal conductance measurements (AP4, Delta-t, UK in 2018 and Licor 6800 photosynthesis system, Li-Cor Biosciences, IL, USA in 2019) were performed for two months covering vegetative growth, transfer to, and start of the reproductive growth stage. An AP4 porometer was calibrated with a calibration plate before the measurement sets according to the manufacturer's instructions. In the case of Li-Cor 6800, an open chamber was used in order to receive sunlight at a natural intensity via a 3 × 3 cm window. Environmental conditions were as follows: the temperature and the relative humidity were set according to the meteorological station parameters (Table 1). CO 2 level was always set to 400 ppm. Measurements were recorded when assimilation and conductance reached a plateau, approximately 30 s to 1 min after closing the chamber on the leaf. In each date, measurements were taken twice (once before noon and then again at noon), where ten measurements per pot quad were taken each season (Table 1 describes meteorological conditions during the measurements received from a local meteorological station). In both years, the meteorological station was about~85 m southwest of the experimental plot.

Pre-Processing of Data for Important Features Selection
We followed Rinnan et al. [18] protocols to prepare the spectra for analysis ( Figure 2). First, we used a box-car averaging technique in order to create an even, nominal pace between the two sides of the stitched spectra. That is, we eventually received 231 wavelengths of the 1222 total wavelengths of the raw data. Then, outlier spectra were identified by Cochran's test [19] and removed from the data set, and a standard normal variate (SNV) test followed in order to correct for multiple scatter [20]. Next, spectra were corrected for their additive dispersion effect with a baseline correction and were normalized to the maximum peak within each spectrum. The procedure resulted in a negative local minimum at the water absorption band (between 1380-1450 nm). Therefore, these wavelengths were omitted from further usage and analysis. Finally, the spectra were mean-centered and standardized before performing further analysis on it. Stomatal conductance data was also curated for outliers per irrigation treatments, where if sample values were far from the average (as in, more than one standard deviation distance), the whole sample was omitted from the analysis.

Wavelength Selection
A Normalized Difference Index (NDI) combinations technique, also termed contour-contour map [21], was used as the first step to find a simple relationship between the spectra and the stomatal conductance. A coefficient of determination was the examined value. NDI formulations states:

Pre-Processing of Data for Important Features Selection
We followed Rinnan et al. [18] protocols to prepare the spectra for analysis ( Figure 2). First, we used a box-car averaging technique in order to create an even, nominal pace between the two sides of the stitched spectra. That is, we eventually received 231 wavelengths of the 1222 total wavelengths of the raw data. Then, outlier spectra were identified by Cochran's test [19] and removed from the data set, and a standard normal variate (SNV) test followed in order to correct for multiple scatter VPD = Vapor Pressure Deficit.

Pre-Processing of Data for Important Features Selection
We followed Rinnan et al. [18] protocols to prepare the spectra for analysis ( Figure 2). First, we used a box-car averaging technique in order to create an even, nominal pace between the two sides of the stitched spectra. That is, we eventually received 231 wavelengths of the 1222 total wavelengths of the raw data. Then, outlier spectra were identified by Cochran's test [19] and removed from the data set, and a standard normal variate (SNV) test followed in order to correct for multiple scatter [20]. Next, spectra were corrected for their additive dispersion effect with a baseline correction and were normalized to the maximum peak within each spectrum. The procedure resulted in a negative local minimum at the water absorption band (between 1380-1450 nm). Therefore, these wavelengths were omitted from further usage and analysis. Finally, the spectra were mean-centered and standardized before performing further analysis on it. Stomatal conductance data was also curated for outliers per irrigation treatments, where if sample values were far from the average (as in, more than one standard deviation distance), the whole sample was omitted from the analysis.

Pre-Processing of Data for Important Features Selection
We followed Rinnan et al. [18] protocols to prepare the spectra for analysis ( Figure 2). First, we used a box-car averaging technique in order to create an even, nominal pace between the two sides of the stitched spectra. That is, we eventually received 231 wavelengths of the 1222 total wavelengths of the raw data. Then, outlier spectra were identified by Cochran's test [19] and removed from the data set, and a standard normal variate (SNV) test followed in order to correct for multiple scatter [20]. Next, spectra were corrected for their additive dispersion effect with a baseline correction and were normalized to the maximum peak within each spectrum. The procedure resulted in a negative local minimum at the water absorption band (between 1380-1450 nm). Therefore, these wavelengths were omitted from further usage and analysis. Finally, the spectra were mean-centered and standardized before performing further analysis on it. Stomatal conductance data was also curated for outliers per irrigation treatments, where if sample values were far from the average (as in, more than one standard deviation distance), the whole sample was omitted from the analysis.

Wavelength Selection
A Normalized Difference Index (NDI) combinations technique, also termed contour-contour map [21], was used as the first step to find a simple relationship between the spectra and the stomatal conductance. A coefficient of determination was the examined value. NDI formulations states: where ρi,j are the reflectance values at a designated wavelength i,j in units of the measurement (in our study nanometer (nm)), and a, b, c, d are rational number coefficients selected by the computer during minimization of sum of squares. The coefficient values (Equation (1)) were either set to "1"

Wavelength Selection
A Normalized Difference Index (NDI) combinations technique, also termed contour-contour map [21], was used as the first step to find a simple relationship between the spectra and the stomatal conductance. A coefficient of determination was the examined value. NDI formulations states: where ρ i,j are the reflectance values at a designated wavelength i,j in units of the measurement (in our study nanometer (nm)), and a, b, c, d are rational number coefficients selected by the computer during minimization of sum of squares. The coefficient values (Equation (1)) were either set to "1" or were determined using the optimum minimize package in Python Software (Python, DE, USA) [22], operated with Constrained Optimization BY Linear Approximation process (COBYLA) [23]. A random forest ensemble method followed [24]; it allows the division of the dataset selected into nodes within a tree. Each node contains a subset of the original dataset. This subset comprises a subset of samples and a subset of the wavelengths. Each node is then averaged, and the number obtained is compared via root mean squared error (RMSE) to the measure stomatal conductance values. The target of the ensemble is to lower the RMSE as much as possible where important node dividers are flagged and kept. In this study, the divider was an arbitrary wavelength. Therefore, the selected architecture of the trees was the one with the lowest achieved RMSE and a set of dividers (wavelengths) that enabled this to occur. During the divisions into nodes, essentially deep in the tree, there could be left only one wavelength. In order to avoid it, pruning of the forest, i.e., stopping the tree from being developed [25], was performed in three levels: (a) numbers of trees in the forest (50, 100, 250, 500); (b) constant percentage of attributes within each branch split (10%, 20%, 30%, 40% of the attributes); (c) maximum depth of the regression tree (selected to be to two-thirds of an average maximum depth reached in preliminary iterations on the dataset).

Validation of Features Selected
The dataset was divided into 75% calibration sub-set and 25% testing sub-set. Initially, a standard multiple linear regression model was searched [26]. However, due to a violation of one of the model's prerequisites, namely the inability to obtain a partial linear regression between each of the predictors and the predictand, a partial least square regression process was used instead [27]. This algorithm projects the original dataset into a latent structure where each of the new variables is independent of each other and therefore minimizes the co-linearity problem. It lowers the dimensionality of the problem, and the user has a greater chance to create a statistical model of spectra calculating stomatal conductance. In this study, the dimension of the problem was reduced from 23 selected features into a set of 4 new latent predictors. The model was cross-validated using a venetian blinds sub-sampling cross validation. With unimproved coefficient of determination, an artificial neural network (ANN) was used afterwards. This model accounts for the nonlinear relationship between the spectra and the physiological process [28]. The ANN architecture [29] included one hidden layer and a standard back-propagation process with Adam's optimizer containing a loss function [30]. Performance of the ANN was validated as well, with a suite of statistical tests as suggested by Sousa et al. [28].

Statistical Packages
Dataset pre-processing, partial least squares regression, and k-means clustering were performed in UnscrumblerX software (Camo Analytics, AS, Norway). Data preparation for multiple linear regression was performed in SPSS (IBM SPSS, Armonk, NY, USA). Random forest of regression trees and artificial neural network were coded in Python programming language [22].

Results
The two most irrigated treatments had similar stomatal conductance ( Figure 3A) and decreased with further reduction in water potential, as expected for plants experiencing water stress. For each stomatal conductance measurement, a spectral acquisition followed ( Figure 3B-D). A portion of the short-wave infra-red (SWIR) region was removed from the analysis (see the region between 1380-1450 nm in Figure 2B) due to its negative appearance during pre-processing of the dataset. This happened because there was a strong absorption of the water band in this spectral region [31].

of 18
The two most irrigated treatments had similar stomatal conductance ( Figure 3A) and decreased with further reduction in water potential, as expected for plants experiencing water stress. For each stomatal conductance measurement, a spectral acquisition followed ( Figure 3B-D). A portion of the short-wave infra-red (SWIR) region was removed from the analysis (see the region between 1380-1450 nm in Figure 2B) due to its negative appearance during pre-processing of the dataset. This happened because there was a strong absorption of the water band in this spectral region [31]. The only visible difference could be seen between the blue and the yellow colored curves: -2.5 and −1.4 MPa water potential, respectively. These were the two extreme conditions during the study: either over irrigating or wilting until it reached a very low water potential (−2.5 MPa) and 33% of the maximum stomatal conductance. When magnifying the spectral regions, it could also be seen that the trends in the reflectance values were very similar to the relations between the actual values of the stomatal conductance ( Figure 3C and 3D compared to 3A). The magnitude of the spectrum in the VISible-Near Infra Red (VIS-NIR) range increased with an increase in stomatal conductance. The spectra in the Short-Wave Infra-Red (SWIR) spectral range presented an opposite behavior to that in the VIS range (compare Figure 3C to Figure 3D). Here, the SWIR region was affected by water absorption bands, and therefore it was opposite because the plant retained more water within its tissues. As such, the water absorbed more light, and the reflectance in this region decreased. Finally, although there was a gradual decrease between the four irrigation treatments, it was visible that the extreme water stress (blue color) affected the spectrum more severely in the VIS region. Its magnitude was much smaller when compared to the other water treatments.
In the following sections, we discuss the wavelength selection process and the validation of that selection as a regression model between spectral and stomatal conductance measurements. The only visible difference could be seen between the blue and the yellow colored curves: -2.5 and −1.4 MPa water potential, respectively. These were the two extreme conditions during the study: either over irrigating or wilting until it reached a very low water potential (−2.5 MPa) and 33% of the maximum stomatal conductance. When magnifying the spectral regions, it could also be seen that the trends in the reflectance values were very similar to the relations between the actual values of the stomatal conductance ( Figure 3C,D compared to Figure 3A). The magnitude of the spectrum in the VISible-Near Infra Red (VIS-NIR) range increased with an increase in stomatal conductance. The spectra in the Short-Wave Infra-Red (SWIR) spectral range presented an opposite behavior to that in the VIS range (compare Figure 3C to Figure 3D). Here, the SWIR region was affected by water absorption bands, and therefore it was opposite because the plant retained more water within its tissues. As such, the water absorbed more light, and the reflectance in this region decreased. Finally, although there was a gradual decrease between the four irrigation treatments, it was visible that the extreme water stress (blue color) affected the spectrum more severely in the VIS region. Its magnitude was much smaller when compared to the other water treatments.
In the following sections, we discuss the wavelength selection process and the validation of that selection as a regression model between spectral and stomatal conductance measurements.

Wavelength Selection
Normalized Difference Index technique relates the chemical quotients of plants, such as nitrogen [32], water content [7], and photosynthetic pigments contents [33], to linear combinations of every two wavelengths within the acquired spectra [34]. Using this technique on 231 wavelengths in our dataset, we found four hot-spot regions that were more correlated to the stomatal conductance values than the rest of the wavelength combinations ( Figure 4A): 693-703 nm, 780-890 nm, 1007-1120 nm, 1500-1560 nm.
Remote Sens. 2020, 12, 2213 8 of 18 [32], water content [7], and photosynthetic pigments contents [33], to linear combinations of every two wavelengths within the acquired spectra [34]. Using this technique on 231 wavelengths in our dataset, we found four hot-spot regions that were more correlated to the stomatal conductance values than the rest of the wavelength combinations ( Figure 4A): 693 nm-703 nm, 780 nm-890 nm, 1007-1120 nm, 1500 nm-1560 nm. However, the coefficient of determination was quite low (R 2 <0.2). NDI index was calculated with the combination that reached the highest R 2 -1094 nm and 1096 nm ( Figure 4B). However, it was very poor due to natural clustering of the data within the correlation plot.
We therefore hypothesized that the detection of hot spot regions within the contour map could be improved if a generalization of the method were to take place. This means that coefficients a, b, c, and d in Equation (1) would receive an arbitrary rational number. In such a case, the NDI method shown in ( Figure 4A) is a private case of the generalized equation, where all the coefficients equal one. Such a step increases the capability of the statistical process to minimize the RMSE and reach a higher coefficient of determination. This model was run on each of the 231 wavelength combinations. Using this approach, we found many more bi-wavelength combinations, which almost tripled the coefficient of determination values found in the private case of NDI ( Figure 4C,D). The contour map in this setup revealed that the range between 634 nm-843 nm related to stomatal conductance both within itself (parallel bold red line in Figure 4C) and within two main regions in the spectrum: 843 nm-1047 nm and 1455-1661 nm. A maximum correlation was achieved for the two-wavelength combination of 734 nm and 830 nm with R 2 = 0.19, as in the private case. Still, it was visible that the correlation to the actual regression line was poor ( Figure 4D). The reason here is similar to the private case, where the generalization of the NDI technique managed to lower the number of natural clusters to two within the correlation plot. However, the coefficient of determination was quite low (R 2 < 0.2). NDI index was calculated with the combination that reached the highest R 2 -1094 nm and 1096 nm ( Figure 4B). However, it was very poor due to natural clustering of the data within the correlation plot.
We therefore hypothesized that the detection of hot spot regions within the contour map could be improved if a generalization of the method were to take place. This means that coefficients a, b, c, and d in Equation (1) would receive an arbitrary rational number. In such a case, the NDI method shown in ( Figure 4A) is a private case of the generalized equation, where all the coefficients equal one. Such a step increases the capability of the statistical process to minimize the RMSE and reach a higher coefficient of determination. This model was run on each of the 231 wavelength combinations. Using this approach, we found many more bi-wavelength combinations, which almost tripled the coefficient of determination values found in the private case of NDI ( Figure 4C,D). The contour map in this setup revealed that the range between 634-843 nm related to stomatal conductance both within itself (parallel bold red line in Figure 4C) and within two main regions in the spectrum: 843-1047 nm and 1455-1661 nm. A maximum correlation was achieved for the two-wavelength combination of 734 nm and 830 nm with R 2 = 0.19, as in the private case. Still, it was visible that the correlation to the actual regression line was poor ( Figure 4D). The reason here is similar to the private case, where the generalization of the NDI technique managed to lower the number of natural clusters to two within the correlation plot.
We searched for a better and more advanced technique to correlate between spectra and stomatal conductance. Specifically, we searched for a technique that can isolate a lower number of features from the spectral range (one feature in this sense is one wavelength) and consider the non-linear relationship between the physiological phenomenon and the spectral information. This can be expected from researching biological activity such as stomatal conductance, which encompasses both chemical and structural attenuations of the sample measured. We selected a supervised machine learning algorithm random forest (RF) of regression trees [24] after the work of Abdel-Rahman and colleagues implementing RF for analysis of spectral measurements in general [35] (Figure 5). The algorithm performs an implicit feature selection (in this case, the features are the spectral wavelengths) and decides during the iterations which features are the best data selectors per a final dataset. In other words, the RF algorithm selects a subset of wavelengths and a subset of samples each time and compares the calculated stomatal conductance coming from the spectra to that of the measured stomatal conductance of that specific data subset (one regression node within one tree within the forest). Thus, it seeks to minimize an RMSE function between the two measures. The wavelengths that were able to obtain the lowest RMSE in each iteration over all the nodes and the trees in the forest are then flagged and regarded as "important" by the algorithm. This algorithm is highly resilient to over-fitting problems that may arise due to noise within the dataset by generating of a large number of data subsets (trees) [36].
features from the spectral range (one feature in this sense is one wavelength) and consider the nonlinear relationship between the physiological phenomenon and the spectral information. This can be expected from researching biological activity such as stomatal conductance, which encompasses both chemical and structural attenuations of the sample measured. We selected a supervised machine learning algorithm random forest (RF) of regression trees [24] after the work of Abdel-Rahman and colleagues implementing RF for analysis of spectral measurements in general [35] (Figure 5). The algorithm performs an implicit feature selection (in this case, the features are the spectral wavelengths) and decides during the iterations which features are the best data selectors per a final dataset. In other words, the RF algorithm selects a subset of wavelengths and a subset of samples each time and compares the calculated stomatal conductance coming from the spectra to that of the measured stomatal conductance of that specific data subset (one regression node within one tree within the forest). Thus, it seeks to minimize an RMSE function between the two measures. The wavelengths that were able to obtain the lowest RMSE in each iteration over all the nodes and the trees in the forest are then flagged and regarded as "important" by the algorithm. This algorithm is highly resilient to over-fitting problems that may arise due to noise within the dataset by generating of a large number of data subsets (trees) [36]. , and the maximum depth to which a tree may divide the dataset ( Figure 5 panels a-k), as suggested by [25]. Each point within the panels represented a different set of flagged wavelengths. The general behavior of the algorithm shows that, with an increase in branch number, the RMSE decreased in a logarithmic trend. Only after passing one-third of the total depth of the trees, the trend became noisy between the various numbers of trees within the forest (compare Figure 5E to Figure 5A). When the depth of the tree passed half of its maximum depth (especially in forests with the minimum number of trees, see the blue colored curves in Figure 5), the RMSE trend revealed a local minima and complete dispersion of the trend. For example, in Figure 5J, each trend obtained a different form with a local minima for each size of the forests. The lowest RMSE was not necessarily reached for the maximum number of trees, depth, and samples. Instead, the best architecture of the random forest was found to be at about two-thirds of the maximum depth, with 20 percent selection of attributes within each branch split and 250 trees in the forest ( Figure 5G, bold black arrow).

Validation of Feature Selection
The random forest algorithm can predict parameters by its non-linear regression algorithm [37], yet, its prediction capability is limited to the calibrated dataset [38]. Therefore, in search of a viable equation or model that can relate actual stomatal conductance to the selected wavelengths, we first tried to build a multi-linear regression model (Supplementary Figure S1). The model could not be assembled due to a violation of the predominant assumption that each of the predictors obtains a partial linear relationship with the dependent variable. This was probably because some of the features selected by the RF mechanism were correlated. In order to neutralize that co-linearity, we projected the data into a latent structure and assembled a partial least squares regression [27] (Figure 6).
behavior of the algorithm shows that, with an increase in branch number, the RMSE decreased in a logarithmic trend. Only after passing one-third of the total depth of the trees, the trend became noisy between the various numbers of trees within the forest (compare Figure 5E to Figure 5A). When the depth of the tree passed half of its maximum depth (especially in forests with the minimum number of trees, see the blue colored curves in Figure 5), the RMSE trend revealed a local minima and complete dispersion of the trend. For example, in Figure 5J, each trend obtained a different form with a local minima for each size of the forests. The lowest RMSE was not necessarily reached for the maximum number of trees, depth, and samples. Instead, the best architecture of the random forest was found to be at about two-thirds of the maximum depth, with 20 percent selection of attributes within each branch split and 250 trees in the forest ( Figure 5G, bold black arrow).

Validation of Feature Selection
The random forest algorithm can predict parameters by its non-linear regression algorithm [37], yet, its prediction capability is limited to the calibrated dataset [38]. Therefore, in search of a viable equation or model that can relate actual stomatal conductance to the selected wavelengths, we first tried to build a multi-linear regression model (Supplementary Figure S1). The model could not be assembled due to a violation of the predominant assumption that each of the predictors obtains a partial linear relationship with the dependent variable. This was probably because some of the features selected by the RF mechanism were correlated. In order to neutralize that co-linearity, we projected the data into a latent structure and assembled a partial least squares regression [27] ( Figure  6).  Partial least squares regression is a regularized regression method able to lower the dimensionality of the spectra through projection into a latent structure, where the new set of variables, each based on a combination of the original wavelengths, are now independent of each other. These latent variables preserve most of the variance of the original dataset. The data were divided into a 75:25 training/testing subset. During calibration of the model, the model presented a similar underlying representation of the stomatal conductance between the two years of experimentation-note the general shape of the red dots compared to the blue dots ( Figure 6A). We tried to identify meaningful clusters within the scores plot by using the discrete dates during each year and water potentials (Supplementary Figure S2). However, we did not have success in associating either classification with stomatal conductance values. We did manage to divide the calibration set into four clusters using a squared Euclidian distance range within a k-means clusters technique, however, this was not correlated with any of our predictor classes. The optimum model was selected across four components, where both the calibration and the cross-validation tests reached about 80% of the explained variance at the fourth principal component ( Figure 6C). However, the projected model only succeeded to predict at R 2 = 0.23 of the behavior of the stomatal conductance in the test subset. Thus, it was not much better than the simpler method of NDI, but it did manage to cancel out the natural clustering seen in the NDI technique correlation plot. Specifically, the model had difficulties in predicting the stomatal conductance values at the range of 400-600 mmol H 2 O m −2 s −1 . These were the treatments with a high stomatal conductance, where the cotton plant was adequately irrigated, and the explanation for this difficulty can be seen clearly upon inspecting Figure 3A. While the -1.4 MPa treatment should have been at least at the same level of the -1.8 MPa treatment, it declined, even if not statistically significant, and this was probably the reason for the failure of the predictive model. The PLS-R cannot decide which of the original features are the most important to the model, and thus we implemented an ANN instead (Figure 7). A standard ANN architecture consists of input nodes that receive the important wavelengths for prediction of stomatal conductance. They are then transferred towards a hidden layer where, along the way, they receive new values via a non-linear transfer function (usually a sigmoid function). Then, they are transferred again into the output node, this time with a linear function. Finally, the ANN calculates a value of stomatal conductance and compares it to the original value. The weight matrix that was created along the way is changed repeatedly until the bias term is minimized and the weights function is optimized [30]. The desired result should be a linear correlation between the spectral information and the stomatal conductance values.
other. These latent variables preserve most of the variance of the original dataset. The data were divided into a 75:25 training/testing subset. During calibration of the model, the model presented a similar underlying representation of the stomatal conductance between the two years of experimentation-note the general shape of the red dots compared to the blue dots ( Figure 6A). We tried to identify meaningful clusters within the scores plot by using the discrete dates during each year and water potentials (Supplementary Figure S2). However, we did not have success in associating either classification with stomatal conductance values. We did manage to divide the calibration set into four clusters using a squared Euclidian distance range within a k-means clusters technique, however, this was not correlated with any of our predictor classes. The optimum model was selected across four components, where both the calibration and the cross-validation tests reached about 80% of the explained variance at the fourth principal component ( Figure 6C). However, the projected model only succeeded to predict at R 2 = 0.23 of the behavior of the stomatal conductance in the test subset. Thus, it was not much better than the simpler method of NDI, but it did manage to cancel out the natural clustering seen in the NDI technique correlation plot. Specifically, the model had difficulties in predicting the stomatal conductance values at the range of 400-600 mmol H2O m −2 s -1 . These were the treatments with a high stomatal conductance, where the cotton plant was adequately irrigated, and the explanation for this difficulty can be seen clearly upon inspecting Figure 3A. While the -1.4 MPa treatment should have been at least at the same level of the -1.8 MPa treatment, it declined, even if not statistically significant, and this was probably the reason for the failure of the predictive model. The PLS-R cannot decide which of the original features are the most important to the model, and thus we implemented an ANN instead (Figure 7). A standard ANN architecture consists of input nodes that receive the important wavelengths for prediction of stomatal conductance. They are then transferred towards a hidden layer where, along the way, they receive new values via a non-linear transfer function (usually a sigmoid function). Then, they are transferred again into the output node, this time with a linear function. Finally, the ANN calculates a value of stomatal conductance and compares it to the original value. The weight matrix that was created along the way is changed repeatedly until the bias term is minimized and the weights function is optimized [30]. The desired result should be a linear correlation between the spectral information and the stomatal conductance values.  We constructed a back-propagated one hidden layer standard architecture, which was proven in the past to correctly predict non-linear mathematical relations (Figure 7) [29]. The ANN algorithm succeeded in creating a linear relationship between 20 features out of the 23 features originally found by the RF algorithm. This model can reach stomatal conductance out of spectral information on the test-subset with R 2 = 0.54 accuracy. We also checked the performance of the procedure with various statistical tests (Table 2) [28].

Pre-Processing of Data for Important Features Selection
We followed Rinnan et al. [18] protocols to prepare the spectra for analysis ( Figure 2). First, we used a box-car averaging technique in order to create an even, nominal pace between the two sides of the stitched spectra. That is, we eventually received 231 wavelengths of the 1222 total wavelengths of the raw data. Then, outlier spectra were identified by Cochran's test [19] and removed from the data set, and a standard normal variate (SNV) test followed in order to correct for multiple scatter [20]. Next, spectra were corrected for their additive dispersion effect with a baseline correction and were normalized to the maximum peak within each spectrum. The procedure resulted in a negative local minimum at the water absorption band (between 1380-1450 nm). Therefore, these wavelengths were omitted from further usage and analysis. Finally, the spectra were mean-centered and standardized before performing further analysis on it. Stomatal conductance data was also curated for outliers per irrigation treatments, where if sample values were far from the average (as in, more than one standard deviation distance), the whole sample was omitted from the analysis.

Wavelength Selection
A Normalized Difference Index (NDI) combinations technique, also termed contour-contour map [21], was used as the first step to find a simple relationship between the spectra and the stomatal conductance. A coefficient of determination was the examined value. NDI formulations states: (1) where ρi,j are the reflectance values at a designated wavelength i,j in units of the measurement (in our study nanometer (nm)), and a, b, c, d are rational number coefficients selected by the computer  [18] protocols to prepare the spectra for analysis ( Figure 2). First, we raging technique in order to create an even, nominal pace between the two sides tra. That is, we eventually received 231 wavelengths of the 1222 total wavelengths en, outlier spectra were identified by Cochran's test [19] and removed from the dard normal variate (SNV) test followed in order to correct for multiple scatter were corrected for their additive dispersion effect with a baseline correction and the maximum peak within each spectrum. The procedure resulted in a negative he water absorption band (between 1380-1450 nm). Therefore, these wavelengths further usage and analysis. Finally, the spectra were mean-centered and e performing further analysis on it. Stomatal conductance data was also curated gation treatments, where if sample values were far from the average (as in, more deviation distance), the whole sample was omitted from the analysis. The ANN architecture correlates with a Pearson correlation at 0.7 between the measured and the predicted stomatal conductance. The "error-free" percentage of the model on the test set had 0.82 confidence, decreasing by 0.1 units from the calibration set, which implied the strength of such a model.

Discussion
The NDI technique pinpoints four regions that may be connected to stomatal conductance. The first region "red-edge" [39] has been long known to be related to plant stress spectral range and has been shown to relate to evapotranspiration in general [12]. It happens because this region in the spectrum is layered with both physiological and chemical processes. On one hand, it borders the absorption of the photosynthetic pigments Chl a and Chl b [40] and, on the other hand, it reflects the photosynthetic activity of the crop [41]. Therefore, it is expected that, with a higher wellbeing, the magnitude of the difference between two plateaus in the reflectance spectrum ( Figure 3A) will increase [42].
The second range relates to the reflectance of the mesophyll tissue of the plant and encompasses many different remotely-sensed traits such as nitrogen concentration [43], pest response-related indices [44], and disease-related indices [45]. The third region relates to water content and cellular structures, such as lignin and cellulose, which are part of the water transfer vessel network [46]. The fourth region relates to starch molecules, which relate indirectly to transpiration, in that it is being synthesized as a transient product of photosynthesis within the leaf's tissues. Thus, it can be argued that, with higher stomatal conductance, more starch is created and, hence, the relationship between spectral properties and chemical activity obtains an increased correlation [47,48]. When generalizing the NDI technique, the horizontal hot spot correlations was in the red-edge region, between 639 nm and 800 nm, with an additional vertical hot spot region spanning 740-1047 nm and 1455-1661 nm. This generalization corroborates the importance of the red-edge region as found in the private case of NDI when all coefficients equaled one.
When inspecting the optimal random forest architecture ( Figure 5G, bold black arrow), the selected architecture included 23 flagged wavelengths (corresponds to 10% of the wavelengths in the dataset, (Table 3)) [49]. The Fraunhofer oxygen absorbance band O 2 A region was selected as the most important region to detect stomatal conductance in the spectra. This region overlaps the fluorescence emission during photosynthesis, which affects the magnitude of the reflected spectrum [41]. It corroborates the dependence of stomatal conductance on photosynthetic activity. Interestingly, the RF algorithm also succeeded in pinpointing the fact that lignin is a very important feature for detection of stomatal conductance, even more than water absorption bands. Mutations in lignin-synthesizing enzymes have been shown to lower the turgor pressure of the plant and generally decrease stomatal conductance and transpiration, thus corroborating this finding as well [68]. Finally, the RF algorithm succeeded in eliciting diseases and pathogen response wavelengths, which are related to stomatal conductance activity. The relationship between plant-pathogen interaction and stomatal conductance was highlighted in past studies, e.g., downy mildew on cucumber leaves [69], carbon starvation in poplar stems caused by various fungi [70], and transpiration decrease in bean by anthracnose [56]. Although these studies corroborate our findings, it should be noted that the pathogen or the disease affected mostly photosynthesis within the plant, where attenuation in stomatal conductance is a secondary effect of the reaction of the plant to the stress in order to sustain photosynthetic activity.
In searching for a positive correlation between the important wavelengths and the actual stomatal conductance, a PLS-R model and an ANN were selected. In the PLS-R constructed model, each component was found to describe a different relation within the spectrum (Figure 6B), where: (a) factor 1 related to tissue structural wavelengths (76%); (b) factors 2-3 related to the water content in the plant (5%); (c) factor 4 related to the oxygen absorption band in the reflectance spectrum (2%). However, as the PLS-R model succeeded in relating the physiological phenomenon to the spectrum only moderately, it implies that a PLS-R model can be used to detect a stressed plant because it differentiates between the maximum and the minimum irrigation treatments well, but it cannot be used to detect more subtle attenuations in the data. Therefore, a more robust statistical engine is required in order to take into account the correlation between the various wavelengths isolated, a task which was eventually achieved by the ANN algorithm in this study.
There are two main sources of error in this study that need to be addressed in future studies. First, the ground truth technique for measurement of stomatal conductance differed in the two seasons (AP4, Delta-t, UK in 2018 and Licor 6800 photosynthesis system, Li-Cor Biosciences, IL, USA in 2019) ().
The authors suggest creating a calibration factor between the two units on a pre-determined gradient irrigated plot where both units are used. Secondly, the two most irrigated treatments exhibited a mixed reaction in their stomatal conductance ( Figure 4A). Although this behavior is expected, it affects the model's construction by the artificial neural network, as it is sensitive to noise in the dataset [71]. It can explain the underestimation of stomatal conductance at the 400-600 mmole H 2 O m −2 s −1 (Figure 7). In order to avoid this in the future, the authors suggest a wider irrigation treatment gradient when calibrating a regression model. An additional step that can be taken is to use herbicides that attenuate the plant's water budget artificially, thus a synthetic effect of stomatal conductance is achieved as well as better control of the physiological phenomenon.
While this study introduces advanced statistical algorithms to relate physiological activity such as stomatal conductance to spectral measurements, this study is limited to the meteorological conditions of Israel and cotton plants only. Several steps are needed in order to generalize the study and validate it in different meteorological conditions: (a) corroborate the measured behavior of the crop with radiative transfer models that account for atmosphere attenuation and its effect on the electromagnetic spectrum, for example, MODTRAN [72], 6S [73], and radiative transfer models of crops, e.g., PROSPECT, which returns a simulated reflectance profile that can be compared with the actual measurements [10], and SCOPE [74] for a simulation of photosynthesis in the given conditions; (b) develop additional formulation that takes into account the geometry of the crop such that, for example, the leaf angle distribution is taken into account [11]; and finally, (c) the authors suggest repeating this study on a large scale for various crops in different phenological stages in order to isolate important key features in the reflectance spectrum that are shared across species and that can improve prediction of stomatal conductance.
Being able to define stomatal conductance at 20 wavelengths is a promising step towards direct plant physiology assessment with remote sensing techniques. However, both the equipment needed for SWIR measurement, which is expensive, as well as the computer power, which is needed in order to compute a 20 wavelength model, render this study as purely academic. An additional analysis is required in order to construct a simpler model that relates fewer wavelengths to physiological activity.

Conclusions
This study introduced an improvement to current chemometric methodologies and intended to detect and extract physiological activity from reflectance spectra. Specifically, this study found the regions that responded mostly to irrigation treatments in cotton plants were photosynthesis and red-edge in the NIR and lignin and water-content in the SWIR spectral regions. When generalizing coefficient of determination contour maps and selecting bi-wavelengths indices, we showed that a three-fold improvement was achieved in correlating the spectral information to stomatal conductance when compared to the classic approach. However, only two wavelengths were not enough to relate spectral measurements to stomatal conductance due to the diverse biophysical and biochemical processes that affect the conductance. Random forest of regression trees isolated 23 out of the overall 231 wavelengths. These related to water vessel structural constituents, photosynthesis, and plant health wavelengths. These wavelengths revealed new traits that relate to stomatal conductance and are important in designing a remote sensing model. They also corroborated past findings about known wavelengths and their relation to stomatal conductance, hence a physiological activity. An artificial neural network succeeded in relating the spectral data within the 23 wavelengths to stomatal conductance at R 2 = 0.54 with 82% error free confidence on the test sub-set, therefore validating the RF findings. While this study showed the advantage of using machine learning algorithms in identifying direct physiological activity such as stomatal conductance, it should be noted that this study was only the first step towards a robust general remote sensing model of stomatal conductance and plant physiological wellbeing. The examination of the effects of different crop varieties and types is needed in order to generalize the relationship between spectral and stomatal conductance measurements. It will contribute to the understanding of which wavelengths are related to all cases and why there are differences in spectral measurements when extracting information across species on stomatal conductance specifically and plant physiology in general.