Open AccessArticle

Seasonal Mapping of Irrigated Winter Wheat Traits in Argentina with a Hybrid Retrieval Workflow Using Sentinel-2 Imagery

Gabriel Caballero

^1,2,*

Alejandro Pezzola

³,

Cristina Winschel

³,

Alejandra Casella

⁴

Paolo Sanchez Angonova

³,

Juan Pablo Rivera-Caicedo

⁵

Katja Berger

^2,6

Jochem Verrelst

and

Jesus Delegido

Agri-Environmental Engineering, Technological University of Uruguay (UTEC), Av. Italia 6201, Montevideo 11500, Uruguay

Image Processing Laboratory (IPL), University of Valencia, C/Catedrático José Beltrán 2, 46980 Paterna, Spain

Remote Sensing and SIG Laboratory, Hilario Ascasubi Agricultural Experimental Station, National Institute of Agricultural Technology (INTA), Hilario Ascasubi 8142, Argentina

⁴

Permanent Observatory of Agro-Ecosystems, Climate and Water Institute-National Agricultural Research Centre (ICyA-CNIA), National Institute of Agricultural Technology (INTA), Nicolás Repetto s/n, Hurlingham 1686, Argentina

⁵

Secretary of Research and Graduate Studies, CONACYT-UAN, Tepic 63155, Mexico

⁶

Mantle Labs GmbH, Grünentorgasse 19/4, 1090 Vienna, Austria

Author to whom correspondence should be addressed.

Remote Sens. 2022, 14(18), 4531; https://doi.org/10.3390/rs14184531

Submission received: 9 August 2022 / Revised: 31 August 2022 / Accepted: 7 September 2022 / Published: 10 September 2022

(This article belongs to the Special Issue Cropland Phenology Monitoring Based on Cloud-Computing Platforms)

Download

Browse Figures

Figure 1
Overview of the Bonaerense Valley of Colorado River study site with test fields for the winter wheat campaign of the year 2020. True colour S2 image (R = B4, G = B3, B = B2) of 27 December 2020. Reference system: WGS84 (EPSG 4326). "> Figure 2
Layout of measurement design for the 2020 campaign at the BVCR study site: three ESUs were defined per wheat paddock and sampling approach for each elementary sampling unit, partly adapted from [<a href="#B69-remotesensing-14-04531" class="html-bibr">69</a>]. Reference system: WGS84 (EPSG 4326). "> Figure 3
Photographic documentation of wheat crop growing period from August to November 2020, BVCR campaign, corresponding to LAI, FVC, C<math display="inline"><semantics> <msub> <mrow/> <mrow> <mi>a</mi> <mi>b</mi> </mrow> </msub> </semantics></math>, and AGFB sampling dates. With (a) seedling stage at 10 August 2020; (b) tillering stage at 4 September 2020; (c) tillering stage at 17 September 2020; (d) tillering stage at 2 Ocotber 2020; (e) ear emergence from boot at 19 October 2020; (f) anthesis stage at 2 November 2020; (g) dough development at 30 November 2020, first appearance of senescence; (h) ripening stage at 16 December 2020, complete senescence. "> Figure 4
Illustration of the hybrid retrieval workflow using the coupled PROSAIL-PRO models to establish a training database for the GPR, partly adapted from [<a href="#B82-remotesensing-14-04531" class="html-bibr">82</a>]. The output maps show our vegetation traits modeling over the BVCR area in Argentina. "> Figure 5
Goodness-of-fit results (RMSE, R<math display="inline"><semantics> <msup> <mrow/> <mn>2</mn> </msup> </semantics></math>) using AL (EBD) against validation data. (a) LAI model; (b) CCC model; (c) VWC model. "> Figure 6
Measuredvs. estimated wheat traits along 1:1-line including uncertainty intervals, using the EBD-optimised training data set. (a) LAI model estimates; (b) CCC model estimates; (c) VWC model estimates. "> Figure 7
LAI [m<math display="inline"><semantics> <msup> <mrow/> <mn>2</mn> </msup> </semantics></math> m<math display="inline"><semantics> <msup> <mrow/> <mrow> <mo>−</mo> <mn>2</mn> </mrow> </msup> </semantics></math>] of different crops in the study site, retrieved by the GPR model using simulated vegetation reflectance spectra, in-situ traits measurements and S2 surface multi-spectral reflectance data. Data obtained from BVCR wheat campaign 2020. (a) wheat tillering stage at 29 August 2020; (b) wheat booting stage at 28 September 2020; (c) wheat anthesis-flowering stage at 2 November 2020; (d) wheat dough development stage at 27 November 2020; (e) wheat ripening stage at 7 December 2020; (f) harvested wheat at 21 January 2021. "> Figure 8
CCC [g m<math display="inline"><semantics> <msup> <mrow/> <mrow> <mo>−</mo> <mn>2</mn> </mrow> </msup> </semantics></math>] of different crops in the study site, retrieved by the GPR model using simulated vegetation reflectance spectra, in-situ biophysical/biochemical measurements and S2 surface multi-spectral reflectance data. Data obtained from BVCR wheat campaign 2020. (a) wheat tillering stage at 29 August 2020; (b) wheat booting stage at 28 September 2020; (c) wheat anthesis-flowering stage at 2 November 2020; (d) wheat dough development stage at 27 November 2020; (e) wheat ripening stage at 7 December 2020; (f) harvested wheat at 21 January 2021. "> Figure 9
VWC [g m<math display="inline"><semantics> <msup> <mrow/> <mrow> <mo>−</mo> <mn>2</mn> </mrow> </msup> </semantics></math>] of different crops in the study site, retrieved by the GPR model using simulated vegetation reflectance spectra, in-situ biophysical/biochemical measurements and S2 surface multi-spectral reflectance data. Data obtained from BVCR wheat campaign 2020. (a) wheat tillering stage at 29 August 2020; (b) wheat booting stage at 28 September 2020; (c) wheat anthesis-flowering stage at 2 November 2020; (d) wheat dough development stage at 27 November 2020; (e) wheat ripening stage at 7 December 2020; (f) harvested wheat at 21 January 2021. "> Figure 10
Seasonal evolution of wheat cropland over the three paddocks at the BVCR study sites described by LAI, CCC, and VWC mean values of nice ESUs within the crop limits and the associated uncertainty, plotted as vertical bars. (a) LAI estimates; (b) LAI uncertainty (SD); (c) CCC estimates; (d) CCC uncertainty (SD); (e) VWC estimates; (f) VWC uncertainty (SD). "> Figure 11
Temporal evolution of wheat crop over the three paddocks at the BVCR study sites described by FVC mean measured values of nine ESUs within the crop limits and the associated SD, which is plotted as shadowed areas. "> Figure 12
Temporal evolution of wheat cropland over the three paddocks at the BVCR study sites described by LAI, CCC and VWC mean values of nine ESUs within the crop limits and the associated SD, which is plotted as vertical bars. (a) LAI vs. CCC temporal estimates; (b) LAI temporal estimates vs. FVC in-situ measured values; (c) LAI vs. VWC temporal estimates. "> Figure A1
Statistics(mean, standard deviation, min–max) of EBD-reduced final training dataset (blue) vs. validation dataset (red). Training data base was simulated with PROSAIL-PRO.). ">

Review Reports Versions Notes

Abstract

Earth observation offers an unprecedented opportunity to monitor intensively cultivated areas providing key support to assess fertilizer needs and crop water uptake. Routinely, vegetation traits mapping can help farmers to monitor plant development along the crop’s phenological cycle, which is particularly relevant for irrigated agricultural areas. The high spatial and temporal resolution of the Sentinel-2 (S2) multispectral instrument leverages the possibility to estimate leaf area index (LAI), canopy chlorophyll content (CCC), and vegetation water content (VWC) from space. Therefore, our study presents a hybrid retrieval workflow combining a physically-based strategy with a machine learning regression algorithm, i.e., Gaussian processes regression, and an active learning technique to estimate LAI, CCC and VWC of irrigated winter wheat. The established hybrid models of the three traits were validated against in-situ data of a wheat campaign in the Bonaerense valley, South of the Buenos Aires Province, Argentina, in the year 2020. We obtained good to highly accurate validation results with LAI: R

^{2}

= 0.92, RMSE = 0.43 m

^{2}

^{- 2}

, CCC: R

^{2}

= 0.80, RMSE = 0.27 g m

^{- 2}

and VWC: R

^{2}

= 0.75, RMSE = 416 g m

^{- 2}

. The retrieval models were also applied to a series of S2 images, producing time series along the seasonal cycle, which reflected the effects of fertilizer and irrigation on crop growth. The associated uncertainties along with the obtained maps underlined the robustness of the hybrid retrieval workflow. We conclude that processing S2 imagery with optimised hybrid models allows accurate space-based crop traits mapping over large irrigated areas and thus can support agricultural management decisions.

Keywords:

leaf area index; vegetation water and chlorophyll content; Gaussian processes regression; hybrid retrieval workflow; dimensionality reduction; active learning

1. Introduction

Wheat is a worldwide cultivated grain crop providing nearly 20% of all calories consumed due to its strong adaptability to various temperature and water conditions [1]. As an important global food crop, wheat yield information is essential in terms of food security [2]. Good agronomic practices and natural resource usage become fundamental for winter wheat grain development and crop yield. These practices have an important significance for modern precision agriculture [3]. Irrigated valleys, although representing only 20% of the world’s croplands, produce 40% of the global crop harvest [4]. In arid and semi-arid intensively cultivated areas, irrigation improves economic returns and can boost food production by up to 400% [4,5]. Hence, accurate and timely mapping of winter wheat biophysical and biochemical variable, or crop traits, results in being indispensable for enhancing crop management, food security, and agriculture structure adjustment [6,7,8]. The estimation and monitoring of crucial crop variables, such as leaf area index (LAI), canopy chlorophyll content (CCC), and vegetation water content (VWC), are key for agricultural applications including crop growth modelling and yield estimation supporting sustainable management of cultivated areas [6,9]. Definitions of these variables are given as follows.

LAI is defined as the total one-sided leaf area per soil unit (m

^{2}

leaf area per m

^{2}

soil). It is strongly related to canopy photosynthesis and evapotranspiration and plays a key role in the exchange of energy and water between the biosphere and atmosphere [10]. Varying definitions of LAI have been presented, such as effective plant area index [11], which includes the area from all plant organs and assumes random distributions of leaves. The green LAI (GAI) is probably the most relevant term describing the radiation transfer in vegetated canopies, i.e., the green photosynthetically active elements of the canopy [12,13]. Here, we refer mainly to GAI being most closely related to the signal actually recorded by satellite instruments. Nonetheless, for the sake of simplicity, we use the term LAI throughout this study. CCC is defined as the product of LAI and leaf chlorophyll content (C

_{a b}

) expressed in grams per unit leaf area (g/m

^{2}

). Hence, CCC represents the optical path in the vegetated canopy where absorption by chlorophyll pigments dominates the radiometric signal [14]. In addition, due to its strong relationship with nitrogen (N) content, often C

_{a b}

is considered a proxy for leaf N status. Therefore, CCC can be considered for quantifying canopy-level N content [15,16], although restricted to vegetative growth stages [17]. Due to its strong association with plant transpiration, vegetation stress, and biomass productivity, VWC has been widely considered a crucial variable for crop physiological status [18,19,20]. VWC of wheat crops is also a significant growth indicator during different development stages [6]. Water availability not only affects wheat photosynthesis, but also the grain filling rate and, ultimately, yield [21,22]. Traditionally, wheat crop water content has been determined by manually sampling plants to obtain fresh weight and dry weight and then calculating crop water content [2]. However, conventional field measurements are destructive and labour-intensive, especially across large areas with significant within-field spatial variability in soil infiltration, drainage, and conductivity characteristics [23]. The same holds true for destructive sampling of other crop variables. Earth observation (EO) technology can be used as an appealing alternative, providing a synoptic view of the Earth’s surface by making use of the complex interactions between radiation and the environment [24]. Nowadays, an attractive operational EO system with the capability to capture the visible (VIS), near-infrared (NIR), red-edge and short-wave infrared (SWIR) spectral domains involves the Sentinel-2 (S2) mission [25]. The recent S2 constellation that combines a high revisit frequency (5-day) and spatial resolution (10–20 m) with systematic global acquisition and an open access policy is promising in the development of operational farming services in near-real-time [16]. In recent years, the rapid development of EO technology has produced an unprecedented amount of research applying satellite data for monitoring crop traits across large areas, in a quick and accurate way e.g., [13,19,26,27,28,29,30].

With respect to the retrieval of crop traits from EO data, a multitude of techniques have been developed, which can be categorized into four main categories [30,31]: (1) Parametric regression approaches typically consist of relations between the traits and spectral data transformed into spectral indices e.g., [32,33]. (2) The second category of regression approaches obtains these variables from reflectance data using linear or nonlinear nonparametric approaches, i.e., chemometrics or machine learning regression algorithms (MLRA) e.g., [34,35,36]. (3) The category of physically-based retrieval methods refers to the inversion of radiative transfer models (RTMs) e.g., [37,38,39,40,41]. RTMs apply physical laws to explain the cause–effect relationships between radiation–photon interactions and plant constituents. The most prominent one-dimensional, (1D) turbid medium RTM for simulating vegetation reflectance is PROSAIL [42]. PROSAIL is composed of a leaf optical properties model, typically from the PROSPECT family [43,44], and, (b) the Scattering of Arbitrarily Inclined Leaves canopy architecture model (SAIL, 4SAIL) [45,46]. The last category of retrieval methods (4) is presented by hybrid approaches, being perhaps of most interest within current research lines and also for operational contexts [26,27,30,47,48]. Hybrid methods blend physics described by RTMs with the speed and efficiency of MLRAs. Within hybrid scenarios, training datasets are simulated by varying the input parameter space of coupled leaf-canopy RTMs. The nonlinear relationships between simulated reflectance and vegetation properties are then learned by an MLRA. Ideally, these training databases represent the canopy states realistically and at the same time are small enough to avoid exhaustive processing times. Regarding these strategies, a diversity of MLRAs have been successfully applied for retrieval tasks during the last few decades, e.g., artificial neural networks, decision trees or kernel-based methods [30,31]. As part of the latter family, especially Gaussian processes, regression (GPR) [49] emerged as a competitive retrieval algorithm e.g., [50,51,52], and has been widely adopted in studies inferring traits from EO data e.g., [48,53,54,55,56,57,58,59,60,61,62,63,64].

