Leaf water potential of field crops estimated using NDVI in ground-based remote sensing—opportunities to increase prediction precision

Application and construction of a ground-based phenotyping tool

To assist in field crop phenotyping efficiently and in a consistent manner, we developed a portable proximal sensing cart based on the model of White & Conley (2013). We added wheel assembly improvements to the cart frame, facilitating easy maneuverability by one person in field conditions, as well as adding nested threaded clamps to adjust the cart width, allowing crop measurement with different row spaces (Fig. 1). Attached to the cart are three highly accurate canopy temperature sensors, one ultrasonic crop displacement sensor, one active light NDVI sensor, and a marine grade GPS receiver (A–F in Fig. 2). We also constructed a backpack frame that allows selected measurement of NDVI in tall crops, such as corn, sorghum and sesame (G in Fig. 2). The ACS-430 Crop Circle sensor measures canopy reflectance from three optical channels (red–670 nm, red-edge–730 nm, and NIR–780 nm), and NDVI was calculated as (Gong et al., 2003):

(1) $$\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}={\displaystyle \frac{{\rho }_{780}-{\rho }_{670}}{{\rho }_{780}+{\rho }_{670}},}$$

where ρ₇₈₀ and ρ₆₇₀ represent reflectance at NIR and red band, respectively. To allow us standardized accurate customizable environmental data sampling, we utilized a Campbell Scientific CR3000 voltage data recorder including a CF memory card and with its CRBasic language data collection program execution functions. We subsequently acquired precise raw sensor potentials and their associated GPS coordinate strings, recorded to table structured digital records of 200 ms granularity (companies and products are not endorsed by the authors and are presented for educational purposes only). A full list of the CRBasic program driving the CR3000 datalogger is available as Computer Code S1.

Field data acquisition and collection

Field crops of wheat, corn, and cotton were grown under deficit- and full irrigation regimes in a center pivot field at the Texas A&M AgirLife Research and Extension Center at Uvalde, TX in the 2017–2018 growing season. The full irrigation was determined based on full-replenishment of the crop evapo-transpiration (ET) measured at the Research Center facility, while deficit irrigation was determined as 60% of the full irrigation. Diurnal changes in canopy NDVI and leaf water potentials for the three crops were measured on five clear days (see Table 1 for further information). Canopy NDVI for wheat and cotton was measured using the sensing cart, while that of corn was measured using the backpack sensing frame. Leaf water potential was measured using a PMS-615 Pressure Chamber (PMS Instrument Company, Albany, OR, USA), with the high pressure generated by pure nitrogen gas. For wheat, one transect encompassing a full and a deficit irrigation field was used, while for corn and cotton, a circular transect within a quarter of the center pivot field was used for the measurement. Within the circular transect, six equal divisions were delineated—each at least 30-m long, covering 15° angle in the circle and randomly assigned to receive either deficit or full irrigation treatment, with the six sections designated as six plots, i.e., deficit-1, full-1, full-2, deficit-2, deficit-3 and full-3, starting at the west side and ending at the east side of the circular path. On each of the five days selected for making field measurement, NDVI was scanned within the transect repeatedly once every hour from early morning to late afternoon. On each day, the datalogger was turned on before start of the first scan, and kept on until completing the day’s measurement in late afternoon. Each hour when NDVI was measured, 4–6 healthy and representative leaves of corn, cotton or wheat were chosen for measuring leaf water potentials (LWP) under field conditions (Fig. 3). For wheat, the flag leaves were selected to make the LWP measurement; for corn and cotton, the fourth leaf counting down from the fully expanded youngest leaf were used for the LWP measurement. To minimize water loss from the leaves enclosed in the chamber during pressurization, one moistened, 15-cm diameter filter paper was lined against the interior of the chamber wall and maintained moistened throughout the measurements.

To document changes in green biomass during the growing season, leaf area index (LAI) for each of the crops was measured using a LI-3100 Area Meter (Li-Cor Inc., Lincoln, NE, USA) through destructive sampling. The measurement was made weekly or once every 2 weeks during the periods when there were significant changes in green biomass. For wheat, three random shoot samples were collected by harvesting a 20-cm long of a representative row from both deficit and full irrigation segments of the transect. For corn and cotton, one representative plant per plot was harvested to measure total leaf area on each measurement period. LAI was estimated based on total green leaf area per plant, row spacing and population density (shown in Table 1).

