Estimating Processing Tomato Water Consumption, Leaf Area Index, and Height Using Sentinel-2 and VENµS Imagery

: Crop monitoring throughout the growing season is key for optimized agricultural pro-duction. Satellite remote sensing is a useful tool for estimating crop variables, yet continuous high spatial resolution earth observations are often interrupted by clouds. This paper demonstrates over-coming this limitation by combining observations from two public-domain spaceborne optical sensors. Ground measurements were conducted in the Hula Valley, Israel, over four growing seasons to monitor the development of processing tomato. These measurements included continuous water consumption measurements using an eddy-covariance tower from which the crop coefficient (K c ) was calculated and measurements of Leaf Area Index (LAI) and crop height. Satellite imagery acquired by Sentinel-2 and VENµS was used to derive vegetation indices and model K c , LAI, and crop height. The conjoint use of Sentinel-2 and VENµS imagery facilitated accurate estimation of K c (R 2 = 0.82, RMSE = 0.09), LAI (R 2 = 0.79, RMSE = 1.2), and crop height (R 2 = 0.81, RMSE = 7 cm). Additionally, our empirical models for LAI estimation were found to perform better than the SNAP biophysical processor (R 2 = 0.53, RMSE = 2.3). Accordingly, Sentinel-2 and VENµS imagery was demonstrated to be a viable tool for agricultural monitoring.