GPR received special attention also thanks to associated uncertainty estimates provided along with the predictions. This special feature enables assessing the fidelity of the models when transferring them into other spaces and times [28]. Moreover, to enhance mapping performance and processing speed, active learning (AL) techniques can be incorporated into hybrid workflows. AL aims to optimize training datasets through intelligent sampling using an iterative procedure devoted to exploiting samples during the design of the regression model [65,66,67].

Time-series of multiple quantitative traits over irrigated winter wheat present a captivating tool to address crop management questions and improve agricultural practises. Multiple authors have developed satellite-based retrieval models for wheat traits mapping using MLRA or parametric regressions [3,13,20,59,68,69]. Hybrid retrieval workflows may outperform these simpler empirical approaches due to their appealing property of combining physical awareness (RTMs) with flexibility of machine learning algorithms. Hence, they are an appealing alternative to be explored in our context. Several experimental studies have already demonstrated the feasibility of producing accurate vegetation traits maps from multispectral data, including S2 [59,70,71,72,73]. However, there is still the need to further test and adapt hybrid methods under diverse environmental conditions, at different locations and to multiple crop types.

In an attempt to demonstrate the possibilities of space-based cropland trait monitoring, this study presents the development of an optimised hybrid retrieval workflow dedicated to irrigated winter wheat monitoring in Argentina. Given the above-sketched general framework, this study aims to reach the following objectives: (1) to develop independent hybrid models optimised with AL and GPR for an explicit quantification of winter wheat LAI, CCC, and VWC from S2 data; (2) to generate accurate S2-based maps of the wheat traits with the inclusion of associated uncertainties over an intensive irrigated agroecosystem; and (3) to evaluate LAI, CCC, and VWC time series identifying seasonal trends over the selected study site.

2. Materials and Methods

2.1. Generation of Training Data Sets

Aiming to develop optimised hybrid models, the workflow initiates with simulating a training database using the coupled PROSPECT-PRO [44] and 4SAIL models, further referred to as PROSAIL-PRO. The models simulate bi-directional reflectance as a function of diverse leaf biochemical input variables, e.g., C

_{a b}

, leaf protein content (C

_{p}

), leaf carotenoid content (C

_{x c}

) or leaf equivalent water thickness (EWT) and biophysical variables, i.e., LAI and average leaf inclination angle (ALIA). The synthetic training database was generated by using Latin Hypercube Sampling (LHS) to randomly select samples from the multi-spectral input space. For this, 1000 combinations were drawn from all PROSAIL-PRO parameters. The selected size of the training dataset may appear rather small compared to classical retrieval approaches, for instance using lookup-table approaches or neural networks. Note that a standard implementation of GPR struggles to cope with thousands of samples within reasonable time intervals. Hence, as suggested by prior studies, e.g., [67,74], successfully training GP models with data sets of rather small sizes, we also restrict here to 1000 samples.

An overview of all included PROSAIL-PRO parameters including sampling, mean values, distributions and ranges is given in Table 1 (see also the study by Berger et al. [26] for full information about the generation of the training database). VWC can then be directly calculated from EWT with a (unit) conversion factor of 10’000 upscaling from leaf to canopy level through LAI. Fractional vegetation cover (FVC) was used to introduce the crop-bare-soil relationship (see Equation (1)). CCC was also obtained through extrapolating from C

_{a b}

to canopy level by means of LAI (see Equation (2)). Finally, simulated VWC and CCC both in [g m

^{- 2}

] were added to the training database:

V W C_{s} = (E W T \times L A I \times F V C) \times 10000 [{g / m}^{2}]

(1)

{C C C}_{s} = (C_{a b} \times L A I) / 100 [{g / m}^{2}]

(2)

Finally, bi-directional canopy reflectance was calculated using S2 spectral configuration and excluding low-resolution bands of 60 m since their focus is on cloud screening and atmospheric corrections. For instance, Estévez et al. [59] analyzed the contribution of the S2 spectral bands covering 10–20 m pixel resolutions for LAI retrieval purposes. The authors concluded to keep all ten bands (with S2 central wavelengths of 493 nm, 560 nm, 665 nm, 704 nm, 740 nm, 783 nm, 833 nm, 865 nm, 1610 nm, and 2190 nm) for further processing. To assure a maximum of spectral information required for the retrieval of multiple crop traits, we also decided to explore the ten S2 bands.

2.2. Gaussian Processes Regression

GPR is used as core MLRA in the hybrid model development. GPR modeling is flexible enough to fit many types of data, including geospatial and time-series data. GP is a Gaussian distribution over functions, which means that, instead of inferring the distribution of a parameter, GP can be then used to infer a distribution over functions directly under the premise that the function values are themselves random variables. At the inference stage, every time a new observation is made, the model hypothesis (prior probability distribution) is updated in light of the new observations. After having observed some function values, prior distribution over functions can be converted into posterior over functions. A Gaussian Process is a collection of random variables any finite number of which have (consistent) joint Gaussian distributions [49].

Notationally, GP can be deployed as follows: we observe training dataset

D = \{(x_{i}, f_{i} (x_{i})) | i = 1, \dots, N\}

, where

x_{i}

∈

R^{B}

(being

B

the number of sensor spectral bands) and

X = {\{x_{i}\}}_{i = 1}^{N}

, GP aims to provide a suitable model to predict the function outputs

f_{*}

given a test dataset

X_{*}

of size

N_{*} \times D

unknown by the model. In a noise-free environment, the joint distribution under GP is:

(\begin{matrix} f \\ f_{*} \end{matrix}) \sim N ([\begin{matrix} μ \\ μ_{*} \end{matrix}], [\begin{matrix} k (X, X) & k (X, X_{*}) \\ k (X_{*}, X) & k (X_{*}, X_{*}) \end{matrix}])

(3)

where

μ

and

μ_{*}

are the means of the functions f and

f_{*}

, which can be assumed as 0,

K = k (X, X)

is the matrix of self-similarities in the training dataset,

K_{*} = k (X, X_{*})

and

K_{*}^{T} = k (X_{*}, X)

represents the similarities between training and test datasets, and

K_{* *} = k (X_{*}, X_{*})

express the self-similarities in the test dataset. The mean

μ_{*}

(expectation of

f_{*}

) and the co-variance

Σ_{*}

of the test dataset

X_{*}

can be calculated as follows:

μ_{*} = E (f_{*}) = K_{*}^{T} K^{- 1} f

(4)

Σ_{*} = K_{* *} - K_{*}^{T} K^{- 1} K_{*}

(5)

To compute the measures of similarity between two points

x_{i}

and

x_{j}

, a positive definite kernel (or covariance function) that describes the covariance of the GP random variables should be implemented. The squared exponential kernel, also known as a radial basis function (RBF), which arises from taking the exponent of the scaled squared Euclidean distance between the data locations, is one of the most popular kernels used in GP modeling. It can be computed as:

k (x_{i}, x_{j}) = σ_{f}^{2} e^{- \frac{1}{2 ℓ^{2}} {∥x_{i} - x_{j}∥}^{2}}

(6)

where

σ_{f}^{2}

is the overall variance (

σ_{f}^{2}

that is also known as amplitude), and the parameter ℓ is the variance of the Gaussians (length scale) that controls the smoothness and confidence of the regression process. Thus, a Gaussian process is a distribution over functions whose smoothness is defined by the kernel k. RBF is a linear combination of basis functions; wherever there is data, we place a basis function scaled by the term

σ_{f}^{2}

. Consequently,

f_{*}

is the sum of basis functions scaled by

σ_{f}^{2}

. There are several manners to find the optimal ℓ value, cross-validation, maximum likelihood estimation (MLE), and Bayesian Learning are examples of different methods that can be implemented in the stage of learning the kernel parameters.

The key concept regarding GP is the fact that given two points

x_{i}

and

x_{j}

that are similar (similarity determined by the kernel k), the function values

f (x_{i})

and

f (x_{j})

, should be expected to be similar too. Straight away, the kernel function parameters have been optimised using a training dataset, and the kernel matrix can be plugged into the predictive equations for the mean

μ_{*}

and co-variance

Σ_{*}

of the test data

x_{*}

to obtain predictions

f_{*}

on the whole test dataset

X_{*}

. The process of making predictions is only computing the posterior distribution

p (f_{*} | X_{*}, X, f)

x_{*}

given f and the previous x, we can condition on the observations, to obtain:

f_{*} | X_{*}, X, f \sim N (f_{*}; k (X, X_{*}) k {(X, X)}^{- 1} f, k (X_{*}, X_{*}) - k (X, X_{*}) k {(X, X)}^{- 1} k (X_{*}, X))

(7)

Conditioning a Gaussian gives another Gaussian:

p (f_{*} | X_{*}, X, f) \sim N (f_{*} | μ_{*}, \sum_{*})

(8)

where

Σ_{*}

models the co-variances and cross-covariances between all combinations of train and test data. Predictions made using GP are not just point predictions: they are whole probability distributions.

More commonly, we have access to noisy observations

{\{x_{i}, y_{i}\}}_{i = 1}^{N}

with

y_{i} = f_{i} (x_{i}) + E

, where

E

is the noise of the observations. Assuming that the noise is Gaussian distributed with zero mean and variance

σ_{n}^{2}

, the noise can be therefore expressed as

E \sim N (0, σ_{n}^{2})

. The value of

σ_{n}^{2}

can be estimated from data utilizing the MLE principle. Assuming independently added Gaussian noise to each observation, the joint distribution of observations and test predictions is given by:

(\begin{matrix} y \\ f_{*} \end{matrix}) \sim N ([0], [\begin{matrix} k (X, X) + σ_{n}^{2} I & k (X, X_{*}) \\ k (X_{*}, X) & k (X_{*}, X_{*}) \end{matrix}])

(9)

To make predictions from a noisy GP, the only thing to do is to establish a conditioning on the observation:

f_{*} | X_{*}, X, y \sim N (k (X, X_{*}) {[k (X, X) + σ_{n}^{2} I]}^{- 1} y, k (X_{*}, X_{*}) - k (X, X_{*}) {[k (X, X) + σ_{n}^{2} I]}^{- 1} k (X_{*}, X))

(10)

The mean is linear in two ways:

μ_{*} = k (X_{*}, X) {[k (X, X) + σ_{n}^{2}]}^{- 1} y

(11)

The predictive co-variance (confidence intervals) is the difference between two terms:

Σ_{*} = k (X_{*}, X_{*}) - k (X_{*}, X) {[k (X, X) + σ_{n}^{2} I]}^{- 1} k (X, X_{*})

(12)

The first term is the prior variance, from which we subtract a (positive) term, telling how much the data X has explained. It is worth mentioning that the variance is independent of the observed outputs y. Although the posterior distribution

p (f_{*} | X_{*}, X, y)

covers noise in training data, it is still a distribution over noise-free predictions

f_{*}

. To additionally include noise

ε

into predictions

y_{*}

, we must add

σ_{n}^{2}

to the diagonal of

Σ_{*}

p (y_{*} | X_{*}, X, y) \sim N (y_{*} | μ_{*}, \sum_{*} + σ_{n}^{2} I)

(13)

Summarizing, as a Bayesian approach, GPR provides a natural and automatic mechanism to construct and calibrate uncertainties.

2.3. Active Learning Principles

GPR as a kernel-based regression method permits an enlarged number of input variables to be scanned and can be particularly convenient to deal with regression uncertainties since the confidence interval for the predictions is provided. Even though GPR offer a suitable alternative to processing large training datasets, processing large simulated training datasets could become computationally inefficient. Consequently, a reduction in the sampling domain, which specifies the size of the training dataset, is needed. A solution to the sampling reduction issue is given by semi-supervised approaches, and these techniques are also known as active learning [66]. AL aims to provide a strategy to efficiently reduce massive training datasets into an optimal dataset by sampling and evaluating the data pool iteratively through an intelligent process [67]. In AL, new informative samples are labeled based on the knowledge acquired during the exploration of the data input space.

Solving regression problems with AL can be addressed through a diversity framework. Diversity-based AL strategies create a reduced training dataset by selecting new samples from the input space, according to their dissimilarities [75]. Euclidean distance-based diversity (EBD) selects those samples out of the pool being distant from the already included ones in the training set, using squared Euclidean distance:

d_{E} = {∥x_{u} - x_{l}∥}_{2}^{2}

(14)

where

x_{u}

is a sample from the candidate set, and

x_{l}

is a sample from the training dataset. All distances between samples are computed and then the farthest are selected.

With the purpose to decrease the GPR training time, the EBD sampling method was performed in order to reduce the size of the input data space. Accordingly, 1% of the 1000 labeled samples (pairs of simulated bi-directional reflectance and vegetation traits) was randomly selected as the initial data set. The GPR training process was iteratively repeated up to 1000 times. At each iteration, the sample from the input data set with the largest distance between the pool of the 1000 labeled samples was selected by the EBD and added to the training data only when performance improved as evaluated by the root mean square error (RMSE) against the in-situ data [66,72].

2.4. Study Site

The Bonaerense Valley of Colorado River (BVCR) is located in the South of the Buenos Aires Province, Argentina, between the 39

^{\circ}

and 40

^{\circ}

S parallels and the 62

^{\circ}

and 63

^{\circ}

W meridians. Two different cropping systems can be well distinguished in the region: irrigated and non-irrigated crops. From the productive point of view, irrigation presents the highest productivity and contribution to the regional socio-economic movement. The area has a surface of 500,000 ha, of which 140,000 ha are irrigated by an extensive irrigation network based on uncoated drainage channels. Gravity irrigation has made the most of the agricultural activities in the area possible. The implanted crops under irrigation conditions occupy 91163 ha in the BVCR, including horticulture, pastures, and cereals. This study focused on three wheat paddocks of the BVCR in the Villarino district (see Figure 1).