Data processing and statistical analysis

In this study, the data analysis is focused on the measured NDVI and its relationship with LWP. Over five days, about 1.3 million lines of data were recorded using the ground-based phenotyping tools on three crops. Some of the data lines were recorded while the sensor was in rest awaiting next round of scans. This was done in order to avoid the datalogger being turned off and back on repeatedly during the day, causing inadvertent errors. These unwanted readings during idle period, as well as other unwanted readings, were deleted during data processing according to recorded GPS coordinates, as well as the Coordinated Universal Time (UTC), in order to match the sensor data with field plots/treatments and time of measurement on different days. This initial data processing was done in Minitab (Version 17.3.1, 2016; Minitab Inc., State College, PA, USA). The datalogger recorded UTC time was translated into the Central Standard Time (CST) during the local day-light-savings time period by the relation CST = UTC − 5. The summarized data of NDVI and LWP were then further analyzed as described below.

Two methods were employed to conduct the least-squares regression analysis describing the relationship between measured values of NDVI (x) and leaf water potential (y) in the form

(2) $$Y=mX+c,$$

where Y and X represent adjusted/calculated values corresponding to measured y and x, m and c represent the slope and y-intercept of the model.

• 1. The first method is default to most statistical software packages where a common uncertainty in y is assumed, and the measurement errors in x are considered negligible as compared to those of y. Specifically, the regression analysis was conducted using GraphPad Prism 6 (Version 6.07 for Windows, 2015, GraphPad Software, San Diego, CA, USA). Statistical differences among and between the slopes and the y-intercepts of the linear regressions representing different treatments were compared using ANCOVA in Prism 6. The uncertainties in the coefficients m and c were calculated according to Taylor (1982) and implemented using a Minitab macro ‘lsq.mac’ (see Computer Code S2 and Dataset S1). The maximum relative errors in predicted LWP given a specific NDVI observation was estimated according to rules of error propagation (Taylor, 1982):

(3) $${\displaystyle \frac{\delta y}{|\overline{y}|}\le {\displaystyle \frac{1}{|\overline{y}|}\left[\left|{\displaystyle \frac{\mathrm{\partial }y}{\mathrm{\partial }m}}\right|\delta m+\left|{\displaystyle \frac{\mathrm{\partial }y}{\mathrm{\partial }c}}\right|\delta c\right]={\displaystyle \frac{1}{|\overline{y}|}(|\overline{x}|\delta m+\delta c),}}}$$

• where δy, δm, and δc are estimated uncertainties in y, m, and c, and $\overline{x}$ and $\overline{y}$ are the mean values of x and y, respectively. The first method is subsequently referred to as the method of OLS (Ordinary Least Squares).

• 2. The second method is based on the principle of weighted linear least-squares considering measurement uncertainties in both coordinates (York, 1966). In this method, the scheme of Deming (1943) was used to give different data points different weighting factors ω (x_i) and ω (y_i) that were defined as the inverse squares of the uncertainties δ(x_i) and δ(y_i) (i.e., ω (x_i) = 1/δ² (x_i) and ω (y_i) = 1/δ² (y_i)), and the best values of m and c of Eq. (2) were identified by minimizing the sum of the weighted squared residuals S as:

(4) $$S=\sum _{i=1}^N\left[\omega (x_i)(x_i-X_i{)}^2+\omega (y_i)(y_i-Y_i{)}^2\right],$$

• where x_i and y_i are observed values and X_i and Y_i and calculated ones with i running from 1 to N. York (1966), McIntyre et al. (1966) and Williamson (1968) developed the full solution of the dual-uncertainties least-squares problem, which has received wide recognition in geosciences and physical sciences (Cantrell, 2008). However, the same method appears to be less-utilized by biologists who may equally be faced with the same type of statistics problems as do geoscientists, physicists and chemists. For this reason, it is necessary to summarize the main steps for finding the model parameters and their uncertainties in this second method, based on equations of York (1969) and Reed (1992).

• With the “overall weight” for the ith data point defined as

(5) $$Z_i={\displaystyle \frac{\omega (x_i)\omega (y_i)}{m^2\omega (y_i)+\omega (x_i)},}$$