Introduction
Agriculture accounts for 70% of global freshwater usage [1,2], and therefore, increasing the agricultural water-use efficiency will improve agricultural sustainability. Where water is a limited resource, optimal water management is vital for food security. Crop coefficient (Kc)-based estimation of crop water consumption is one of the most commonly used irrigation management methods [3,4]. Kc is defined as the ratio between the actual evapotranspiration from a crop field and the environmental evaporative demand [3]. One of Kc estimation's most reliable sources is vegetation indices (VIs) derived from optical remote sensing [5][6][7][8][9][10][11][12][13][14]. Until recently, this method's application was hampered by the insufficient amount of public domain imagery at a high revisit time with fine spatial resolution. Since 2017 the Sentinel-2 constellation consists of two satellites and serves as a reliable satellite imagery source with high spatial (10, 20, or 60 meters; depending on the band) and temporal resolution (5 days). Despite that, in cloudy regions, even such a high temporal resolution might not be sufficient [15]. For example, despite the Sentinel-2 five-day revisit time, no cloud-free images were acquired for one and a half months in February and March 2018 over one of our experimental sites in Israel. Optical imagery from one satellite system could supplement the imagery from another system to address this problem. Previous studies have analyzed the performance of such conjunction of imagery from different platforms, for example, Landsat-7 and Landsat-8 [16], MODIS and Landsat-8 [17], as well as Landsat-8 and Sentinel-2 [18][19][20][21], and finally, Landsat-7, Landsat-8, and Sentinel-2 combined [22]. Similarly, the present study exploits the possibility of conjoint use of imagery acquired by the Sentinel-2 and the new Vegetation and Environment monitoring on a New MicroSatellite (VENµS) satellite, which has similar spectral bands in the visual, near infrared spectral region, and a 5-10 m spatial resolution (depending on the Collection) as Sentinel-2 in addition to a very high temporal resolution of two days [23].
Tomatoes are grown in many regions around the world. Previously, several studies were devoted to estimating tomato Kc based on lysimeters [24,25] or eddy covariance measurements [26,27] without the correlation to the satellite remote sensing data. Another approach previously used a mechanistic crop model to derive the crop evapotranspiration and correlate it with optical remote sensing data. In this way, previous work [28] used the EPIC model [29], which, in turn, used variables derived from Sentinel-2 imagery.
Additionally, satellite imagery was previously used to estimate other vegetation variables such as LAI and height [11,[30][31][32][33][34]. Much like with Kc, VIs are good surrogates for other crop variables since there are similarities in the temporal change dynamics of VIs with LAI and height [35,36]. LAI is a good proxy of the vegetation state [37][38][39] and a good yield predictor [40]. Similarly, vegetation height estimation is useful for crop management [41]. Therefore accurate estimations of LAI and height from satellite imagery are desired.
Recently, the use of machine learning algorithms has become widespread in remote sensing. In the present study, the LAI biophysical processor [42] implemented in the ESA SNAP (Sentinel Application Platform) 7.0 software (http://step.esa.int/main/download/snap-download/, accessed on 21 February 2021) was tested. The LAI biophysical processor is a "black-box" module developed for Sentinel-2 imagery that cannot currently be used with other imagery.
Therefore, this study's overarching aim was to derive empirical models to estimate vegetation variables based on a combined time series of spaceborne optical imagery from VENµS and Sentinel-2 and field measurements. Specifically, the goal was to develop reliable Kc, LAI, and height estimation models for processing tomato based on Sentinel-2 and VENµS imagery.

Test Sites and Field Measurements
The field data used in this study were collected during four experiments in commercial processing tomato fields in the Hula Valley, Israel ( Figure 1, Table 1). Two experiments took place in Gadash farm in 2018 and 2019, and two more experiments were conducted in Kibutz Gadot in 2019 and 2020. LAI was measured by a SunScan Canopy Analysis System-SS1 manufactured by Delta-T Company (Cambridge, UK) during the two experiments conducted in 2019 and one experiment conducted in 2020. The SunScan is a widely used, accurate, nondestructive LAI measurement system that was successfully employed in many previous studies [31,43,44]. Plant height was measured using a measuring tape during all four experiments conducted in 2018-2020. Each LAI and vegetation height value used in the empirical modelling presented here is an average value of at least 30 field measurements. Both LAI and vegetation height were measured throughout the growing seasons; therefore, they represent the typical range of these variables.
The number of satellite images used for the development of the various models was not uniform because each model was based on the period for which field measurements were available, and therefore, a different number of corresponding satellite images. For example, LAI could not be measured using the SunScan system when the plants were very small, while vegetation height was easily measured at any time. Accordingly, the LAI models were based on shorter time-spans and fewer images than height models.
Each processing tomato field consisted of ridges and furrows. The distance between the rows was 2 m. Even during the vegetation development peak, the plants did not cover the furrows completely; thus, some soil reflectance signal is mixed with vegetation over the entire growing season. This mix of soil and vegetation reflectance hinders the vegetation variables estimation using remote sensing [45]. The Sentinel-2 and VENµS spectral bands used to derive vegetation indices were averaged for an area corresponding with the eddy-covariance footprint. In-field paths and their surrounding area were masked out from analysis polygons to remove bare soil areas and avoid edge effects. These excluded areas consisted of roughly 20% of the overall polygon areas. Therefore, each analysis consisted of either two or four vegetated regions separated by the paths (Figure 1).

Agro-Meteorological Measurements
The reference evapotranspiration, ET0, was calculated based on nearby meteorological stations according to the FAO56 Penman-Monteith method based on meteorological measurements of air temperature, relative humidity, wind speed, and solar irradiance [3]. The actual evapotranspiration (ETc) was measured using eddy covariance systems [26]. Based on these two measurements, the crop coefficient, Kc, was calculated as: Kc = ETc/ET0. Kc is an important variable used to determine the irrigation dose [9]. The resulting Kc time series were smoothed using cubic or second-order splines.

Satellite Imagery
Sentinel-2 is an Earth observation mission and part of the European Space Agency (ESA) Copernicus program. It includes two satellites, each equipped with a Multi-Spectral Instrument (MSI), namely, Sentinel-2A (launched June 2015) and Sentinel-2B (launched March 2017). VENµS is a joint satellite mission of the Israeli and French space agencies (ISA and CNES) launched in August 2017. VENµS has a two-day revisit time over Israel and a multispectral camera with 12 narrow spectral bands in the range of 415-910 nm [46]. VENµS and Sentinel-2 produce 10 and 12-bit radiometric data, respectively. The radiometric correction procedure of VENµS imagery was updated in 2020. The imagery acquired before the update is known as Collection 1; the imagery acquired after the update is known as Collection 2. VENµS captures s imagery with a spatial resolution of 10 m. Sentinel-2 RGB and NIR bands also have a spatial resolution of 10 m, and other bands are coarser: narrow NIR, SWIR, and red edge bands, 20 m; coastal aerosol, water vapour, and SWIR-cirrus bands used mostly for atmospheric correction, 60 m. Atmospherically corrected reflectance products from both sensors were used in this analysis. Level-2 VENµS products, initially distributed at 10 m spatial resolution, were later distributed at a resolution of 5 m when an updated processing procedure was initiated in 2020. This product was used for the analysis of the 2020 experiment in Gadot. Sentinel-2 level-2A data were obtained from the ESA Copernicus Open Access Hub website (https://scihub.copernicus.eu/dhus/#/home, accessed on 21 February 2021). VENµS level-2 products were obtained from the Israel VENµS website maintained by Ben-Gurion University of the Negev (https://venus.bgu.ac.il/venus/, accessed on 21 February 2021). Table 2 lists the overlapping spectral bands of the Sentinel-2 and VENµS sensors used in this study to derive vegetation indices. The LAI and Kc estimation models were derived based on three seasons, and crop height models were based on four seasons. An inventory of the Sentinel-2 and VENµS images used in the present study can be found in Table 3, alongside the number of LAI and height measurements taken during each season and used for model derivation.

Vegetation Indices and Model Validation
All Sentinel-2 and VENµS bands were resampled to 10 m spatial resolution. After that, thirteen vegetation indices (Appendix A) were derived based on the Sentinel-2 and VENµS imagery, including transformed VENµS imagery that utilised a corrective transformation (Table 4) derived for collection 1 VENµS imagery [23]. Since the radiometric processing of VENµS was improved in collection 2, the applicability of the transformation functions to the re-calibrated VENµS imagery was studied by comparing the performance of models based on the imagery transformed for all seasons against the models based on transformed imagery for 2018-2019 seasons (collection 1) and not transformed for 2020 (collection 2). The performance of the former was found to be better than the latter. Therefore, the transformed VENµS imagery models were applied to all seasons. Overall, three types of tomato estimation models were derived: models based on Sentinel-2; models based on Sentinel-2/non-transformed VENµS; models based on the Sentinel-2/transformed VENµS imagery. Hereafter the combined Sentinel-2/transformed VENµS models will be referred to as S2/VT, and combined Sentinel-2/non-transformed VENµS models will be referred to as S2/VNT. Linear regression models were derived for the time series of field-measured Kc, LAI, height, and each spectral index time series. Each model was based on all available field measurements of each vegetation variable collected during all seasons when the variable was measured. For every model, the R 2 and root mean square error (RMSE) values were calculated. RMSE was calculated for each model based on all available data and also for each field experiment separately. In addition to vegetation index-based models, an LAI estimation from the ESA SNAP 7.0 biophysical processor for Sentinel-2 imagery was also produced [42].
The S2/VT and S2/VNT models were compared, and the Steiger variation [47] of the two-tailed Fisher's Z-score tests [48] was performed to determine whether the difference in the models' R 2 is significant (α ≤ 0.05). The same test was also performed to determine whether the difference in R 2 of the LAI Biophyscal processor and DVI was significant.
The field-measured processing tomatoes LAI and height measured in Gadash 2019 and Gadot in 2019 and 2020 were used to calibrate prediction models for Kc, as was done previously [49].  Figure 2E-G shows the three types of the aforementioned variations of the Kc associated with three experiments conducted in Israel. The standard Kc recommendation differs from the measured Kc values. Early in the season, during the crop vegetative development, the standard table recommendation is slightly higher than the measured water consumption. In Gadot 2019, the standard recommendation and measured water consumption are about the same at the peak. In Gadot 2020 and Gadash 2019, the standard recommendation's peak is higher than the measured water consumption. However, from the mid-late season, the measured water consumption drops below the standard recommendation. Interestingly, in Gadash 2019, the crop height and LAI and the Kc were lower compared to the other seasons. Moreover, the changes in LAI and height in Gadash 2019 were different compared to other seasons. These discrepancies in behavior between tomato variables and differences in the variables' values from season to season demonstrate the variance in crop development and water consumption between seasons. Therefore, real-time estimations of those variables are advantageous over the use of standard tables. GEMI and WDVI were found to be the best VIs for the tomato Kc, crop height, and LAI estimation. These results repeated in all three types of models: Sentinel-2-based, S2/VNT, and S2/VT. Tables 5-7 show Sentinel-2, S2/VNT, and S2/VT-based Kc, crop height, and LAI estimation models based on the five best-performing VIs: DVI, GEMI, WDVI, SAVI, and MSAVI. The best combined Sentinel-2/VENµS models in the present study are presented in Figure 3. The data points in Figure 3 are not clustered by sensors or experiments, which is indicative of the models' generality. Therefore, both sensors used in the study can be employed interchangeably. The tomato Kc, height, and LAI estimation models' performance is based on eight other VIs (NDVI, MTCI, IPVI, IRECI, S2REP, REIP, GNDVI, and TNDVI), which can be found in Appendices B-D. Table 5 shows that the RMSE of LAI derived from the biophysical processor is higher and the R 2 is lower than VIs such as GEMI, DVI, WDVI, SAVI, and MSAVI. The biophysical processor's R 2 was found significantly lower than the R 2 of DVI (p = 0.016). It was found that the majority of S2/VT and S2/VNT models do not present significant differences in performance and that the transformation of VENµS imagery is mostly beneficial for the red edge VIs such as MTCI and S2REP (Appendix E). Table 8 shows the difference in performance between S2/VT models and S2/VNT models of the best performing VIs (DVI, GEMI, WDVI, SAVI, and MSAVI). Appendix E shows the difference in performance between S2/VT models and S2/VNT models for eight additional VIs (NDVI, MTCI, IPVI, IRECI, S2REP, REIP, NDVI, and TNDVI).    Figure 4A shows LAI and height measurements (in dm; to fit them to the same Yaxis) recorded during three field campaigns in 2019 and 2020. Interestingly, in Gadot 2019, height continued to increase in the middle of the season, while LAI has already started to decrease. In the other fields measured in this study, LAI and height varied simultaneously. Figure 4B   The performance of processing tomato height and LAI-based Kc estimation models using field measurements is shown in Table 9. Table 9. Kc prediction models based on field measurements of processing tomatoes height and LAI.

Discussion
The field experiments conducted in Israel in 2018-2020 showed that Kc, LAI, and crop height in processing tomato differ from season to season but can be estimated correctly in near-real-time from satellite remote sensing imagery. Consequently, agricultural decisions, including the irrigation dose determination, can rely on remote sensing data rather than standard tabular recommendations until late in the season. During the last stage of the season, deficit irrigation is applied according to the percentage of ripe fruit (the ratio of red to green tomatoes on the plant) to delay ripening or expedite it according to the desired harvest schedule [50]. Thus, the irrigation dose at the end of the season cannot be estimated using the remote sensing approach described here.
The field-measured Kc in this study yielded high correlations with VIs from Sentinel-2 and VENµS. Consequently, this study paves the way for more precise Kc, LAI, and crop height estimations on a local and global scale based on the freely available optical satellite imagery. These crop variable estimations could be used for better irrigation and fertilization management [51], as well as for early detection of crop disease [52,53], waterlogging [54,55], pest management, and biological control [56].
This study's most important result was the demonstration of effectively joining Sentinel-2, and VENµS imagery for agricultural monitoring suggested before the launch of those missions [38]. This was possible because of the close resemblance of Sentinel-2 and VENµS spectral response functions and a good radiometric and atmospheric correction. Application of corrective transformation functions [23] improved the performance of VIs based on the red edge bands (MTCI, S2REP, and REIP), while for the other VIs, the transformation was found unnecessary or provided only marginal performance improvement.
Many VIs showed good Kc estimation performance; the best performing Kc estimation was achieved with the GEMI S2/VNT model (R 2 = 0.82, RMSE = 0.09). In an earlier study, a canopy cover-based Kc estimation model achieved R 2 = 0.96 [27]. In that study, the canopy cover was calculated using cameras installed in the field. Unlike this approach that relies on in-field sensors, the approach suggested in this paper, based on satellite remote sensing, facilitates the estimation of vegetation variables over more extensive areas at a low cost. This study shows that Kc estimation from optical satellite remote sensing can serve as a reliable source for irrigation decisions and potentially for other agricultural activities throughout the whole growing season of processing tomato. The best performing LAI estimation models showed promising results (S2/VT WDVI LAI estimation model: R 2 = 0.79, RMSE = 1.2). This result agrees with a previous study that found WDVI, which takes soil reflectance into account, as a good indicator of LAI [57]. In comparison to the newly-obtained processing tomato LAI models, multi-crop models derived in previous studies demonstrated lower performance, e.g., R 2 = 0.62 [58], R 2 = 0.66 [59], R 2 = 0.72 [60]. A tomato LAI model from previous work [28] showed a lower coefficient of determination (R 2 = 0.69) and lower RMSE = 0.56 compared to this study. However, this model was based on only four days of field measurements. Moreover, that work [28] did not include LAI measurements in the final stage of a growing season, while the LAI models in the present study were based on three full growing seasons. Consequently, the processing tomato LAI estimation models developed in the present study are suitable for general use in precision agriculture applications throughout the growing season. Additionally to LAI estimation based upon the VIs, the performance of the SNAP biophysical processor LAI estimation algorithm was studied (R 2 = 0.53, RMSE = 2.3) and found to be significantly less accurate compared to the empirical model based on DVI, which was found to be the most accurate Sentinel-2-based VI for LAI estimation.
Similar to Kc and LAI estimation models, the tomato height estimation models were found to perform well throughout the processing tomato growing season. The S2/VT WDVI-based height estimation model (R 2 = 0.81, RMSE = 7 cm) was found to be the best, and this approach shows great promise for agricultural crop monitoring. The obtained results confirmed the previously found conclusion that WDVI is a well-suited VI for crop LAI and height estimations [33].
Kc, LAI, and height estimation models based solely on Sentinel-2 data were as accurate as the combined Sentinel-2/VENµS models. Subsequently, a pooled time-series of imagery from both sensors increases the available satellite imagery's temporal resolution. In cloudy regions, either sensor could fill gaps in the acquisitions of the other, and either sensor can efficiently monitor crop development when imagery from the other sensor is not available. For example, during both experiments in 2019, many VENµS images filled in a long gap in Sentinel-2 data in April-May, and one Sentinel-2 image filled a gap in VENµS images in May-June.
Additionally to the Kc estimation based on the remote sensing data, Kc estimation models based on the field measured LAI and height were derived. These models' performance was similar to the remote sensing-based models and might be used on the local scale in the absence of remote sensing imagery. The Kc-height model is of particular interest from a practical viewpoint since farmers can easily and routinely take plant height measurements.
While this study provided useful results from thirteen VIs (including VIs based on the red edge bands and soil adjusted VIs) to estimate Kc, LAI, and height in the processing tomato using Sentinel-2 and VENµS imagery, there is merit in future studies on other crops. Future efforts could follow the procedure suggested in this paper to empirically calibrate and test prediction models for different indices and identify those that achieve the best performance. Studies based on two or more different sensors should make sure to perform a radiometric calibration between sensors.

Conclusions
This work demonstrates the conjoint use of Sentinel-2 and VENµS imagery for estimating Kc, LAI, and height of processing tomato. It was found that red edge VIs should be based on Sentinel-2 and transformed VENµS imagery. At the same time, other VIs can be derived directly from imagery obtained by both systems, and no corrective transformation is required to match the two sensors. In addition, models based solely on Sentinel-2 showed similar performance as the joint Sentinel-2 and VENµS imagery models. The Kc, LAI, and height estimation models derived empirically using field measurements show good performance and are ready for application. The LAI estimation performance from the SNAP biophysical processor was also studied and found inferior to the VI-based LAI estimation models. The irrigation in the early and middle parts of the processing tomato growing season can rely on remote sensing-based models rather than standard table values to best match the actual crop development and capture within-field variability.   Appendix B Table A2. Performance statistics of newly developed Sentinel-2-based LAI, Height, Kc models, and the performance of the SNAP biophysical processor LAI estimation algorithm. Performance statistics of better performing VIs can be found in Table 5.  Table A5. Difference between performance statistics of S2/VT and S2/VNT-based LAI, Height, Kc models. If R 2 is positive and RMSE is negative, it means that this parameter performance of the combined S2/VT model is higher than the equal parameter of the S2/VNT model. Significant differences are marked with (*). Performance statistics of better performing VIs can be found in Table 8.