2.5. Wheat Crop Experimental Design

A total of three repetitions were set up for each of the wheat paddocks. The position of nine elementary sampling units (ESU) was measured with a differential global positioning system (GPS) and recorded as point sample references in the field campaign database. The ESUs were revisited at each date corresponding to cloud-free S2 acquisitions to guarantee solid time-series. If there was a coincidence in ±6 days related to the S2 acquisition, the traits of the wheat crop were sampled.

The measurements were performed following a defined field protocol for the BVCR 2020 wheat campaign; see Figure 2. Hereby, the sampling strategy consisted of defining ESUs of 10 m × 10 m size for each wheat paddock in order to comply with the S2 spatial resolution (10 m). In addition, a minimum distance of 30 m from the parcels’ edges was kept. Sampling points were distributed within the ESUs following a five-measurements square spatial sampling strategy providing a statistically averaged LAI estimate per ESU [69]. The centre of the ESUs (sampling point A) was georeferenced using the S2 pixel grid to assure correspondence to the S2 reflectance data.

The in-situ measurement protocol used to carry out the field campaigns involved the collection of four main traits: LAI, fractional vegetation cover (FVC), C

_{a b}

, and above-ground fresh biomass (AGFB). The crop phenological stage and field data observations were also registered in the in-situ measurements database. The phenology was determined according to secondary growth stages of the Zadoks-scale [76].

Wheat Crop Management and Field Data Collection

The wheat crop was sown on 25 June 2020 at three different paddocks in the study site, a uniform seed density was established at 95 kg ha

^{- 1}

, and the weight of 1000 seeds was 45.8 g on average. At the wheat sowed stage, 80 kg ha

^{- 1}

of diammonium phosphate was used as a fertilizer considering the formulation 18-46-00 (18% N, 46% P

_{2}

_{5}

, and K

_{2}

O 0%). Phosphorus fertilization is a key factor in replenishing soil nutrients and obtaining more vigorous plants. Moreover, a rapid formation and growth of the root system are promoted making plants more resistant to water deficit. A system of ground edges separated by 14 m was implemented at the pre-sowing labour stage to ensure water embankments during the irrigation processes. Sunflower hybrid seed was the predecessor crop. At the beginning of the tillering stage, a value of 243 wheat plants per square meter (plant density) was measured. The distance between rows (rows spacing) was set as 17.5 cm according to wheat crop management recommendations in the region of the Colorado river valley. Wheat fertilization took place in two instances throughout the phenological cycle, first on 31 August 2020 and second by the middle of September. The dose of nitrogen (Urea) supplied was distributed uniformly and ranged between 150–200 kg ha

^{- 1}

. It is reasonable to expect that wheat plants start to receive nutrients from Urea up to one week later of fertilization. Three gravity irrigations were performed at different instances of the crop cycle between late August and November.

A collection of field data was obtained from the study area of the BVCR. In the year 2020, the test site was visited regularly during the growing periods from August to December. Data were collected on wheat fields belonging to the communal farmlands of BVCR by a group of expert professionals and technicians of HAEE-INTA, Argentina. The exact sampling sites were the same across the wheat 2020 field campaign. Each site was confined to a 10 × 10 m area, corresponding to the average S2 pixel size.

LAI measurements were taken using the PocketLAI R Smart-App [77]. Six observations were made-up per ESU and averaged giving a total of 54 LAI values for the entire campaign.

The FVC was measured utilizing the “Canopeo^®” App [78]. The app uses the RGB camera of the smart devices applying the relations R/G, B/G, and 2G-R-B to determine the crop canopy coverage percentage related to the soil.

To determine the AGFB value, five wheat plants per ESU were cut at soil level at each sampling date throughout the growing season. A sowing track or an area of 0.02 m

^{2}

was considered to obtain the AGFB samples. Each fresh biomass sample made up of wheat leaf, stalk, and fruit was enclosed in a sealed plastic bag and brought to the laboratory inside of a cooler for the fresh weight (FW) and dry weight (DW) measurements (see Table 2). The AGFB samples were weighed in the fresh state and oven-dried until reached constant weight for 24 h at 60 °C before dry weight was determined [79]. The DW corresponds to the above-ground dry biomass (AGDB) vegetation biophysical variable expressed in g units.

The SPAD-502 instrument (Minolta r) was used to perform C

_{a b}

measurements. Six leaves per wheat plant were selected randomly at each ESU, and five SPAD values were registered per wheat leaf and then averaged giving a total of 54 C

_{a b}

values. To extrapolate from SPAD non-dimensional values to C

_{a b}

values expressed in µg cm

^{- 2}

, the following Equation (15) was implemented:

C_{a b} [\frac{μ g}{{cm}^{2}}] = 12.23 e^{0.0279 . S P A D}

(15)

LAI was used to upscale from C

_{a b}