• Equation (4) can be written as

(6) $$S=\sum _{i=1}^N{Z}_i(y_i-m{x}_i-c{)}^2.$$

• The pursuit of minimization of Eq. (6) leads to a “least-squares cubic” in the form of

(7) $$m^3\sum _{i=1}^N{\displaystyle \frac{Z_i^2{U}_i^2}{\omega (x_i)}-2m^2\sum _{i=1}^N{\displaystyle \frac{Z_i^2{U}_i{V}_i}{\omega (x_i)}-m\left[\sum _{i=1}^N{Z}_i{U}_i^2-\sum _{i=1}^N{\displaystyle \frac{Z_i^2{V}_i^2}{\omega (x_i)}}\right]+\sum _{i=1}^N{Z}_i{U}_i{V}_i=0,}}$$

• where Z_i is from Eq. (5), $U_i=x_i-{\sum }_{i=1}^N{Z}_i{x}_i({\sum }_{i=1}^N{Z}_i{)}^{-1}$ and $V_i=y_i-{\sum }_{i=1}^N{Z}_i{y}_i({\sum }_{i=1}^N{Z}_i{)}^{-1}$. Equation (7) can be solved for m to obtain the best value of the slope. However, since Z_i, U_i and V_i are themselves all functions of m, the equation is actually a pseudo-least-squares cubic, and the best value of m is obtained through successive iterations starting from a crude estimation, such as from an unweighted least-squares fitting. A least-squares cubic similar to Eq. (7) with the associated iterative solution was developed independently by McIntyre et al. (1966) in analyzing the Rb-Sr isochrons from geological samples. Equation (7) may also be formulated in pseudo-quadratic- (York, 1969; Reed, 1992) or pseudo-linear form (Williamson, 1968; Reed, 2015), but they all need to be solved iteratively, similar to the cubic case. Using multiple data sets, Reed (1989) demonstrated that finding the correct solution to Eq. (7) can sometimes be tricky and confusing, and provided a direct searching method to quickly “pin down” the true value of m starting from a relatively good “seed” value. Once the best value of m is obtained, the like value of c can be found from

(8) $$c=\sum _{i=1}^N{Z}_i{y}_i{\left(\sum _{i=1}^N{Z}_i\right)}^{-1}-m\sum _{i=1}^N{Z}_i{x}_i{\left(\sum _{i=1}^N{Z}_i\right)}^{-1}.$$

• Built on the work of York (1966), Reed (1992) provided a corrected expression of the variance of m:

(9) $${\delta }_m^2={\displaystyle \frac{S}{N-2}\sum _{i=1}^N\left[{\displaystyle \frac{1}{\omega (y_i)}{\left({\displaystyle \frac{\mathrm{\partial }m}{\mathrm{\partial }{y}_i}}\right)}^2+{\displaystyle \frac{1}{\omega (x_i)}{\left({\displaystyle \frac{\mathrm{\partial }m}{\mathrm{\partial }{x}_i}}\right)}^2}}\right],}$$

• where S is from Eq. (6), ${\displaystyle \frac{\mathrm{\partial }m}{\mathrm{\partial }{y}_i}}$ and ${\displaystyle \frac{\mathrm{\partial }m}{\mathrm{\partial }{x}_i}}$ can be computed from x_i, y_i, ω (x_i), ω (y_i) and the best value of m (see Appendix of Reed (1992) for the rather involved explicit expressions). The expression for ${\delta }_c^2$ (variance of c) is identical to that of ${\delta }_m^2$ but with derivative of c in place of m in Eq. (9).

• Several researchers (McIntyre et al., 1966; York, 1969; Lybanon, 1984; Reed, 1989; Press et al., 1992), and a few others as discussed by Cantrell (2008), have developed computer programs (mainly in Fortran) to implement the algorithms of finding both the best parameter values (i.e., m and c) and their variances for the least-squares fitting problem with errors in both coordinates. Notably, Cantrell (2008) developed an Excel spreadsheet program to implement the Williamson-York method of bivariate linear fit. Reed (2010) translated his Fortran program into an easy-to-use Excel spreadsheet, and later updated it (Reed, 2015) to account for the situation with correlated x-y uncertainties, one that occurs widely in the study of isotope ratios in geosciences (York, 1969; York et al., 2004). However, in this paper, we only consider the situation with non-correlated x-y uncertainties, because the interested quantities (i.e., NDVI and LWP) were derived from measurements using different types of equipment independently, and thus their uncertainties are assumed to be uncorrelated.

• Our second method relied on the use of the Excel spreadsheet ‘LLS(SIGMAS).xls’ of Reed (2010). Specifically, the Excel function ‘Goal seek’ was used to obtain the best estimation of coefficient m that ensured a function g(m) (such as Eq. (7)) taking the value of zero. Then the best value of m was used to determine the best value of c according to Eq. (8). To facilitate the determination of m, it was important to first use the unweighted least-squares routine builtin in ‘LLS(SIGMAS).xls’ to get a seed value of m which was then input to the ‘Goal seek’ function to find the best value of m (Reed, 2010). Since the values in both coordinates are adjusted in this method, the uncertainties δm, and δc can be calculated either at the observed points or at the calculated/adjusted points (Reed, 2010). In ‘LLS(SIGMAS).xls’, this is achieved using a switch: to evaluate the derivatives at the observed points, set PASS = 0, or to evaluate them at the calculated/adjusted points, set PASS = 1. We used the second option (i.e., at the adjusted points) but the differences in results using the two options were generally small for well-correlated data. Equipped with the best estimates of m and c and their variances (see Eq. (9)), which are part of the outputs of the above spreadsheet program, the maximum relative errors in predicted LWP using the second method were calculated according to Eq. (3). This method is subsequently referred to as the method of WLS (Weighted Least Squares).

We also used statistical resampling to investigate the effect of measurement errors in LWP on the regression model. In particular, we wanted to see how the reduced sample size of LWP would affect the estimated regression coefficients and their uncertainties. The resampling was done using a Minitab macro ‘sem_a.mac’ (see Computer Code S3 and Dataset S2), in which multiple sets of resampled 2 or 3 replicates of LWP values were drawn randomly (with replacement) from the original 4–6 measurements made during the field surveys. The resampled values of LWP were then used as inputs, along with the originally measured values of NDVI, to ‘lsq.mac’ and ‘LLS(SIGMAS).xls’, in order to obtain the best estimates of the regression coefficients and their uncertainties using the method of OLS and WLS, respectively. Differences in the regression coefficients (and associated uncertainties) between the situation with resampled LWP values and that with the originally measured LWP values were compared using the One-sample t-test available in Minitab.

NDVI in relation to leaf area index and leaf water potential

Average values of LAI of wheat, corn and cotton were significantly lower under deficit irrigation than under full irrigation (p < 0.05, Fig. 4), especially in the mid- to late growth stages when canopy NDVI and leaf water potentials were measured. These trends in LAI paralleled those in measured NDVI, as seen in corn measured on June 7, 2018 (Fig. S1). In the case of the latter, it is evident that the differences between the deficit and full irrigation grew larger with the progression of time during a day. However, the NDVI values of wheat measured on April 12, 2018 exhibit a different trend, one in which the values measured from the deficit irrigation was high than that from the full irrigation, showing an opposite trend with the measured LWP (Fig. 5A). This was an unexpected result, but we did notice that during this day’s measurement there was significant lodging in the wheat field. This happened especially after the most recent irrigation event, possibly related to the Hessian fly infection. The measured NDVI values on April 19, however, show an expected trend, being higher in full than deficit irrigation (Fig. 5B). As a result, except for wheat measured on April 12, a higher LWP was associated with a higher NDVI in the three field crops.

In drawing Fig. 5A, the NDVI values measured at 8 am (NDVI = 0.55 and NDVI = 0.63 for full and deficit irrigation, respectively; see raw data file NDVI_wheat_april_12.xls in the Supplemental Material) were not used, since these values were collected when the Crop Circle sensor was not fully warmed up in the relatively cold morning of April 12, 2018, with an air temperature as low as 15.8 °C at 8 am (see Table 1 and Fig. S2). The peculiar observation that the wheat NDVI under deficit irrigation was higher than that under full irrigation as measured on April 12 (Fig. 5) prompted us to inspect the raw data of canopy reflectance measured at the red and NIR bands. As seen in Fig. 6, the higher wheat NDVI on April 12 in the deficit irrigated plots was due both to a lower reflectance at the red band and the higher reflectance at the NIR band, as compared with the plots under full irrigation (see the definition of NDVI in Eq. (1)). The observation that the wheat NDVI in plots of deficit irrigation on April 19 was higher than that of the full irrigation plots was, however, due only to the higher red band reflectance associated with the deficit irrigation, because the NIR reflecance was similar between the two irrigation regimes (again see Fig. 6).