to CCC. Factor 100 led to the correct ground surface unit by upscaling from leaf [µg cm

^{- 2}

] to canopy level [g m

^{- 2}

C C C [\frac{g}{m^{2}}] = C_{a b} [\frac{μ g}{{cm}^{2}}] \times L A I \times \frac{1}{100}

(16)

Plant leaf water content is commonly expressed as equivalent water thickness corresponding to a hypothetical thickness of a single layer of water averaged over the whole leaf area [80]. EWT relates the leaf’s water content to the leaf’s area (Al), so it is usually expressed in µg cm

^{- 2}

as follows:

E W T = \frac{(F W - D W)}{A l} [μ g {cm}^{- 2}]

(17)

Canopy water content (CWC) is commonly derived through extrapolation by means of the LAI:

C W C = E W T \times L A I [g m^{- 2}]

(18)

Due to the linkage of LAI to the area of leaves, CWC may neglect the water content contained in other organs, such as stalks and fruits. To monitor the total amounts of water stored in wheat plants, leaves, stalks, and fruits were included in our analysis. The FW and DW values determined in the laboratory per each ESU of the study site (consisting of wheat fresh and dry organic matter) were used as a proxy for the determination of VWC expressed in g m

^{- 2}

units. The VWC values were then obtained by calculating the difference between the FW [g] and DW [g] and referring it to the sowing area A, implicated in the field data collection process:

{V W C}_{h} = \frac{(F W - D W)}{A} [g m^{- 2}]

(19)

Equation (19) assumes a wheat plant homogeneously distributed in one square meter of soil, hence we refer here to VWC

_{h}

, with the subscript h standing for homogeneous. However, this approach can lead to errors in the estimation of total VWC values. During the first growing stages, plants are separated by rows leading to a proportion of the soil visible to the sensor. Since the S2 sensor acquires surface reflectance in a 10 m × 10 m ground sampling distance (GSD), the spectral value of each band summarizes both the soil and vegetation contributions. Only when the wheat crops reach the maximum greenness stage is the soil completely covered by vegetation. Hence, to be more precise in the determination of the VWC, we introduced a modification considering in-situ measured FVC in the calculation of plant water content per square meter of sowed wheat:

{V W C = V W C}_{h} \times F V C = \frac{(F W - D W)}{A} \times F V C [g m^{- 2}]

(20)

The calculated wheat CCC and VWC values were recorded as new inputs of the field campaign database as shown in Table 3.

The wheat phenological stage was registered along the BVCR campaign from August to December 2020 (see Figure 3 and Table 4).

2.6. Sentinel-2 Image Acquisitions

A series of S2 images consisting of 15 cloud-free L1C products from 29 August to 11 January 2020 was downloaded from ESA’s web server https://scihub.copernicus.eu/ (accessed on 20 June 2022).The S2 images were atmospherically corrected using the Sen2cor plugin [81] to obtain bottom-of-atmosphere (BOA) reflectance data from top-of-atmosphere (TOA) L1C products. As a result, 15 level 2A (L2A) S2 images were obtained, coinciding with ±6 days of in-situ sampling. The images were resampled to 10 m GSD and cropped according to the study site in the SNAP 7.0 software https://step.esa.int/main/snap-7-0-released/ (accessed on 18 July 2022).Bands B1 (443 nm), B9 (940 nm), and B10 (1375 nm) were excluded from the resulting products due to their low GSD of 60 m, being in correspondence to the simulations. Finally, the S2 products composed of ten spectral bands were used for further analysis within the AL optimization process (see Table 5).

2.7. Delineation of the Hybrid Retrieval Workflow

An overview of the complete hybrid retrieval workflow is given in Figure 4. Three well-distinguished conceptual blocks are detailed, starting with an RTM section, followed by in-situ data collection, and the retrieval of biochemical and biophysical traits. In summary, the implemented retrieval workflow consisted of the following four main steps:

1.: generation of the training database, i.e., simulated TOC reflectance with corresponding traits using the PROSAIL-PRO model;
2.: building the in-situ database containing multitemporal field measurements from the BVCR site and S2 spectra;
3.: optimizing the training database with AL-EBD and GPR, applying retrieval models to obtain wheat LAI, CCC, and VWC; and
4.: seasonal mapping of the three crop traits over irrigated wheat fields and corresponding uncertainties using S2 scenes.

Figure 4. Illustration of the hybrid retrieval workflow using the coupled PROSAIL-PRO models to establish a training database for the GPR, partly adapted from [82]. The output maps show our vegetation traits modeling over the BVCR area in Argentina.

Note that, according to the 2020 wheat campaign observations (see Table 4), the senescence of the crops started on 16 November 2020. In order to prioritize the retrieval models’ response within the vegetation greenness period, the in-situ database was restricted from 3 September to 16 November 2020. Additionally, 12 non-vegetated and 27 senescent crop samples were added to the database aiming to improve the robustness of the crop traits retrieval models after starting yellowing and senescing crop growth stages (see Appendix A for more information on the training database). For model evaluation, the coefficient of determination (R

^{2}

), the RMSE, the normalized root mean square error (NRMSE), the mean absolute error (MAE), and the mean absolute percentage error (MAPE) were registered for each crop trait model.

The hybrid workflow was entirely built within the Automated Radiative Transfer Models Operator (ARTMO) toolbox [83]. ARTMO was developed as a modular graphical user interface in Matlab, to automate the simulation of RTMs [84]. The toolbox brings multiple RTMs together with essential tools required for the retrieval of a diversity of biophysical and biochemical vegetation traits. ARTMO has been expanded over the years with all kinds of RTMs and image processing options, such as the MLRA toolbox [85] with included AL module [74], emulation, sensitivity analysis, and scene generation. More information can be found at: http://artmotoolbox.com/ (accessed on 18 July 2022).

3. Results

3.1. Optimized Sample Selection for LAI, CCC and VWC Modeling

In order to build efficient GPR-based retrieval models for wheat traits, at first, the AL EBD technique was explored for optimizing the training samples. In Figure 5, we demonstrate the results of training database reduction through EBD as a function of RMSE (left y-axis) and R

^{2}

(right y-axis). For all three variables, an initial dataset of 10 samples was used as starting point for the AL procedure.

Figure 5a shows the effect of the training data size on the LAI models’ accuracy, with optimal results obtained after adding 112 samples (RMSE: 0.42 m

^{2}

^{- 2}

; R

^{2}

: 0.92). For CCC, the highest accuracy was obtained after adding 137 samples reaching RMSE: 0.27 g m

^{- 2}

and R

^{2}

: 0.8 (see Figure 5b). Figure 5c, shows that the optimal accuracy of the VWC model was achieved after adding 232 samples (RMSE: 416 g m

^{- 2}

; R

^{2}

: 0.76). Note hereby that lowering the RMSE does not necessarily go along with an improvement of R

^{2}

, as can be read on the right y-axis of Figure 5. Although it follows the same general trends as RMSE, the pattern provided by R

^{2}

is more irregular, indicating it as a less reliable measure than RMSE for AL testing [72].

Figure 6 presents the scatter plots of estimated against in-situ measured LAI, CCC and VWC samples. Overall, values for R

^{2}

and RMSE indicate relatively high agreements between modelled and measured wheat traits. In addition, uncertainties are provided as standard deviation (SD) expressed by the colour table. The SD of LAI, ranging from 0.62 to 1.23 m

^{2}

^{- 2}

, provides a proxy of the model’s uncertainty and is indicated by the colour bar close to the regression graph (see Figure 6a). Figure 6b shows the resulting optimised CCC model, with uncertainty interval ranging from 0.2 to 0.48 g m

^{- 2}

. In Figure 6c, the optimised VWC model are displayed, with uncertainty interval ranging from 370 to 542 g m

^{- 2}

3.2. Lai, CCC, and VWC Mapping

Next, the optimally trained hybrid GPR models were applied to the six selected dates of the 2020 BVCR wheat campaign. The scenes cover a variety of crop types, including the wheat paddocks selected for this study. Non-vegetated land covers, comprising bare soils, water bodies, and man-made surfaces, were removed from the scenes by applying a vectorial mask.

Figure 7 displays the six S2-derived LAI maps over the BVCR site during the growing season of 2020 (ending 01/2021). All grey-coloured areas correspond to fallow or harvested fields or dried-out natural vegetation. In contrast, LAI distribution of the intensively cultivated crops may vary during the whole growing cycles of the BVCR 2020 campaign. The wheat crop was sown in late June 2020, leading to the emergence of the main wheat stem and three tillers in late August. The highest LAI values during the growing season are identified: LAI reaches a maximum greenness stage at the beginning of November 2020 with 4.12 m

^{2}

^{- 2}

. Wheat senescence starts in the middle of November and harvest takes place by the first days of January 2021. Figure 8 shows the multitemporal maps resulting from the CCC estimates for the BVCR study site. The CCC values oscillate between 0 and 2 g m

^{- 2}

being the maximum detected on 2 November. Later in the growing season, from around 16 November 2020 onwards, all wheat parcels exhibited a markedly low CCC driving toward senescence. Due to the linkage between LAI and CCC, both variables are strongly correlated across the whole wheat’s growing cycle.

VWC maps are displayed in Figure 9 for the six selected dates of the 2020 BVCR wheat campaign. These maps show spatio-temporal changes in VWC for all crop parcels at the BVCR study site, which were in parity with expected changes in VWC, as driven by rainfall conditions and irrigation management operations during the whole wheat’s phenological cycle. Since the VWC was calculated with FVC (see Equation (20)), VWC strongly depends on the percentage of bare soil present in the scene’s pixel observed by S2. FVC values widely vary during the vegetation growing period determined primarily by the crop development and structure as well as the implemented crop management practices. Typical VWC values for wheat crop vary from close to 0 g m

^{- 2}

at the start and end of the season to around 2500 g m

^{- 2}

in late November 2020.

3.3. Wheat Phenology Based on Multi-Temporal LAI Maps

According to the Z1.3 state of the Zadoks’ scale and collected field data, three leaves of wheat plants emerged in August 2020. Hence, most of the wheat pixels in the S2-observed scene have LAI < 0.5 m

^{2}

^{- 2}

. This is clearly noticeable in the LAI map produced for the date 29 August 2020 (see Figure 7a and Table 4). Later in the growing season, around 28 September 2020 (Z4.3, wheat at the tillering stage), all crop parcels revealed distinct LAI growing curves due to field irrigation and fertilization which took place by the third week of September. The S2-derived LAI map for the date 28 September 2020 reveals LAI values of about 2.5 m

^{2}

^{- 2}

for the wheat crop (see Figure 7b).

Wheat plants were at the flowering stage on November 2 (Z6.1—few anthers at the middle of the ear) and the plants reached a height of around 80 cm. Hence, LAI ranged from 3.5 to 4.5 m

^{2}

^{- 2}

by this date and almost all the cropped pixels at each paddock have experienced the beginning of the anthesis, resulting in higher LAI values (see Figure 7c). Note that some distinct heterogeneities are visible principally in paddock 322; it had a variety of vegetation vigor. In the anthesis stage, the wheat ear arises from the plant stalk, implying that all wheat paddocks are in the flowering stage.

The milk development stage started once the flowering is completed during the month of November, and wheat plants also reached their maximum greenness on 16 November 2020 (Z7.5—Medium milk). Moreover, senescence starts at this time and young leaves may become yellowish. Accordingly, wheat parcels had LAI > 4 m

^{2}

^{- 2}

, which represents the LAI curve’s inflection point hence from this date onwards LAI decreases until wheat harvest. The wheat traits maps show the early kernel formation stage which occurs one to two weeks after pollination. Once the plants’ dehydration period has begun the dough development stage starts. On 30 November (Z8.7—Hard dough), most of the kernel dry weight (starch and protein content) accumulates and therefore LAI < 3.5 m

^{2}

^{- 2}

on average, this stage is also known as physiological maturity. The LAI map generated for the date 27 November 2020 displays the decrease of this biophysical variable over the BVCR study site (see Figure 7d) indicating the end of plant growth.

Finally, by the middle of December 2020, wheat parcels revealed complete senescence, and the seed moisture decreased down to 13 to 14%, which corresponds to the ripening stage. The wheat seeds could have been harvested at the end of this stage on 16 December 2020; however, wheat plants remained standing until the first days of January 2021. This can be appreciated in the LAI maps produced for the date 7 December 2020 where (green) LAI values oscillate between 1.5 and 2 m

^{2}

^{- 2}

(see Figure 7e). At the end of the ripening stage, the seeds were harvested, which resulted in pronounced decrease in LAI on 4 January 2021 (see Figure 7f).

3.4. Seasonal Analysis of Retrieved Traits

This section explores the seasonal evolution of LAI, CCC and VWC along the wheat 2020 campaign. The trait models were applied to a total of 15 free-cloud L2A S2 scenes to obtain sufficiently dense time series along the wheat seasonal cycle. We also explored the GPR model uncertainty during the course of the season. The uncertainties are mapped in the form of SD to describe the variation of LAI, CCC and VWC values. LAI, CCC and VWC temporal gaps correspond to S2 cloudy acquisitions. Figure 10a displays the temporal profile of retrieved LAI and Figure 10b the LAI SD profiles of the nine in-situ measurement points corresponding to the three paddocks of wheat cropland in the BVCR study site. The figure shows the temporal evolution of the averaged nine ESUs (solid line), and the mean of the in-situ measured LAI values (yellow dots). Figure 10c,d show the temporal profile of the retrieved CCC, the SD, and the mean of the in-situ measured values of all ESUs. Figure 10e,f illustrate the VWC and SD trends. Seasonal patterns of LAI, CCC and VWC reproduce the typical phenological stages of the wheat crop. LAI, CCC and VWC increase steadily, reaching a maximum peak in November 2020. By that time, a plant’s dehydration process also emerges, and senescence starts. During the senescence stage, plant water content drops significantly, leading to the wilting of leaves and thus causing LAI to decrease [86]. Abrupt drops of LAI, CCC and VWC in January 2021 suggest the harvest event.

3.4.1. In-Situ Measured FVC Time Series Analysis of Irrigated Winter Wheat Crops

In the course of wheat tillering, plants reach 10–20 cm height on 17 September and FVC values oscillate between 25–45% (see Figure 3c and Figure 11). The crop leaves remain horizontal during the tillering stage with four leaves per tiller leading to increase of vegetation coverage to FVC about 30–65% on October 2 (see Figure 3d). From the booting stage onwards, the wheat FVC decreases due to plants’ stem elongation up to a minimum of 22% on average on November 2. Figure 12b, showing the seasonal patterns of the retrieved LAI vs. the in-situ measured FVC, reflects this behavior: a strong dropping of FVC at the beginning of November (minimum on 2 November), while LAI values remain high in this period (around 3.7 m

^{2}

^{- 2}

) and only decrease at the end of November with starting senescence. This specific behavior of FVC can be explained by the fact that the Canopeo^® application uses only the RGB spectral bands to determine the crop canopy coverage.

Later on in the growing season, by November 16, the FVC values increase up to 86% (see Figure 11). During the senescence period, the FVC values remain high; wheat plants remain standing, and the FVC in-situ measured values range from 80 to 90.8% (see Table 2).

3.4.2. Seasonal Analysis of S2-Retrieved CCC and LAI

Figure 12a presents the temporal profiles of the mean estimated CCC and LAI values along the wheat’s seasonal cycle. A strong positive correlation between LAI and CCC temporal profiles can be seen as C

_{a b}

values have been upscaled with LAI to obtain CCC. CCC retrieved values along the seasonal cycle of winter wheat range from 0 to 1.7 g m

^{- 2}

approximately. At the start of the crop growing season and during the greenness phase, the CCC curve follows the LAI’s trend until 27 November. After the vegetation has reached the maximum greenness, which typically corresponds to the LAI peak value, the senescence phase begins until crops are ultimately harvested throughout January 2021. From 27 November onwards, during cropland senescence, the difference between LAI and CCC profiles starts to become more noticeable. On 7 December, 1.5 < LAI < 2 m

^{2}

^{- 2}

indicates the presence of vegetation structures, while 0.2 < CCC < 0.4 g m

^{- 2}

suggests that the chlorophyll content of wheat plants is declining until, finally, the crop turns yellowish due to the absence of the pigments. According to the in-situ collected notes (see Table 4), the complete senescence for the BVCR 2020 campaign of winter wheat crop at the study site took place by 16 December. Figure 12a shows how the CCC curve drops to zero from 22 December until 21 January. The agreement between both data sources confirms the capability of S2 to retrieve the CCC for monitoring winter wheat croplands.

3.4.3. Seasonal Analysis of S2-Retrieved VWC and LAI

Figure 12c shows the seasonal evolution of the mean estimated value of VWC against LAI. As the VWC model is correlated to FVC in-situ measured values instead of LAI ones, a notorious time shift can be appreciated between the two curves. While the LAI curve describes the leaves’ surface per square meter of soil, VWC involves not only the leaves but also the stem and spikes of wheat plants. In addition, third wheat irrigation occurred on 17 November to boost the grain development. Consequently, although senescence had started on 16 November, the VWC values remain high at around 2300 to 2400 g m

^{- 2}

until 27 November. Once the wheat plants have reached the complete senescence by the middle of December 2020, the mature stems remain with a certain level of humidity, thus VWC = 750 g m

^{- 2}

approximately on 22 December.

4. Discussion

With the ambition to develop retrieval models for mapping the diversity of Argentinian winter wheat traits, time-series of S2 imagery was processed to obtain LAI, CCC and VWC maps along the phenological crop cycle in an irrigated, intensively cultivated area. This research is built upon a hybrid retrieval workflow using an optimised simulated training dataset and in-situ measurements of vegetation traits from the BVCR 2020 campaign. In the following, advantages and limitations of the RTMs (Section 4.1), performance of GPR models (Section 4.2), potential of the temporal vegetation traits mapping for wheat agronomic management (Section 4.3), study limitations (Section 4.4), and finally operational monitoring opportunities (Section 4.5) are discussed.

4.1. Advantages and Limitations of Coupled RTMs

Physically-based models simulating the influence of water content and pigments on the obtained signal acquired by satellite sensors potentially lead to improved predictions of crop traits compared to empirical approaches. Regarding the leaf level, we chose the PROSPECT model (i.e., the actual PROSPECT-PRO version), being perhaps the most widely applied leaf RTM in EO vegetation studies. PROSPECT had been first introduced in the 1990s by Jacquemoud et al. [87] and was continuously updated, with for instance improved calibration of the specific absorption coefficients of C

_{a b}

(introduced in PROSPECT-D [88]) and separation of dry matter content into leaf protein content and carbon-based constituents in PROSPECT-PRO [44]. Although results may not change significantly, the usage of the most actual model version is always recommended. Still, improvements to the model may be required, such as the inclusion of specular reflection [88].

To upscale from leaf to canopy level, C

_{a b}

and EWT were multiplied by LAI simulated with the SAIL model to obtain CCC and VWC, respectively. The SAIL model assumes that the horizontal direction of the canopy is infinite [45]. Therefore, it is mainly suitable for estimating LAI and other traits of rather continuous and homogeneous crops [89]. Winter wheat is not a pronounced row crop (like maize) but belongs to the grass family (Poaceae). Hence, with its turbid medium assumptions, the choice of SAIL (or PROSAIL when combined with PROSPECT) is justified for our analysis. In addition, we used the simulated FVC data for driving calculations of the VWC training dataset (see Equation (1)). This approach differs from previous studies, e.g., [90,91] where canopy water content was obtained by upscaling EWT by LAI (see Equation (18)).

To optimise the simulated training dataset, we extensively explored the parameterization of winter wheat. Multiple settings were tested to achieve the final VWC model minimising the RMSE. The EWT values ranged from 0.0001 to 0.182 g cm

^{- 2}

, being the most convenient interval presented in Table 1. The in-situ collected variables, LAI, AGFB and AGDB, allowed us to calculate the EWT for winter wheat, involving plant leaves, stalks, and fruits. The calculated EWT values from the campaign data oscillated between 0.019 to 0.175 g cm

^{- 2}

, leading to an uncertainty increase in the final built model (RMSE = 587 g m

^{- 2}

). This can be argued by the fact that EWT implies only water content at the leaf level, whereas field measured data include the water content of all plant organs. Varying the EWT in the PROSAIL-PRO training dataset generation not only affected the VWC model statistics but also the R

^{2}

and RMSE of the LAI and CCC obtained models, which highlighted the sensitivity of PROSAIL-PRO to this leaf variable.

4.2. Performance of Hybrid GPR Models

The proposed workflow is built on the following fundamental factors. First, the careful design of the in-situ measurements methodology led to a collection of high-quality field wheat traits data. Second, the in-situ database covered multi-temporal vegetative stages. In this respect, the GPR models achieved high accuracy for the studied traits, with trustful estimates along the crop greenness phase. Third, we added non-vegetated spectra to the training dataset, i.e., from bare soils, man-made surfaces, and senescent crops as the retrieval models needed to be adapted to the spectral variability of actual land covers of the scenes. Note that, with these additional spectra in the learning process, a slight decrease in the validation may occur; however, more generic wheat traits models can be obtained [58,72].

Active learning provides a powerful and smart sampling method to decrease the redundancy of the training dataset maintaining the spectral variability and thus information content. Finally, from a total of 1000, 112 samples of LAI, 137 samples for CCC, and 232 samples for VWC were intelligently chosen to train the retrieval algorithms. In general terms, the three traits GPR models performed well in learning the nonlinear relationship between the reflectance of S2 spectral bands and LAI, CCC and VWC. This confirms the results of other hybrid studies [59,60,62], where these traits were successfully estimated from S2 data with a GPR model.

The provision of per-pixel associated uncertainty information helps to understand how trustful results are, required to appropriately use the information obtained through a model or measurement [92]. Such information can further indicate the portability of the methods in space and time and is important when the models are being tested across different crops or sites.

The seasonal analysis of LAI, CCC, and VWC revealed differences in the uncertainty for both green and senescent development phases of winter wheat. Generally, low uncertainties refer to spectra that were well represented during the learning phase, whereas retrievals with high uncertainties refer to spectral information that deviates from what has been represented during the training stage of the models [26]. During the growing season, the availability of in-situ measured traits renders the models more robust and confident. Additionally, the reasonably low uncertainty provided by the GPR retrieval models during this period for most paddocks at the study site suggests reliability in the mapping results. From 16 November onwards, the senescence started, and consequently, winter wheat spectra changed drastically. Additionally, field data were absent for this period; therefore, the uncertainty of the GPR models increased considerably (see Figure 10).

4.3. Potential of Seasonal Traits Mapping for Wheat Agronomic Management

Monitoring of winter wheat traits provides valuable information on the life cycle of the plants needed for efficient crop irrigation and fertilization management being critical factors to crop yield. Winter wheat in our study site was irrigated in the transition between phenological stages. Dates of fertilization were close to those of irrigation since the urea is incorporated into the soil with water, and fertilization occurs first followed by irrigation (see Figure 12). Wheat water absorption depends on the previous hydric stress condition as well as plants’ phenological stage. Effects on the vegetation water content can be noticed up to 4 or 5 days after crop irrigation [93,94,95]. Nitrogen fertilization increases the crop water use efficiency (WUE) due to the higher concentration of nitrogen in plants, leading to improved radiation use efficiency. C

_{a b}

is positively correlated to photosynthetic activity [96]; consequently, C

_{a b}

is an indirect indicator of plants’ total nitrogen content [97]. Plant stress conditions may affect the photosynthetic process generating nitrogen assimilation and carbon fixation reduction. The leaf chlorophyll content is a genetic attribute of wheat plants; consequently, it remains almost constant during most crop growth stages [3].

CCC monitoring is crucial to understanding plants’ photosynthetic activity. Repetitive LAI mapping of winter wheat can provide not only valuable information about plants’ leaf area but also a proxy for C

_{a b}

upscaling to canopy level. CCC monitoring provides useful information to infer plants’ salinity and water stress [98,99], as well for the early detection of nitrogen absorption deficiencies and decrease yield associated with plant stress conditions [14,100].

In Figure 12a, the winter wheat fertilization and irrigation dates are presented against the LAI and CCC trends throughout the crop season. The first and second irrigation take place during the tillering stage (Z2.3 and Z2.4 according to the Zadoks growth scale); as a consequence, the slope of the LAI and CCC curves increase. This analysis reveals the potential of the seasonal mapping of irrigated winter wheat traits in Argentina to manage the natural resource usage in arid and semi-arid regions.

The WUE of crops can only be roughly estimated, since, under field conditions, it is arduous to know precisely how much water has been consumed by plants and their growth in accumulated biomass [101]. Consequently, water consumption is usually estimated from indirect data on precipitation and the volume of water lost through runoff, percolation, or direct evaporation from the soil and which has never been consumed by plants [102,103].

Third irrigation occurs during the wheat’s milk development stage to boost grain yield (Z7.5 according to Table 4). Figure 12c shows the VWC curve response to this irrigation while LAI has started to decrease by the middle of November 2020.

Well-adjusted VWC models represent an easy-access tool for helping to obtain a response to the most fundamental questions: when and how much to irrigate winter wheat crops in arid and semi-arid regions. To support this analysis, one of the most tested and widespread tools is the hydric balance, which allows for estimating the water content in the soil and deciding the optimal moments of irrigation.

The data fusion of LAI, CCC, and VWC temporal profiles presented in our research constitutes a differentiating element in the field of remote sensing monitoring of irrigated winter wheat development at the study site. (see Section 3.4.2 and Section 3.4.3). The potential yield generated by irrigated wheat in the arid region of the focused study site can be compared with the grain production of the harvested winter wheat crops from the agricultural core area of Argentina. Consequently, monitoring crop traits results is crucial to reducing fertilizer and irrigation water uptake assuring high yields. All these statements point toward the potential of the vegetation traits retrieval models to address the agronomic analysis of irrigated winter wheat crops.

4.4. Study Limitations

Even though the in-situ collected data process was conscientiously carried out, the uncertainty of the measurements and used instruments were not considered in the analysis. However, uncertainty coming from field data has to be kept in mind when assessing the overall accuracy of the retrieval models.

From early August to mid-November, winter wheat plants remain green; thereby, the in-situ LAI samples correspond to the crop’s growing season. During this period green, LAI (i.e., GAI) increases as the wheat plants grow. Once crop senescence arrives, plants remain standing and the share of brown leaves increases [13,68]. Each of the winter wheat paddocks was revisited periodically from 10 August to 16 December 2020. The in-situ LAI values were taken from 3 September to the end of the BVCR 2020 wheat campaign; nevertheless, only the green LAI samples were considered to train the LAI model (from 3 September to 16 November). From 16 November to 16 December, in-situ LAI values continue to rise although the cropland has become senescent (see Table 2).

Since in-situ LAI values were used to tune the GPR retrieval models with the AL-EBD techniques, the final LAI retrieval model is better adapted toward green vegetation states. Hence, the uncertainty of the developed models will increase during later (mature/senescent) growing stages.

Finally, it must be noted that our study was conducted with data coming from one single growing season. Thus, temporal transferability to other seasons remains an open question. However, since GPR models provide associated uncertainty estimates, the fidelity of temporal and spatial transferability of the retrieval models to other regions with different cropping environments can be assessed to some degree [104].

4.5. Opportunities for Operational Monitoring

The possibility to estimate the irrigated winter wheat traits at 10 m or 20 m GSD would offer important opportunities for operational precision farming applications. The promising line of hybrid retrieval workflows opens some new branches to be explored in the field of agricultural applications using EO data [105]. The following questions arise: (1) Can we extrapolate the traits models for winter wheat crops to other regions and times? (2) Are the crop traits models sufficiently generic to be applied to other crop types in arid and semi-arid intensively cultivated areas? (3) Are the outcomes of this research useful to address the management of the crops more efficiently, for instance determining when to irrigate and to fertilize optimally?

In this context, the availability of the winter wheat in-situ database obtained from the BVCR 2021 campaign represents an excellent opportunity for the validation and refinement of the developed models and will be explored in a future study.

With respect to the implementation of the GPR models into operational processing, recent possibilities are opened thanks to cloud computing technologies. In a pioneering study, [106] presented an approach on how to integrate GPR models into Google Earth Engine (GEE), e.g., for S2 processing. In this way, the GPR models can be applied anywhere and anytime in the world. Pursuing this approach, first applications already emerged regarding multiple crop traits mapping and monitoring of key growing cycle events (e.g., start/end-of-season) [60] or for wide-scale S2 cropland traits mapping, e.g., over the entirety of Germany [62]. GPR models developed here can be likewise implemented into GEE for monitoring applications. For instance, such a straightforward monitoring interface can help users to understand the appropriate rates and timing of fertilization and irrigation.

5. Conclusions

We presented a hybrid retrieval workflow for operational mapping of LAI, CCC, and VWC that was optimised for an irrigated winter wheat cropland.

The implemented hybrid retrieval method used the advantages of the radiation-surface physics interaction being considered by the RTM and at the same time flexibility, scalability, and computational speed that are provided by GPR. The usage of GPR as a solid probabilistic regression algorithm provides the additional advantage of delivering associated uncertainties along with the trait estimates, helping to assess the fidelity of the mapping results.

The hybrid models for the winter wheat traits were validated with relative high accuracies against in-situ data (LAI: RMSE = 0.43 m

^{2}

^{- 2}

; CCC: RMSE = 0.27 g m

^{- 2}

; VWC: RMSE = 416 g m

^{- 2}

). They were subsequently applied to a series S2 imagery during the growing season of 2020 over the South of Buenos Aires Province, Argentina. All produced winter wheat maps showed consistent spatiotemporal performance and relatively low associated uncertainty intervals. The seasonal evolution of the winter wheat crops was quantified by correlating in-situ collected data and interpreting crop growth stages of the S2-derived maps. We conclude that routine mapping of wheat LAI, CCC and VWC constitutes a unique opportunity to monitor irrigated winter wheat development, which will contribute to managing sustainable agricultural production and thus assuring food security.

Author Contributions

G.C., J.D. and J.V. proposed the general objectives and goals of the research; A.P., C.W., P.S.A. and A.C. designed the field campaigns and collected the in-situ data; G.C., J.V. and J.P.R.-C. analyzed the data and obtained the results; G.C. and K.B. wrote the paper; K.B., J.V. and J.D. reviewed the paper and supervised all the procedures. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the European Research Council (ERC) under the ERC-2017-STG SENTIFLEX project (Grant No. 755617) (K.B.) and Ramón y Cajal Contract (Spanish Ministry of Science, Innovation and Universities) (J.V.).

Data Availability Statement

Not applicable.

Acknowledgments

We would especially like to thank the Hilario Ascasubi Agricultural Experimental Station of the National Institute of Agricultural Technology of Argentina for the experimental datasets used in this study. This publication is also the result of the project implementation: “Scientific support of climate change adaptation in agriculture and mitigation of soil degradation” (ITMS2014 + 313011W580) supported by the Integrated Infrastructure Operational Programme funded by the ERDF. The research was also supported by the Action CA17134 SENSECO (Optical synergies for spatiotemporal sensing of scalable ecophysiological traits) funded by COST (European Cooperation in Science and Technology, www.cost.eu (accessed on 22 July 2022)).

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. EBD-Reduced Final Training Dataset versus Validation Dataset

With the purpose of observing the spectral similarity between the final training data set (Composed of couples of wheat traits and vegetation spectral reflectance simulated by the RTM PROSAIL-PRO) and validation dataset (Composed of wheat traits associated with S2 surface spectral reflectance) from BVCR campaign, both data sets were statistically compared. Figure A1 illustrates averaged training and validation spectra, with standard deviation and ranges in partly transparent colours. It can be observed that the training data match closely with the validation data in the red-edge region and from the NIR on, the difference between them increases. In the visible region of the electromagnetic spectrum, a difference in both training and validation means curves can be appreciated, this phenomenon is attributed to the fact that the MSI on board of S2 satellite integrates vegetation and bare soil reflectance in a 10 × 10 m pixel. The proportionality between vegetation and bare soil for each pixel varies with the fractional vegetation cover. The broader range of the validation data is due to the senescent vegetation and bare soil added spectrums, which implies a sufficient degree of generalization to ensure that the model is applicable to S2 time series data and able to retrieve the vegetation traits across the complete phenological cycle.

Figure A1. Statistics(mean, standard deviation, min–max) of EBD-reduced final training dataset (blue) vs. validation dataset (red). Training data base was simulated with PROSAIL-PRO.).