Another important point is if there was a significant difference in NDVI within a single day. For winter wheat, this was not tested statistically because the NDVI measurement was not replicated. By visual inspection (Figs. 5A, 5B), the diurnal trend of NDVI was unclear in both April 12 and 19, which was in contrast to the clear diurnal pattern of LWP in those two days. For corn and cotton, however, there were clear diurnal changes in NDVI (Figs. 5C, 5D, 5E). One-way ANOVA (n = 3, since each of the NDVI scans for corn or cotton was done with three replicates) indicated that, except for corn under full irrigation on June 11 that exhibited a marginally significant diurnal change (p = 0.056), five other diurnal courses of NDVI for corn/cotton all showed significant or highly significant diurnal trends (p < 0.05).

Since the relationship between LAI and NDVI was strong (Fig. S3), and there was significant differences in LAI values between deficit and full irrigation treatments, we analyzed the relationship between NDVI and LWP separately across the different crop and different irrigation treatments (Fig. 7). For wheat, there was an unclear or negative relationship between NDVI and LWP (Fig. 7A), while for corn and cotton, the relationship was positive. Though in the case of corn under full irrigation on June 11, 2018, the slope of the linear regression was not significantly different from zero (Figs. 7B, 7C). For those statistically significant regression lines shown in Fig. 7, further test indicated that, within corn or cotton, the slopes of the regressions lines are not significantly different, but the y-intercepts are significantly different (Table 2). Despite the low number of data points used in fitting the linear regression models, the relationships are generally strong as evidenced in the high R² values in Table 2.

Uncertainty analysis of NDVI and leaf water potential

Table 3 shows the uncertainties of the regression coefficients for the six significant linear regressions using the OLS method. In the last column of this table, the relative errors of the predicted LWP are shown. Since there was no statistical significance in regression coefficients between two cotton datasets (full vs. deficit irrigation), they were combined to form a larger dataset for further analysis. This is shown in the last row of Table 3. We can see that the prediction errors were large, with the best being 34% when the two cotton datasets were combined. However, when the uncertainties in both coordinates were considered in the least-squares regression analysis using the WLS method, the uncertainties in both coefficients were reduced, and also reduced were the prediction errors (Table 4). Further calculation details of the results in Tables 3 and 4 are shown in Dataset S3.

The difference between the results using two regression methods was further illustrated by the relative errors of the regression coefficients m and c. As seen in Table 5, the WLS method yielded lower relative sizes of the uncertainties for the coefficients, as compared with the OLS method, except for corn measured on June 11. The latter case was likely due to the very high relative errors in measured NDVI for corn on this day, as shown in Table 6.

The 20 sets of resampled LWP values and the corresponding uncertainties for cases of 2 and 3 resampled values of LWP are shown in columns 1–40 of Datasets S4 and S5, respectively. The resultant model coefficients and their uncertainties after these resampled LWP values were used as inputs to the OLS procedure of ‘lsq.mac’ are shown in columns 41–44 of Datasets S4 and S5, respectively. The resultant model coefficients and their uncertainties using the WLS procedure of ‘LLS(SIGMAS).xls’ are shown in Datasets S6 and S7, respectively. The results of statistical resampling indicate that reducing the sample size of LWP to 2 or 3 significantly increased the uncertainties of the estimated coefficients m (0.92 vs. 0.77) and c (0.48 vs. 0.39) as compared with the respective uncertainties when the original 4 replicates of LWP were used as inputs to the WLS method (Fig. 8). Reducing the sample size of LWP also increased the prediction uncertainty of LWP to 0.32 (from 0.27 using the WLS method as seen in Table 4). However, there was no significant difference in the uncertainty of either coefficient when the OLS method was used, as seen by the dashed lines passing through or touching the error bars of the resampled values using the OLS method (Fig. 8). This result was anticipated since no errors from the individual observations of NDVI or LWP were considered in the OLS method.