References

Shiferaw, B.; Smale, M.; Braun, H.J.; Duveiller, E.; Reynolds, M.; Muricho, G. Crops that feed the world 10. Past successes and future challenges to the role played by wheat in global food security. Food Secur. 2013, 5, 291–317. [Google Scholar] [CrossRef]
Han, D.; Liu, S.; Du, Y.; Xie, X.; Fan, L.; Lei, L.; Li, Z.; Yang, H.; Yang, G. Crop Water Content of Winter Wheat Revealed with Sentinel-1 and Sentinel-2 Imagery. Sensors 2019, 19, 4013. [Google Scholar] [CrossRef]
Zhang, J.; Han, W.; Huang, L.; Zhang, Z.; Ma, Y.; Hu, Y. Leaf Chlorophyll Content Estimation of Winter Wheat Based on Visible and Near-Infrared Sensors. Sensors 2016, 16, 437. [Google Scholar] [CrossRef]
Caballero, G.R.; Platzeck, G.; Pezzola, A.; Casella, A.; Winschel, C.; Silva, S.S.; Ludueña, E.; Pasqualotto, N.; Delegido, J. Assessment of Multi-Date Sentinel-1 Polarizations and GLCM Texture Features Capacity for Onion and Sunflower Classification in an Irrigated Valley: An Object Level Approach. Agronomy 2020, 10, 845. [Google Scholar] [CrossRef]
Fernández-Cirelli, A.; Arumí, J.L.; Rivera, D.; Boochs, P.W. Environmental effects of irrigation in arid and semi-arid regions. Chil. J. Agric. Res. 2009, 69 (Suppl. S1), 27–40. [Google Scholar] [CrossRef]
Hanes, J. Biophysical Applications of Satellite Remote Sensing; Springer: Berlin, Germany, 2013. [Google Scholar]
Dong, Q.; Chen, X.; Chen, J.; Zhang, C.; Liu, L.; Cao, X.; Zang, Y.; Zhu, X.; Cui, X. Mapping Winter Wheat in North China Using Sentinel 2A/B Data: A Method Based on Phenology-Time Weighted Dynamic Time Warping. Remote Sens. 2020, 12, 1274. [Google Scholar] [CrossRef]
Ingram, J.S.I.; Gregory, P.J.; Izac, A.M. The role of agronomic research in climate change and food security policy. Agric. Ecosyst. Environ. 2008, 126, 4–12. [Google Scholar] [CrossRef]
Hank, T.B.; Bach, H.; Mauser, W. Using a Remote Sensing-Supported Hydro-Agroecological Model for Field-Scale Simulation of Heterogeneous Crop Growth and Yield: Application for Wheat in Central Europe. Remote Sens. 2015, 7, 3934–3965. [Google Scholar] [CrossRef]
Weiss, M.; Frederic, B.; Smith, G.; Jonckheere, I.; Coppin, P. Review of methods for in situ leaf area index (LAI) determination: Part II. Estimation of LAI, errors and sampling. Agric. For. Meteorol. 2004, 121, 37–53. [Google Scholar] [CrossRef]
Leblanc, S.G.; Chen, J.M.; Fernandes, R.; Deering, D.W.; Conley, A. Methodology comparison for canopy structure parameters extraction from digital hemispherical photography in boreal forests. Agric. For. Meteorol. 2005, 129, 187–207. [Google Scholar] [CrossRef] [Green Version]
Duveiller, G.; Weiss, M.; Baret, F.; Defourny, P. Retrieving wheat Green Area Index during the growing season from optical time series measurements based on neural network radiative transfer inversion. Remote Sens. Environ. 2011, 115, 887–896. [Google Scholar] [CrossRef]
Amin, E.; Verrelst, J.; Rivera-Caicedo, J.P.; Pipia, L.; Ruiz-Verdú, A.; Moreno, J. Prototyping Sentinel-2 green LAI and brown LAI products for cropland monitoring. Remote Sens. Environ. 2021, 255, 112168. [Google Scholar] [CrossRef] [PubMed]
Schlemmer, M.; Gitelson, A.; Schepers, J.; Ferguson, R.; Peng, Y.; Shanahan, J.; Rundquist, D. Remote estimation of nitrogen and chlorophyll contents in maize at leaf and canopy levels. Int. J. Appl. Earth Obs. Geoinf. 2013, 25, 47–54. [Google Scholar] [CrossRef]
Baret, F.; Houles, V.; Guérif, M. Quantification of plant stress using remote sensing observations and crop models: The case of nitrogen management. J. Exp. Bot. 2007, 58, 869–880. [Google Scholar] [CrossRef] [PubMed]
Delloye, C.; Weiss, M.; Defourny, P. Retrieval of the canopy chlorophyll content from Sentinel-2 spectral bands to estimate nitrogen uptake in intensive winter wheat cropping systems. Remote Sens. Environ. 2018, 216, 245–261. [Google Scholar] [CrossRef]
Berger, K.; Verrelst, J.; Féret, J.B.; Wang, Z.; Wocher, M.; Strathmann, M.; Danner, M.; Mauser, W.; Hank, T. Crop nitrogen monitoring: Recent progress and principal developments in the context of imaging spectroscopy missions. Remote Sens. Environ. 2020, 242, 111758. [Google Scholar] [CrossRef]
Peñuelas, J.; Filella, I.; Biel, C.; Serrano, L.; Savé, R. The reflectance at the 950–970 nm region as an indicator of plant water status. Int. J. Remote Sens. 1993, 14, 1887–1905. [Google Scholar] [CrossRef]
Clevers, J.G.P.W.; Kooistra, L.; Schaepman, M.E. Estimating canopy water content using hyperspectral remote sensing data. Int. J. Appl. Earth Obs. Geoinf. 2010, 12, 119–125. [Google Scholar] [CrossRef]
Zhang, C.; Pattey, E.; Liu, J.; Cai, H.; Shang, J.; Dong, T. Retrieving Leaf and Canopy Water Content of Winter Wheat Using Vegetation Water Indices. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2017, 11, 112–126. [Google Scholar] [CrossRef]
Fensholt, R.; Sandholt, I. Derivation of a shortwave infrared water stress index from MODIS near- and shortwave infrared data in a semiarid environment. Remote Sens. Environ. 2003, 87, 111–121. [Google Scholar] [CrossRef]
Raj, R.; Walker, J.P.; Vinod, V.; Pingale, R.; Naik, B.; Jagarlapudi, A. Leaf water content estimation using top-of-canopy airborne hyperspectral data. Int. J. Appl. Earth Obs. Geoinf. 2021, 102, 102393. [Google Scholar] [CrossRef]
Rud, R.; Cohen, Y.; Alchanatis, V.; Levi, A.; Brikman, R.; Shenderey, C.; Heuer, B.; Markovitch, T.; Dar, Z.; Rosen, C.; et al. Crop water stress index derived from multi-year ground and aerial thermal images as an indicator of potato water status. Precis. Agric. 2014, 15, 273–289. [Google Scholar] [CrossRef]
Schott, J.W. Remote Sensing Of The Earth: A Synoptic View. Phys. Today 1989, 42, 72. [Google Scholar] [CrossRef]
Drusch, M.; Del Bello, U.; Carlier, S.; Colin, O.; Fernandez, V.; Gascon, F.; Hoersch, B.; Isola, C.; Laberinti, P.; Martimort, P.; et al. Sentinel-2: ESA’s Optical High-Resolution Mission for GMES Operational Services. Remote Sens. Environ. 2012, 120, 25–36. [Google Scholar] [CrossRef]
Berger, K.; Verrelst, J.; Féret, J.B.; Hank, T.; Wocher, M.; Mauser, W.; Camps-Valls, G. Retrieval of aboveground crop nitrogen content with a hybrid machine learning method. Int. J. Appl. Earth Obs. Geoinf. 2020, 92, 102174. [Google Scholar] [CrossRef]
Danner, M.; Berger, K.; Wocher, M.; Mauser, W.; Hank, T. Efficient RTM-based training of machine learning regression algorithms to quantify biophysical & biochemical traits of agricultural crops. ISPRS J. Photogramm. Remote Sens. 2021, 173, 278–296. [Google Scholar]
Verrelst, J.; Rivera, J.; Camps-Valls, G.; Moreno, J. Recent advances in biophysical parameter retrieval methods—Opportunites for Sentinel-2. In Proceedings of the ESA Living Planet Symposium 2013, Edinbrugh, UK, 9–13 September 2013. [Google Scholar]
Verrelst, J.; Rivera, J.; Tol, C.; Magnani, F.; Mohammed, G.; Moreno, J. Global sensitivity analysis of the SCOPE model: What drives simulated canopy-leaving sun-induced fluorescence? Remote Sens. Environ. 2015, 166, 8–21. [Google Scholar] [CrossRef]
Verrelst, J.; Malenovskỳ, Z.; Van der Tol, C.; Camps-Valls, G.; Gastellu-Etchegorry, J.P.; Lewis, P.; North, P.; Moreno, J. Quantifying vegetation biophysical variables from imaging spectroscopy data: A review on retrieval methods. Surv. Geophys. 2019, 40, 589–629. [Google Scholar] [CrossRef]
Verrelst, J.; Camps-Valls, G.; Muñoz Marí, J.; Rivera, J.; Veroustraete, F.; Clevers, J.; Moreno, J. Optical remote sensing and the retrieval of terrestrial vegetation bio-geophysical properties - A review. ISPRS J. Photogramm. Remote Sens. 2015, 108, 273–290. [Google Scholar] [CrossRef]
Gitelson, A.; Zur, Y.; Chivkunova, O.; Merzlyak, M. Assessing carotenoid content in plant leaves with reflectance spectroscopy. Photochem. Photobiol. 2002, 75, 272–281. [Google Scholar] [CrossRef]
Broge, N.H.; Leblanc, E. Comparing prediction power and stability of broadband and hyperspectral vegetation indices for estimation of green leaf area index and canopy chlorophyll density. Remote Sens. Environ. 2001, 76, 156–172. [Google Scholar] [CrossRef]
Atzberger, C.; Guérif, M.; Baret, F.; Werner, W. Comparative analysis of three chemometric techniques for the spectroradiometric assessment of canopy chlorophyll content in winter wheat. Comput. Electron. Agric. 2010, 73, 165–173. [Google Scholar] [CrossRef]
Mountrakis, G.; Im, J.; Ogole, C. Support vector machines in remote sensing: A review. ISPRS J. Photogramm. Remote Sens. 2011, 66, 247–259. [Google Scholar] [CrossRef]
Verger, A.; Baret, F.; Camacho, F. Optimal modalities for radiative transfer-neural network estimation of canopy biophysical characteristics: Evaluation over an agricultural area with CHRIS/PROBA observations. Remote Sens. Environ. 2011, 115, 415–426. [Google Scholar] [CrossRef]
Jacquemoud, S.; Baret, F.; Andrieu, B.; Danson, F.M.; Jaggard, K. Extraction of vegetation biophysical parameters by inversion of the PROSPECT+SAIL models on sugar beet canopy reflectance data. Application to TM and AVIRIS sensors. Remote Sens. Environ. 1995, 52, 163–172. [Google Scholar] [CrossRef]
Jacquemoud, S.; Verhoef, W.; Baret, F.; Bacour, C.; Zarco-Tejada, P.; Asner, G.; François, C.; Ustin, S. PROSPECT + SAIL models: A review of use for vegetation characterization. Remote Sens. Environ. 2009, 113, S56–S66. [Google Scholar] [CrossRef]
Kimes, D.S.; Nelson, R.F.; Manry, M.T.; Fung, A.K. Attributes of neural networks for extracting continuous vegetation variables from optical and radar measurements. Int. J. Remote Sens. 1998, 19, 2639–2662. [Google Scholar] [CrossRef]
Houborg, R.; Boegh, E. Mapping leaf chlorophyll and leaf area index using inverse and forward canopy reflectance modeling and SPOT reflectance data. Remote Sens. Environ. 2008, 112, 186–202. [Google Scholar] [CrossRef]
Verrelst, J.; Rivera, J.; van der Tol, C.; Magnani, F.; Mohammed, G.; Moreno, J. Global sensitivity Analysis of the A-SCOPE model in support of future FLEX fluorescence retrievals. In Proceedings of the 5th International Workshop on Remote Sensing of Vegetation Fluorescence, Paris, France, 22–24 April 2014. [Google Scholar]
Jacquemoud, S. PROSPECT + SAIL models: A review of use for vegetation characterization. Remote Sens. Environ 2009, 113, 56–66. [Google Scholar] [CrossRef]
Feret, J.B.; François, C.; Asner, G.P.; Gitelson, A.A.; Martin, R.E.; Bidel, L.P.R.; Ustin, S.L.; le Maire, G.; Jacquemoud, S. PROSPECT-4 and 5: Advances in the leaf optical properties model separating photosynthetic pigments. Remote Sens. Environ. 2008, 112, 3030–3043. [Google Scholar] [CrossRef]
Féret, J.B.; Berger, K.; de Boissieu, F.; Malenovskỳ, Z. PROSPECT-PRO for estimating content of nitrogen-containing leaf proteins and other carbon-based constituents. Remote Sens. Environ. 2021, 252, 112173. [Google Scholar] [CrossRef]
Verhoef, W. Light scattering by leaf layers with application to canopy reflectance modeling: The SAIL model. Remote Sens. Environ. 1984, 16, 125–141. [Google Scholar] [CrossRef]
Verhoef, W.; Bach, H. Coupled soil–leaf-canopy and atmosphere radiative transfer modeling to simulate hyperspectral multi-angular surface reflectance and TOA radiance data. Remote Sens. Environ. 2007, 109, 166–182. [Google Scholar] [CrossRef]
Brede, B.; Verrelst, J.; Gastellu-Etchegorry, J.P.; Clevers, J.G.; Goudzwaard, L.; den Ouden, J.; Verbesselt, J.; Herold, M. Assessment of workflow feature selection on forest LAI prediction with sentinel-2A MSI, landsat 7 ETM+ and Landsat 8 OLI. Remote Sens. 2020, 12, 915. [Google Scholar] [CrossRef]
Abdelbaki, A.; Udelhoven, T. A Review of Hybrid Approaches for Quantitative Assessment of Crop Traits Using Optical Remote Sensing: Research Trends and Future Directions. Remote Sens. 2022, 14, 3515. [Google Scholar] [CrossRef]
Rasmussen, C.E.; Williams, C.K.I. Gaussian Processes for Machine Learning; The MIT Press: New York, NY, USA, 2006. [Google Scholar]
Verrelst, J.; Muñoz, J.; Alonso, L.; Delegido, J.; Rivera, J.P.; Camps-Valls, G.; Moreno, J. Machine learning regression algorithms for biophysical parameter retrieval: Opportunities for Sentinel-2 and-3. Remote Sens. Environ. 2012, 118, 127–139. [Google Scholar] [CrossRef]
Verrelst, J.; Rivera, J.; Veroustraete, F.; Muñoz Marí, J.; Clevers, J.; Camps-Valls, G.; Moreno, J. Experimental Sentinel-2 LAI estimation using parametric, non-parametric and physical retrieval methods—A comparison. ISPRS J. Photogramm. Remote Sens. 2015, 108, 260–272. [Google Scholar] [CrossRef]
Sinha, S.K.; Padalia, H.; Dasgupta, A.; Verrelst, J.; Rivera, J.P. Estimation of leaf area index using PROSAIL based LUT inversion, MLRA-GPR and empirical models: Case study of tropical deciduous forest plantation, North India. Int. J. Appl. Earth Obs. Geoinf. 2020, 86, 102027. [Google Scholar] [CrossRef]
Verrelst, J.; Alonso, L.; Camps-Valls, G.; Delegido, J.; Moreno, J. Retrieval of Vegetation Biophysical Parameters Using Gaussian Process Techniques. IEEE Trans. Geosci. Remote Sens. 2012, 50, 1832–1843. [Google Scholar] [CrossRef]
Van Wittenberghe, S.; Verrelst, J.; Rivera, J.P.; Alonso, L.; Moreno, J.; Samson, R. Gaussian processes retrieval of leaf parameters from a multi-species reflectance, absorbance and fluorescence dataset. J. Photochem. Photobiol. B Biol. 2014, 134, 37–48. [Google Scholar] [CrossRef]
Verrelst, J.; Rivera, J.P.; Gitelson, A.; Delegido, J.; Moreno, J.; Camps-Valls, G. Spectral band selection for vegetation properties retrieval using Gaussian processes regression. Int. J. Appl. Earth Obs. Geoinf. 2016, 52, 554–567. [Google Scholar] [CrossRef]
Xie, R.; Darvishzadeh, R.; Skidmore, A.K.; Heurich, M.; Holzwarth, S.; Gara, T.W.; Reusen, I. Mapping leaf area index in a mixed temperate forest using Fenix airborne hyperspectral data and Gaussian processes regression. Int. J. Appl. Earth Obs. Geoinf. 2021, 95, 102242. [Google Scholar] [CrossRef]
Camps-Valls, G.; Verrelst, J.; Munoz-Mari, J.; Laparra, V.; Mateo-Jimenez, F.; Gomez-Dans, J. A survey on Gaussian processes for earth-observation data analysis: A comprehensive investigation. IEEE Geosci. Remote Sens. Mag. 2016, 4, 58–78. [Google Scholar] [CrossRef]
Pascual-Venteo, A.B.; Portalés, E.; Berger, K.; Tagliabue, G.; Garcia, J.L.; Pérez-Suay, A.; Rivera-Caicedo, J.P.; Verrelst, J. Prototyping Crop Traits Retrieval Models for CHIME: Dimensionality Reduction Strategies Applied to PRISMA Data. Remote Sens. 2022, 14, 2448. [Google Scholar] [CrossRef]
Estévez, J.; Berger, K.; Vicent, J.; Rivera-Caicedo, J.P.; Wocher, M.; Verrelst, J. Top-of-Atmosphere Retrieval of Multiple Crop Traits Using Variational Heteroscedastic Gaussian Processes within a Hybrid Workflow. Remote Sens. 2021, 13, 1589. [Google Scholar] [CrossRef]
Salinero-Delgado, M.; Estévez, J.; Pipia, L.; Belda, S.; Berger, K.; Paredes Gómez, V.; Verrelst, J. Monitoring Cropland Phenology on Google Earth Engine Using Gaussian Process Regression. Remote Sens. 2021, 14, 146. [Google Scholar] [CrossRef]
Reyes-Muñoz, P.; Pipia, L.; Salinero-Delgado, M.; Belda, S.; Berger, K.; Estévez, J.; Morata, M.; Rivera-Caicedo, J.P.; Verrelst, J. Quantifying Fundamental Vegetation Traits over Europe Using the Sentinel-3 OLCI Catalogue in Google Earth Engine. Remote Sens. 2022, 14, 1347. [Google Scholar] [CrossRef]
Estévez, J.; Salinero-Delgado, M.; Berger, K.; Pipia, L.; Rivera-Caicedo, J.P.; Wocher, M.; Reyes-Muñoz, P.; Tagliabue, G.; Boschetti, M.; Verrelst, J. Gaussian processes retrieval of crop traits in Google Earth Engine based on Sentinel-2 top-of-atmosphere data. Remote Sens. Environ. 2022, 273, 112958. [Google Scholar] [CrossRef]
Adeluyi, O.; Harris, A.; Verrelst, J.; Foster, T.; Clay, G.D. Estimating the phenological dynamics of irrigated rice leaf area index using the combination of PROSAIL and Gaussian Process Regression. Int. J. Appl. Earth Obs. Geoinf. 2021, 102, 102454. [Google Scholar] [CrossRef]
Perich, G.; Aasen, H.; Verrelst, J.; Argento, F.; Walter, A.; Liebisch, F. Crop Nitrogen Retrieval Methods for Simulated Sentinel-2 Data Using In-Field Spectrometer Data. Remote Sens. 2021, 13, 2404. [Google Scholar] [CrossRef]
Settles, B. Active Learning Literature Survey; University of Wisconsin-Madison Department of Computer Sciences: Madison, WI, USA, 2009. [Google Scholar]
Verrelst, J.; Dethier, S.; Rivera, J.P.; Munoz-Mari, J.; Camps-Valls, G.; Moreno, J. Active learning methods for efficient hybrid biophysical variable retrieval. IEEE Geosci. Remote Sens. Lett. 2016, 13, 1012–1016. [Google Scholar] [CrossRef]
Berger, K.; Rivera Caicedo, J.P.; Martino, L.; Wocher, M.; Hank, T.; Verrelst, J. A Survey of Active Learning for Quantifying Vegetation Traits from Terrestrial Earth Observation Data. Remote Sens. 2021, 13, 287. [Google Scholar] [CrossRef]
Delegido, J.; Verrelst, J.; Rivera, J.P.; Ruiz-Verdú, A.; Moreno, J. Brown and green LAI mapping through spectral indices. Int. J. Appl. Earth Obs. Geoinf. 2015, 35, 350–358. [Google Scholar] [CrossRef]
Pasqualotto, N.; D’Urso, G.; Bolognesi, S.F.; Belfiore, O.R.; Van Wittenberghe, S.; Delegido, J.; Pezzola, A.; Winschel, C.; Moreno, J. Retrieval of Evapotranspiration from Sentinel-2: Comparison of Vegetation Indices, Semi-Empirical Models and SNAP Biophysical Processor Approach. Agronomy 2019, 9, 663. [Google Scholar] [CrossRef]
Locherer, M.; Hank, T.; Danner, M.; Mauser, W. Retrieval of Seasonal Leaf Area Index from Simulated EnMAP Data through Optimized LUT-Based Inversion of the PROSAIL Model. Remote Sens. 2015, 7, 10321–10346. [Google Scholar] [CrossRef]
Upreti, D.; Huang, W.; Kong, W.; Pascucci, S.; Pignatti, S.; Zhou, X.; Ye, H.; Casa, R. A comparison of hybrid machine learning algorithms for the retrieval of wheat biophysical variables from sentinel-2. Remote Sens. 2019, 11, 481. [Google Scholar] [CrossRef]
Verrelst, J.; Rivera-Caicedo, J.P.; Reyes-Muñoz, P.; Morata, M.; Amin, E.; Tagliabue, G.; Panigada, C.; Hank, T.; Berger, K. Mapping landscape canopy nitrogen content from space using PRISMA data. ISPRS J. Photogramm. Remote Sens. 2021, 178, 382–395. [Google Scholar] [CrossRef]
Liang, L.; Geng, D.; Yan, J.; Qiu, S.; Di, L.; Wang, S.; Xu, L.; Wang, L.; Kang, J.; Li, L. Estimating Crop LAI Using Spectral Feature Extraction and the Hybrid Inversion Method. Remote Sens. 2020, 12, 3534. [Google Scholar] [CrossRef]
Verrelst, J.; Berger, K.; Rivera-Caicedo, J.P. Intelligent sampling for vegetation nitrogen mapping based on hybrid machine learning algorithms. IEEE Geosci. Remote Sens. Lett. 2020, 18, 2038–2042. [Google Scholar] [CrossRef]
Douak, F.; Melgani, F.; Benoudjit, N. Kernel ridge regression with active learning for wind speed prediction. Appl. Energy 2013, 103, 328–340. [Google Scholar] [CrossRef]
Zadoks, J.C.; Chang, T.T.; Konzak, C.F. A decimal code for the growth stages of cereals. Weed Res. 1974, 14, 415–421. [Google Scholar] [CrossRef]
Confalonieri, R.; Foi, M.; Casa, R.; Aquaro, S.; Tona, E.; Peterle, M.; Boldini, A.; De Carli, G.; Ferrari, A.; Finotto, G.; et al. Development of an app for estimating leaf area index using a smartphone. Trueness and precision determination and comparison with other indirect methods. Comput. Electron. Agric. 2013, 96, 67–74. [Google Scholar] [CrossRef]
Patrignani, A.; Ochsner, T.E. Canopeo: A Powerful New Tool for Measuring Fractional Green Canopy Cover. Agron. J. 2015, 107, 2312–2320. [Google Scholar] [CrossRef]
Gamiely, S.; Randle, W.M.; Mills, H.A.; Smittle, D.A. A Rapid and Nondestructive Method for Estimating Leaf Area of Onions. HortScience 1991, 26, 206. [Google Scholar] [CrossRef]
Danson, F.M.; Steven, M.D.; Malthus, T.J.; Clark, J.A. High-spectral resolution data for determining leaf water content. Int. J. Remote Sens. 1992, 13, 461–470. [Google Scholar] [CrossRef]
Louis, J.; Debaecker, V.; Pflug, B.; Main-Knorn, M.; Bieniarz, J.; Mueller-Wilm, U.; Cadau, E.; Gascon, F. SENTINEL-2 SEN2COR: L2A Processor for Users. In Proceedings of the ESA Living Planet Symposium 2016, Prague, Czech Republic, 9–13 May 2016; Ouwehand, L., Ed.; Spacebooks Online; ESA Special Publications (on CD). Volume SP-740, pp. 1–8. [Google Scholar]
Berger, K.; Atzberger, C.; Danner, M.; Wocher, M.; Mauser, W.; Hank, T. Model-Based Optimization of Spectral Sampling for the Retrieval of Crop Variables with the PROSAIL Model. Remote Sens. 2018, 10, 2063. [Google Scholar] [CrossRef]
Verrelst, J.; Rivera, J.; Alonso, L.; Moreno, J. ARTMO: An Automated Radiative Transfer Models Operator toolbox for automated retrieval of biophysical parameters through model inversion. In Proceedings of the EARSeL 7th SIG-Imaging Spectroscopy Workshop, Edinburgh, UK, 11–13 April 2011. [Google Scholar]
Verrelst, J.; Romijn, E.; Kooistra, L. Mapping Vegetation Density in a Heterogeneous River Floodplain Ecosystem Using Pointable CHRIS/PROBA Data. Remote Sens. 2012, 4, 2866–2889. [Google Scholar] [CrossRef]
Caicedo, J.; Verrelst, J.; Munoz-Mari, J.; Moreno, J.; Camps-Valls, G. Toward a semiautomatic machine learning retrieval of biophysical parameters. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2014, 7, 1249–1259. [Google Scholar] [CrossRef]
Smit, M.A.; Singels, A. The response of sugarcane canopy development to water stress. Field Crops Res. 2006, 98, 91–97. [Google Scholar] [CrossRef]
Jacquemoud, S.; Ustin, S.L.; Verdebout, J.; Schmuck, G.; Andreoli, G.; Hosgood, B. Estimating leaf biochemistry using the PROSPECT leaf optical properties model. Remote Sens. Environ. 1996, 56, 194–202. [Google Scholar] [CrossRef]
Féret, J.B.; Gitelson, A.A.; Noble, S.D.; Jacquemoud, S. PROSPECT-D: Towards modeling leaf optical properties through a complete lifecycle. Remote Sens. Environ. 2017, 193, 204–215. [Google Scholar] [CrossRef] [Green Version]
Ma, X.; Lu, L.; Ding, J.; Zhang, F.; He, B. Estimating Fractional Vegetation Cover of Row Crops from High Spatial Resolution Image. Remote Sens. 2021, 13, 3874. [Google Scholar] [CrossRef]
Danner, M.; Berger, K.; Wocher, M.; Mauser, W.; Hank, T. Fitted PROSAIL parameterization of leaf inclinations, water content and brown pigment content for winter wheat and maize canopies. Remote Sens. 2019, 11, 1150. [Google Scholar] [CrossRef]
Wocher, M.; Berger, K.; Danner, M.; Mauser, W.; Hank, T. Physically-based retrieval of canopy equivalent water thickness using hyperspectral data. Remote Sens. 2018, 10, 1924. [Google Scholar] [CrossRef]
Povey, A.C.; Grainger, R.G. Known and unknown unknowns: Uncertainty estimation in satellite remote sensing. Atmos. Meas. Tech. 2015, 8, 4699–4718. [Google Scholar] [CrossRef]
Ali, Z.; Merrium, S.; Habib-ur Rahman, M.; Hakeem, S.; Saddique, M.A.B.; Sher, M.A. Wetting mechanism and morphological adaptation; leaf rolling enhancing atmospheric water acquisition in wheat crop—A review. Environ. Sci. Pollut. Res. 2022, 29, 30967–30985. [Google Scholar] [CrossRef] [PubMed]
Gil, H.M.; Salort, J.B.; Llompart, J.C.; Sans, J.F.; Carbó, M.R.; León, J.G. Eficiencia en el uso del agua por las plantas. Investig. Geográficas 2007, 63–84. [Google Scholar] [CrossRef]
Sánchez, R.; Pezzola, N.; Cepeda, J. Caracterización edafoclimátia del área de influencia del INTA E.E.A Hilario Ascasubi. Boletín Divulgación 1998, 18, 72. [Google Scholar]
Delegido, J.; Verrelst, J.; Alonso, L.; Moreno, J. Evaluation of sentinel-2 red-edge bands for empirical estimation of green LAI and chlorophyll content. Sensors 2011, 11, 7063–7081. [Google Scholar] [CrossRef]
Salisbury, F.B.; Ross, C.W. Fisiología Vegetal; Grupo Editorial Iberoamérica: Mexico City, Mexico, 1994. [Google Scholar]
Delegido, J.; Alonso, L.; González, G.; Moreno, J. Estimating chlorophyll content of crops from hyperspectral data using a normalized area over reflectance curve (NAOC). Int. J. Appl. Earth Obs. Geoinf. 2010, 12, 165–174. [Google Scholar] [CrossRef]
Elarab, M.; Ticlavilca, A.M.; Torres-Rua, A.F.; Maslova, I.; McKee, M. Estimating chlorophyll with thermal and broadband multispectral high resolution imagery from an unmanned aerial system using relevance vector machines for precision agriculture. Int. J. Appl. Earth Obs. Geoinf. 2015, 43, 32–42. [Google Scholar] [CrossRef]
Cammarano, D.; Fitzgerald, G.; Basso, B.; O’Leary, G.; Chen, D.; Grace, P.; Fiorentino, C. Use of the Canopy Chlorophyl Content Index (CCCI) for Remote Estimation of Wheat Nitrogen Content in Rainfed Environments. Agron. J. 2011, 103, 1597–1603. [Google Scholar] [CrossRef]
Irmak, S. Evapotranspiration. In Encyclopedia of Ecology; Academic Press: Cambridge, MA, USA, 2008; pp. 1432–1438. [Google Scholar] [CrossRef]
Zeleke, K.T.; Wade, L.J. Evapotranspiration Estimation Using Soil Water Balance, Weather and Crop Data. In Evapotranspiration-Remote Sensing and Modeling; INTECH: Rijeka, Croatia, 2012; pp. 41–58. [Google Scholar]
Darouich, H.; Ramos, T.B.; Pereira, L.S.; Rabino, D.; Bagagiolo, G.; Capello, G.; Simionesei, L.; Cavallo, E.; Biddoccu, M. Water Use and Soil Water Balance of Mediterranean Vineyards under Rainfed and Drip Irrigation Management: Evapotranspiration Partition and Soil Management Modelling for Resource Conservation. Water 2022, 14, 554. [Google Scholar] [CrossRef]
Verrelst, J.; Rivera, J.; Moreno, J.; Camps-Valls, G. Gaussian processes uncertainty estimates in experimental Sentinel-2 LAI and leaf chlorophyll content retrieval. ISPRS J. Photogramm. Remote Sens. 2013, 86, 157–167. [Google Scholar] [CrossRef]
Machwitz, M.; Pieruschka, R.; Berger, K.; Schlerf, M.; Aasen, H.; Fahrner, S.; Jiménez-Berni, J.; Baret, F.; Rascher, U. Bridging the Gap Between Remote Sensing and Plant Phenotyping—Challenges and Opportunities for the Next, Generation of Sustainable Agriculture. Front. Plant Sci. 2021, 2334. [Google Scholar] [CrossRef] [PubMed]
Pipia, L.; Amin, E.; Belda, S.; Salinero-Delgado, M.; Verrelst, J. Green LAI Mapping and Cloud Gap-Filling Using Gaussian Process Regression in Google Earth Engine. Remote Sens. 2021, 13, 403. [Google Scholar] [CrossRef] [PubMed]

Figure 1. Overview of the Bonaerense Valley of Colorado River study site with test fields for the winter wheat campaign of the year 2020. True colour S2 image (R = B4, G = B3, B = B2) of 27 December 2020. Reference system: WGS84 (EPSG 4326).

Figure 2. Layout of measurement design for the 2020 campaign at the BVCR study site: three ESUs were defined per wheat paddock and sampling approach for each elementary sampling unit, partly adapted from [69]. Reference system: WGS84 (EPSG 4326).

Figure 3. Photographic documentation of wheat crop growing period from August to November 2020, BVCR campaign, corresponding to LAI, FVC, C

_{a b}

, and AGFB sampling dates. With (a) seedling stage at 10 August 2020; (b) tillering stage at 4 September 2020; (c) tillering stage at 17 September 2020; (d) tillering stage at 2 Ocotber 2020; (e) ear emergence from boot at 19 October 2020; (f) anthesis stage at 2 November 2020; (g) dough development at 30 November 2020, first appearance of senescence; (h) ripening stage at 16 December 2020, complete senescence.

Figure 3. Photographic documentation of wheat crop growing period from August to November 2020, BVCR campaign, corresponding to LAI, FVC, C

_{a b}

Figure 5. Goodness-of-fit results (RMSE, R

^{2}

) using AL (EBD) against validation data. (a) LAI model; (b) CCC model; (c) VWC model.

Figure 5. Goodness-of-fit results (RMSE, R

^{2}

) using AL (EBD) against validation data. (a) LAI model; (b) CCC model; (c) VWC model.

Figure 6. Measuredvs. estimated wheat traits along 1:1-line including uncertainty intervals, using the EBD-optimised training data set. (a) LAI model estimates; (b) CCC model estimates; (c) VWC model estimates.

Figure 7. LAI [m

^{2}

^{- 2}

] of different crops in the study site, retrieved by the GPR model using simulated vegetation reflectance spectra, in-situ traits measurements and S2 surface multi-spectral reflectance data. Data obtained from BVCR wheat campaign 2020. (a) wheat tillering stage at 29 August 2020; (b) wheat booting stage at 28 September 2020; (c) wheat anthesis-flowering stage at 2 November 2020; (d) wheat dough development stage at 27 November 2020; (e) wheat ripening stage at 7 December 2020; (f) harvested wheat at 21 January 2021.

Figure 7. LAI [m

^{2}

^{- 2}

Figure 8. CCC [g m

^{- 2}

] of different crops in the study site, retrieved by the GPR model using simulated vegetation reflectance spectra, in-situ biophysical/biochemical measurements and S2 surface multi-spectral reflectance data. Data obtained from BVCR wheat campaign 2020. (a) wheat tillering stage at 29 August 2020; (b) wheat booting stage at 28 September 2020; (c) wheat anthesis-flowering stage at 2 November 2020; (d) wheat dough development stage at 27 November 2020; (e) wheat ripening stage at 7 December 2020; (f) harvested wheat at 21 January 2021.

Figure 8. CCC [g m

^{- 2}

Figure 9. VWC [g m

^{- 2}

Figure 9. VWC [g m

^{- 2}

Figure 10. Seasonal evolution of wheat cropland over the three paddocks at the BVCR study sites described by LAI, CCC, and VWC mean values of nice ESUs within the crop limits and the associated uncertainty, plotted as vertical bars. (a) LAI estimates; (b) LAI uncertainty (SD); (c) CCC estimates; (d) CCC uncertainty (SD); (e) VWC estimates; (f) VWC uncertainty (SD).

Figure 11. Temporal evolution of wheat crop over the three paddocks at the BVCR study sites described by FVC mean measured values of nine ESUs within the crop limits and the associated SD, which is plotted as shadowed areas.

Figure 12. Temporal evolution of wheat cropland over the three paddocks at the BVCR study sites described by LAI, CCC and VWC mean values of nine ESUs within the crop limits and the associated SD, which is plotted as vertical bars. (a) LAI vs. CCC temporal estimates; (b) LAI temporal estimates vs. FVC in-situ measured values; (c) LAI vs. VWC temporal estimates.

Table 1. Parameter ranges for PROSPECT-PRO and 4SAIL (PROSAIL-PRO). Specified ranges are uniformly distributed in LHS sampling.

Leaf Optical Properties			Canopy Reflectance Model
PROSPECT-PRO Parameters	Notation [Unit]	Range	4SAIL Parameters	Notation [Unit]	Range
Leaf chlorophyll a+b content	Cab [µg cm $^{- 2}$ ]	5–75	Leaf area index	LAI [m $^{2}$ m $^{- 2}$ ]	0.1–7.0
Leaf structure parameter	Nstruct, no dim.	1.0–2.0	Average leaf inclination angle	ALIA [ $^{\circ}$ ]	30–70
Leaf carotenoid content	Cxc [µg cm $^{- 2}$ ]	0–15	Soil brightness	soil, no dim.	0–1
Leaf equivalent water thickness	EWT [cm]	0.0002–0.05	Sun zenith angle	SZA [ $^{\circ}$ ]	20–40
Carbon-based constituents	CBC [g cm $^{- 2}$ ]	0.001–0.01	Hot spot effect	Hot [m m $^{- 1}$ ]	0.01
Leaf anthocyanin content	Canth [µg cm $^{- 2}$ ]	0–2	Observer zenith angle	OZA [ $^{\circ}$ ]	0
Leaf protein content	Cp [µg cm $^{- 2}$ ]	0.001–0.0025	Diffuse/direct radiation	DDR [%]	80
Leaf mass per area	Cm [µg cm $^{- 2}$ ]	0.0001–0.03	Relative azimuth angle	rAA [ $^{\circ}$ ]	0
Brown pigment content	Cbrown, no dim.	0

Table 2. Database of in-situ wheat traits with dates, ranges, mean values as well as standard deviations (SD) of the measurements.

Wheat Variable	Date	Range	Mean	SD
LAI (m $^{2}$ m $^{- 2}$ )	03-09-2020	0.16–0.30	0.23	0.05
	17-09-2020	0.56–1.54	0.94	0.29
	02-10-2020	1.59–3.81	2.57	0.66
	19-10-2020	1.53–3.27	2.62	0.51
	02-11-2020	2.78–5.05	4.12	0.63
	16-11-2020	3.31–5.39	4.02	0.80
	30-11-2020	3.29–4.75	4.08	0.50
	16-12-2020	3.97–5.64	4.68	0.43
FVC (%)	10-08-2020	6.2–9.1	7.63	1.00
	03-09-2020	23.0–48.0	34.94	8.32
	17-09-2020	22.1–80.2	44.21	22.73
	02-10-2020	32.7–69.2	48.30	13.28
	19-10-2020	23.6–69.5	46.15	12.97
	02-11-2020	11.1–29.3	20.92	5.07
	16-11-2020	74.4–92.0	87.44	4.83
	30-11-2020	80.0–90.8	88.22	3.06
C $_{a b}$ (µg cm $^{- 2}$ )	03-09-2020	38.23–44.82	41.96	2.23
	17-09-2020	36.49–52.61	42.21	4.91
	02-10-2020	38.12–52.02	45.33	4.46
	19-10-2020	39.64–45.12	43.08	1.84
	02-11-2020	33.63–42.44	38.29	2.89
	16-11-2020	35.72–44.92	39.41	2.85
	30-11-2020	13.93–48.31	35.32	9.59
AGFB (g)	03-09-2020	15–25	19.67	3.59
	17-09-2020	31–54	42.00	6.83
	02-10-2020	76–175	112.67	27.23
	19-10-2020	47–94	66.00	13.61
	02-11-2020	131–296	213.67	42.41
	16-11-2020	57–101	81.78	15.84
	30-11-2020	73–184	121.60	38.45
AGDB (g)	03-09-2020	2.00–6.00	3.56	1.17
	17-09-2020	9.00–15.00	11.67	2.00
	02-10-2020	23.00–48.00	33.22	8.23
	19-10-2020	8.00–15.00	12.11	2.42
	02-11-2020	38–62	45.67	7.86
	16-11-2020	17–32	26.00	4.22
	30-11-2020	23–70	46.67	14.04

Table 3. Dates, ranges, mean values, and standard deviations (SD) for calculated variables (CCC and VWC).

Wheat Variable	Date	Range	Mean	SD
CCC (g m $^{- 2}$ )	03-09-2020	0.06–0.12	0.10	0.02
	17-09-2020	0.22–0.71	0.41	0.16
	02-10-2020	0.74–1.61	1.16	0.29
	19-10-2020	0.60–1.44	1.13	0.23
	02-11-2020	1.18–1.74	1.56	0.17
	16-11-2020	1.22–2.10	1.59	0.34
	30-11-2020	0.64–1.70	1.40	0.29
VWC (g m $^{- 2}$ )	03-09-2020	207–455	284	70
	17-09-2020	315–1554	666	360
	02-10-2020	868–3021	1996	693
	19-10-2020	494–2240	1315	554
	02-11-2020	944–3083	1777	629
	16-11-2020	1481–3275	2399	589
	30-11-2020	2016–5079	3287	1113

Table 4. Winter wheat field campaign observations with in-situ sampling dates, crop growth stages and details of the field observations.

In-Situ Measurements Date	Wheat Growth Stage	Field Observations
10-08-2020	Seedling growth Z1.3—Three leaves emerged	Plant density: 248 plants m $^{- 2}$ (on average), previous crop: sunflower for seed
03-09-2020	Tillering, 2–3 tillers, Z2.3—Main stem and three tillers	Leaves per tiller 2 + 1 flag leaf
17-09-2020	Tillering, 4 tillers Z2.4—Main stem and four tillers	Leaves per tiller: 3 Irrigation date: 17/09/2020 Fertilization date: 16/09/2020
02-10-2020	Tillering, 5 tillers 4—Booting Z4.3—Boots just visible swollen	Plants height 22 cm from the base to the second node Leaves per tiller: 4
19-10-2020	Ear emergence from boot Z5.5—Ear half emerged	Plants height 71 cm (on average) Plants stem nodes: 5
02-11-2020	Anthesis (flowering) Z6.1—Beginning of anthesis (few anthers at the middle of ear)	Plants height 80 cm (on average)
16-11-2020	Milk development Z7.5—Medium milk	Plants height 80 cm (on average) Start of the senescence
30-11-2020	Dough development Z8.7—Hard dough	Senescence process
16-12-2020	Ripening Z9.7—Seed not dormant	Complete senescence Distance between rows: 19 cm Number of ears at 0.50 cm: 72 on average

Table 5. Field campaign and Sentinel-2 acquisition dates.

In-Situ Measurements Date	S2 Acquisition	± $Δ$ Days
10-08-2020	NA	NA
03-09-2020	29-08-2020	−5
17-09-2020	18-09-2020	+1
02-10-2020	28-09-2020	−4
19-10-2020	13-10-2020	−6
02-11-2020	02-11-2020	0
16-11-2020	17-11-2020	+1
30-11-2020	NA	NA
16-12-2020	22-12-2020	+6

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Caballero, G.; Pezzola, A.; Winschel, C.; Casella, A.; Sanchez Angonova, P.; Rivera-Caicedo, J.P.; Berger, K.; Verrelst, J.; Delegido, J. Seasonal Mapping of Irrigated Winter Wheat Traits in Argentina with a Hybrid Retrieval Workflow Using Sentinel-2 Imagery. Remote Sens. 2022, 14, 4531. https://doi.org/10.3390/rs14184531

AMA Style

Caballero G, Pezzola A, Winschel C, Casella A, Sanchez Angonova P, Rivera-Caicedo JP, Berger K, Verrelst J, Delegido J. Seasonal Mapping of Irrigated Winter Wheat Traits in Argentina with a Hybrid Retrieval Workflow Using Sentinel-2 Imagery. Remote Sensing. 2022; 14(18):4531. https://doi.org/10.3390/rs14184531

Chicago/Turabian Style

Caballero, Gabriel, Alejandro Pezzola, Cristina Winschel, Alejandra Casella, Paolo Sanchez Angonova, Juan Pablo Rivera-Caicedo, Katja Berger, Jochem Verrelst, and Jesus Delegido. 2022. "Seasonal Mapping of Irrigated Winter Wheat Traits in Argentina with a Hybrid Retrieval Workflow Using Sentinel-2 Imagery" Remote Sensing 14, no. 18: 4531. https://doi.org/10.3390/rs14184531

APA Style

Caballero, G., Pezzola, A., Winschel, C., Casella, A., Sanchez Angonova, P., Rivera-Caicedo, J. P., Berger, K., Verrelst, J., & Delegido, J. (2022). Seasonal Mapping of Irrigated Winter Wheat Traits in Argentina with a Hybrid Retrieval Workflow Using Sentinel-2 Imagery. Remote Sensing, 14(18), 4531. https://doi.org/10.3390/rs14184531

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Seasonal Mapping of Irrigated Winter Wheat Traits in Argentina with a Hybrid Retrieval Workflow Using Sentinel-2 Imagery

Abstract

1. Introduction

2. Materials and Methods

2.1. Generation of Training Data Sets

2.2. Gaussian Processes Regression

2.3. Active Learning Principles

2.4. Study Site

2.5. Wheat Crop Experimental Design

Wheat Crop Management and Field Data Collection

2.6. Sentinel-2 Image Acquisitions

2.7. Delineation of the Hybrid Retrieval Workflow

3. Results

3.1. Optimized Sample Selection for LAI, CCC and VWC Modeling

3.2. Lai, CCC, and VWC Mapping

3.3. Wheat Phenology Based on Multi-Temporal LAI Maps

3.4. Seasonal Analysis of Retrieved Traits

3.4.1. In-Situ Measured FVC Time Series Analysis of Irrigated Winter Wheat Crops

3.4.2. Seasonal Analysis of S2-Retrieved CCC and LAI

3.4.3. Seasonal Analysis of S2-Retrieved VWC and LAI

4. Discussion

4.1. Advantages and Limitations of Coupled RTMs

4.2. Performance of Hybrid GPR Models

4.3. Potential of Seasonal Traits Mapping for Wheat Agronomic Management

4.4. Study Limitations

4.5. Opportunities for Operational Monitoring

5. Conclusions

Author Contributions

Funding

Data Availability Statement

Acknowledgments

Conflicts of Interest

Appendix A. EBD-Reduced Final Training Dataset versus Validation Dataset

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI