Exploitation of SAR and Optical Sentinel Data to Detect Rice Crop and Estimate Seasonal Dynamics of Leaf Area Index

This paper presents and evaluates multitemporal LAI estimates derived from Sentinel-2A data on rice cultivated area identified using time series of Sentinel-1A images over the main European rice districts for the 2016 crop season. This study combines the information conveyed by Sentinel-1A and Sentinel-2A into a high-resolution LAI retrieval chain. Rice crop was detected using an operational multi-temporal rule-based algorithm, and LAI estimates were obtained by inverting the PROSAIL radiative transfer model with Gaussian process regression. Direct validation was performed with in situ LAI measurements acquired in coordinated field campaigns in three countries (Italy, Spain and Greece). Results showed high consistency between estimates and ground measurements, revealing high correlations (R2 > 0.93) and good accuracies (RMSE < 0.83, rRMSEm < 23.6% and rRMSEr < 16.6%) in all cases. Sentinel-2A estimates were compared with Landsat-8 showing high spatial consistency between estimates over the three areas. The possibility to exploit seasonally-updated crop mask exploiting Sentinel-1A data and the temporal consistency between Sentinel-2A and Landsat-7/8 LAI time series demonstrates the feasibility of deriving operationally high spatial-temporal decametric multi-sensor LAI time series useful for crop monitoring.


Introduction
Monitoring Earth's vegetated surfaces is a challenging issue in which the green leaf area index (LAI) has been recognized as an essential climate variable (ECV) by the Global Terrestrial Observation System (GTOS) and the Global Climate Observing System (GCOS) [1].Defined as half of the total green leaf area per unit ground surface area [2], LAI plays a key role in land surface models [3] and can be used to address various agricultural issues, such as rice crop monitoring, yield forecasting and crop management [4].
Crop monitoring is necessary to identify the onset of stress conditions, which require actions for reducing their impact on crop yield.Anomaly drops in canopy LAI are indicators of plant stress conditions, which could lead to a decrease of plant production and to an increased senescence rate [5].
Rice yield forecasting is a crucial task for management and planning, and it can be addressed with both statistical and mathematical modeling approaches.Statistical approaches directly link crop biophysical variables, such as LAI, to crop yield.LAI is indeed recognized as the major morphological parameter for determining crop growth, and it is strongly correlated with crop productivity [6].Crop models are able to simulate rice growing and are used to provide indications on crop status and to predict yields over large areas [7,8].However, crop models require information on soil, meteorological variable, crop parameters and management practices, which are not always available or practical to be obtained in a spatially distributed way and during the season [9].Hence, the best way to simulate the real spatial-temporal differences in crop development of fields sowed the same day with the same variety is to assimilate exogenous observation, such as LAI maps, in the modeling solution [10].Accurate estimation of LAI has been shown to improve the accuracy of grain yield estimates [11], and an operational application of this workflow for rice was successfully demonstrated in Asia in the framework of the RIICE (Remote sensing-based Information and Insurance for Crops in Emerging economies) project (http://www.riice.org/)where rice yield is estimated from the Oryza2000 model by assimilating LAI maps derived from synthetic aperture radar (SAR) images [12].
As far as precise crop management is concerned, LAI data have been found to be useful in new approaches for the determination of nitrogen concentration dilution curves, which are traditionally based on plant dry matter (PDM) estimation [13,14].Remotely-sensed data at decametric resolution (e.g., Sentinel-2A) are the sole source of information available to provide high-resolution (HR) LAI estimation on wide areas to exploit the nitrogen concentration dilution curve approach for optimal crop fertilization.
In the last few years, the scientific community has made big efforts with the goal of providing reliable and accurate LAI estimates at local scales taking advantage of unmanned aerial vehicles and high-resolution sensors, such as Landsat [15] and SPOT5 [16][17][18].The recently launched Sentinel-2A satellite [19] provides well-suited spectral and temporal data for LAI retrieval at high-resolution in near real time, useful for assessing crop status and providing support in agro-practices at the parcel level.Many methods have been proposed and implemented in retrieval chains from Earth observation (EO) data to derive LAI estimates [20].Empirical parametric algorithms have been developed based on calibrated relationships between vegetation indices and canopy biophysical variables [21].On the other hand, non-parametric algorithms estimate biophysical variables using a training database containing pairs of the biophysical parameter and the associated spectral data [22].In statistical approaches, concomitant in situ measurements of the biophysical parameter of interest and the associated spectral data from remote sensing platforms are used as a training database, whereas the hybrid approaches rely on a database generated by radiative transfer models (RTMs).RTMs describe the interactions between the incoming radiation, canopy elements and the background soil surface.The PROSAIL RTM has become one of the most popular and widely-used RTMs due to its consistency and robustness [23].Hybrid approaches retrieve LAI by inverting RTM through machine learning techniques with a large number of methods [24,25].Among them, neural networks (NN) [26] and kernel-based methods [27,28] are the most popular and used regression methods in remote sensing.State-of-the-art kernel-based methods, such as Gaussian process regression (GPR) [29], provided encouraging results in the framework of biophysical parameter estimation outperforming other kernel methods and NN [18,30].
This study shows Sentinel-2A LAI estimates produced in the framework of the ERMES (an Earth obseRvation Model based RicE information Service) project (http://www.ermes-fp7space.eu/),which aims to develop a prototype of Copernicus down-stream services assimilating EO and in situ data on crop modeling dedicated to the rice sector.In this framework, the main objectives of this paper are: (1) to derive multitemporal HR LAI maps based on Sentinel-2A real data acquired during an entire European rice season for supporting agro-practices; (2) to validate the Sentinel-2A LAI estimates by direct comparison with ground LAI measurements conducted in three countries (Italy, Spain and Greece) during the season; and (3) to assess the feasibility of jointly using Sentinel-2A and Landsat-7/8 to assess a decametric high resolution multitemporal time series with an unprecedented frequency.It is worth mentioning that in this study, multitemporal synthetic aperture radar (SAR) data (Sentinel-1A) were used to obtain seasonally-updated rice maps used for focusing optical-based LAI retrievals on rice areas.Sentinel-2A and Landsat-7/8 data were processed in near real time in the three countries, and the corresponding LAI maps were provided through the web-based geo-portal in the project framework to the users and modelers.
The next sections describe the study areas and the remote sensing data used in this study (Section 2), also describing Sentinel-1A rice mapping and Sentinel-2A LAI retrieval methodology (Section 3).Section 4 validates the LAI estimates through direct comparison with in situ measurements and discusses the spatial-temporal consistency between Sentinel-2A and Landsat LAI estimates as a proof of concept for the operational generation of multisensor decametric LAI time series.Finally, Section 5 highlights the main conclusions.

Study Areas
The study areas have been selected in three Mediterranean countries, which are responsible for 85% of total European rice production: Italy (51.9%),Spain (25.4%) and Greece (7.0%) (see Figure 1).The Spanish study area is located in the rice district of Valencia (east of Spain) belonging to the Albufera Natural Park, which is a special protection area, thus allowing only rice crop practices.The Italian study area belongs to the Lomellina rice district (south-western Lombardy region), where rice is the dominant crop (>90%).The Greek study area is located in the rice district of Thessaloniki, which is the main cultivation area for Greece.Within each study area, rice is a common crop with a long tradition and economic value.

Field Campaigns
In the framework of the local ERMES field activities, LAI ground measurements were conducted over the study areas.In Spain, Italy and Greece, 32, 16 and 10 ESUs (elementary sampling units) were selected.The temporal frequency of the campaigns was approximately 7-10 days starting from the very beginning of rice emergence (early June) up to the maximum green rice LAI development (mid-August).The field campaigns were planned in such a way as LAI measurements within every ESU were either temporally coincident or 1-2 days apart with respect to the Sentinel-2A images.The same sampling scheme was used over each ESU, following the guidelines and recommendations of the Validation of Land European Remote sensing Instruments (VALERI) protocol (http://w3.avignon.inra.fr/valeri/).In order to characterize the spatial variability within each ESU, a range of 18-24 measurements was taken allowing one to obtain a statistically significant mean LAI estimate per ESU.The center of the ESU was geo-located using a GPS for later matching and associating the mean LAI estimate with the corresponding satellite estimation.
LAI measurements were acquired using a dedicated smartphone app (PocketLAI), which uses both smartphone's accelerometer and camera to acquire images at 57.5 • below the canopy and compute LAI estimates through an internal segmentation algorithm [31].PocketLAI estimates can reproduce destructive LAI measurements with acceptable results in terms of both reliability and accuracy [31].PocketLAI allows the acquisition of in situ LAI measurements at an affordable cost both in computational and human resources and aligned well with estimates obtained using plant canopy analyzers and DHP (digital hemispherical photography) techniques [32].

Sentinel-2A
The Sentinel-2A mission of the Copernicus program provides enhanced continuity of data so far provided by SPOT-4/5 and Landsat-7/8 missions.In this framework, the European Space Agency (ESA) provides free and open access to the Sentinels' data from the Sentinels Scientific Data Hub (https://scihub.copernicus.eu/).In this study, Sentinel-2A Level 1C data (top-of-atmosphere reflectances in cartographic geometry) were downloaded from early May up to the end of September in 2016 for the three study areas providing information every 10 days in 13 bands in the visible, near infra-red and short wave infra-red spectrum at a 10-, 20-and 60-m spatial resolution depending on the spectral band.Table S1 in the Supplementary Materials section lists the Sentinel-2A images used over the three areas.

Landsat-7/8
The Landsat Data Continuity Mission (LDCM) provides free access to the Landsat-7/8 archive.In this study, Landsat-8 Operational Land Imager (OLI) and Landsat-7 Enhanced Thematic Mapper (ETM+) surface reflectance data at 30-m spatial resolution were downloaded through the United States Geological Survey (USGS) Earth Resources Observation and Science (EROS) Center Science Processing Architecture (ESPA) during the 2016 rice season over the three study areas.The provisional Landsat-8 Surface Reflectance (LaSRC) [33] and the Landsat-7 ETM+ LEDAPS (Landsat Ecosystem Disturbance Adaptive Processing System) products were used as inputs to derive Landsat--7/8 LAI estimates.Images were available every 16 days in Italy and Greece.On the other hand, since the Spanish rice area lies in two Landsat paths within the same row, the temporal resolution of the images increased up to seven and nine days.See Table S1 in the Supplementary Materials section for further details.

Sentinel-1A
The Sentinel-1A SAR sensor operates at 5.405 GHz (C-band, corresponding to a microwave wavelength of about 5.6 cm) with VH and VV polarizations offering 12 days of revisit time.Sentinel-1A data over the three study areas were downloaded from the ESA Sentinels Scientific Data Hub during the 2016 European rice season.A fully-automated processing chain was followed in order to convert the multi-temporal SAR single look complex (SLC) data into terrain geocoded σ 0 values.The processing chain is a module within MAPscape-RICE, which included the following steps: strip mosaicking when SAR data are zero-Doppler focused, image co-registration, time-series speckle filtering, terrain geocoding, radiometric calibration and normalization, anisotropic non-linear diffusion (ANLD) filtering and the removal of atmospheric attenuation.The ANLD filter performs a smoothing and a high level of regularization in homogenous areas, while preserving inhomogeneous signal variations (discontinuities), which are typical of features, such as field edges, roads, rivers and other image information that we want to preserve [34].On the other hand, although microwave signals are able to penetrate clouds, it is possible to observe sudden changes in the dielectric constant at shorter wavelengths (X-and C-band) due to atmospheric moisture changes, dew, rain and severe storms.These events affect the temporal signature of σ 0 in two ways: (i) the thick layer of water vapor generates a strong decrease in σ 0 during the event, followed by a strong increase after the event; (ii) the intense rainfall generates a strong increase in σ 0 during the event, followed by a strong decrease after the event.These effects were removed by identifying anomalous peaks or troughs and correcting the σ 0 values by means of an interpolator.

Retrieval Methodology
The PROSAIL radiative transfer model was ran in forward mode for building a database, which was used for training the retrieval model and for mimicking canopy reflectance using the turbid medium assumption, which is particularly well suited for homogeneous canopies like rice [18,35].It assumes the canopy as a turbid medium for which leaves are randomly distributed.PROSAIL simulates leaf reflectance for the optical spectrum from 400-2500 nm with a 1-nm spectral resolution and as a function of biochemistry and structure of the canopy, its leaves, the background soil reflectance and the sun-view geometry.The mesophyll structural parameter (N), leaf chlorophyll (C ab ), dry matter (C m ) and water (C w ) contents are the parameters that determine the leaf optical properties.C w was tied to the dry matter content (C w = C m × C wREL /(1 − C wREL )) assuming that green leaves have a relative water content (C wREL ) varying within a relatively small range [36].The canopy structure is characterized by the average leaf angle inclination (ALA), the LAI and the hot-spot parameter (hotspot).A multiplicative brightness parameter (β s ) was introduced and applied to spectral flooded and dry soil signatures to represent different background reflectance types [36,37].The system geometry was described by the solar zenith angle (θ s ), view zenith angle (θ v ) and the relative azimuth angle between both angles (∆Θ).
Sub-pixel rice mixed conditions [18] were represented by a linear mixture of vegetation (vCover) and bare/flooded soil (1-vCover) spectra.In this study, vCover was assumed to follow a truncated Gaussian distribution (see Table 1).The leaf and canopy variables, as well as the soil brightness parameter were randomly generated following the parametrization in [18] in order to constrain the behavior of the model to Mediterranean rice areas (see Table 1).In particular, a spectral library of underlying rice background (flooded and dry) signatures was used to obtain multitemporal LAI retrievals robust against changes in background condition related to water management.The PROSAIL was run in forward mode to build a database of 2000 pairs of multispectral data and LAI.
The inversion of PROSAIL was done using standard regression.In particular, we used a Gaussian process regression (GPR) model [29].GPR has been recently introduced as a powerful regression tool in remote sensing studies [30].The GPR model exploits the relations between the input (e.g., bands spectra) x = [x 1 , . . ., x B ] ∈ R B and the output variable (i.e., LAI) y ∈ R in a training phase, and the predictive mean function for a new test point x * is expressed in the form: where {x i } N i=1 are the N spectra used in the training phase, α i is the weight assigned to each one of them, α o is the bias in the regression function and k θ is a kernel or covariance function (parametrized by a set of hyperparameters θ) that evaluates the similarity between the test spectrum and all N training spectra.In order to generate a kernel regression model, one needs to specify a covariance/kernel function k θ , to infer its hyperparameters θ and model weights α = [α o , α 1 , . . ., α N ].For the GPR prediction model, we used the so-called automatic relevance determination (ARD) kernel function: where ν is a scaling factor, The inversion process was conducted using 70% of the samples for training and the rest for testing with our SimpleR MATLAB toolbox, freely available at the Image Processing Laboratory website http://isp.uv.es/.The toolbox is intended for practitioners with little expertise in machine learning and that may want to assess advanced methods in their problems easily.

Sentinel-2A Processing
In the case of Sentinel-2A, top of canopy reflectances were obtained from Level 1C data after the atmospheric correction (AC) made using the Sen2Cor (Sentinel-2 atmospheric Correction) toolbox [38].Sen2cor produces an orthoimage of top of canopy corrected reflectance for all bands, excluding the cirrus band, which does not contain surface information.
The atmospheric correction performed by Sen2cor is based on the inversion of the radiative transfer equation through the use of a large set of look-up tables (LUT) of sensor-specific functions, already integrated in the processor [39].The LUT has been generated using the libRadtran, which contains a wide variety of atmospheric conditions, solar geometries and ground elevations.This database is generated with a high spectral resolution (0.6 nm) and then resampled using Sentinel-2A spectral responses.The Sen2cor algorithm allows the detection of clouds, water, snow and cloud shadows.The algorithm is based on a series of threshold tests over the Sentinel-2A Level 1C data, band ratios and spectral indexes, which produces a probabilistic cloud mask and snow mask quality indicators [40].
Six Sentinel-2A surface reflectance spectral bands were used during the retrieval process: blue (Band 2), green (Band 3), red (Band 4), near infrared (Band 8) and the two short wave infrared channels (Bands 11 and 12).These channels were selected to enhance the consistency with Landsat-7/8 data [18], allowing thus the creation of a robust multi-sensor retrieval.We have disregarded red-edge bands, despite their potential usefulness for the retrieval of important biophysical parameters, such as chlorophyll, nitrogen content and brown LAI [21,41].
The combined use of Sentinel-2A and Landsat-7/8 poses a number of technical challenges, due to the differences in their orbital, spatial, spectral response functions and image processing chains.In fact, although the radiometric characteristics of these sensors are similar, they are not identical and can lead to slight differences in surface reflectance and retrieval quantities [42], as confirmed by the inspection of the data.In order to reduce this bias in the multi-sensor dataset, the Sentinel-2A images were radiometrically normalized by assuming a linear relationship between the reflectance of Sentinel-2A and Landsat-8 images.For each band, slope and intercept parameters were obtained using orthogonal linear regression.We applied Iteratively Reweighted Multivariate Alteration Detection (IRMAD) transformation [43] to select invariant pixels, which proved to be invariant to changes in sensor gain and offset and radiometric and atmospheric corrections.

Sentinel-1A Rice Maps
Rice maps derived from Sentinel-1A data were used as masking layer for LAI retrieval.The multi-temporal rule-based rice detection (MSRD) algorithm proposed by [44] for X-band HH SAR data was adapted to C-band VV/VH data and applied to the Sentinel-1A datasets over the three study areas.Processing has been carried out in MAPscape-RICE software.The MRSD algorithm relies on rules applied to the temporal profile of Sentinel-1A σ 0 and defined based on a priori knowledge on the rice calendar, crop practices and agro-ecological conditions of the study areas.The following rules are implemented for each pixel: 1.
Rice exclusion condition: average σ 0 lower/higher than expected or average σ 0 below a minimum value for longer than expected or variation of σ 0 larger than expected; 2.
Presence of flooding conditions at the start of season (SoS): the temporal series is analyzed starting from the first image acquisition to identify when σ 0 drops below a maximum values for SoS flooding; this time is set as SoS; 3.
Confirmation of rice presence: after flooding detection; rice presence is confirmed if σ 0 increase after SoS is consistent with expected value for rice crops (rapid growth of rice biomass after flooding); 4.
Late rice condition: when the length of the S1A time series after detected SoS is shorter than expected, rice is classified as being sown late in the season; 5.
New rice/crop season: unexpected drop in σ 0 after SoS, which suggest a new flood or a new crop season.Notice that if Rule 2 is not met during the entire temporal series, another rule allows one to determine whether flooding occurred before the first Sentinel-1A acquisition.If σ 0 variation is above a pre-defined value, which is expected for rice and σ 0 decreases after the maximum value, the pixel is classified as early rice.
Crop maps for Italy were produced for the entire Rice district of Piedmont and Lombardy (∼210,000 ha).The Lomellina test site (∼50,000 ha) where LAI was retrieved belongs to this area.

Rice Detection from Sentinel-1A
2016 rice maps were produced over the three study areas: Greece (GR), Italy (IT), Spain (SP) from Sentinel-1A (see Figure 2).Validation has been performed by comparison of the rice mapping products with reference observations collected during the season over the three study areas.In particular, in Italy, point data were used as provided by field survey exploiting a Voluntary Geographic Information App (S4A Smart App) [45], while for Greece and Spain, polygons were digitalized over HR optical data based on ground inspection during field campaigns.A confusion matrix [46] was performed at the pixel level (IT: 1833 rice and 1394 no rice pixels; GR: 8181 rice and 7507 no rice pixels; ES: 21792 rice and 7376 no rice pixels), and accuracy measures (overall accuracy and Kappa coefficient (K)) were derived.Results show that the maps have an OA and K greater than 89 and 0.78 respectively in all cases (IT: OA = 89.1,K = 0.78, GR: OA = 96.7,K = 0.93 and ES: OA = 99.3,K = 0.99).The corresponding confusion matrices are shown in Table S2 of the Supplementary Materials section.In general, the obtained results are more than satisfactory, in particular in Spain and Greece, where an overall accuracy higher than 95% has been achieved.In Italy, a lower overall accuracy, but still more than satisfactory, is reported.The reason is due to the different rice practices (wet and dry sowing) and the presence of other crops (in particular maize and soya), while showing, at the beginning of the crop season, a similar backscattering coefficient.

Sentinel-2A LAI Estimates
High-resolution (10 m) LAI retrievals were obtained over the three rice areas during the 2016 rice season applying the trained model to Sentinel-2A data.Figure 3 shows the LAI map derived over the Spanish study area using a Sentinel-2A image acquired on 9 August 2016.The high spatial detail of estimates can be seen, which allows identifying significant different values within the same rice field.Those intra-field LAI differences are mainly due to the heterogeneity of the rice field caused by non-homogenous seeding and agro-practices [18,47].This kind of information provided by LAI estimates at such a resolution is indeed very useful for rice crop users and modelers.In fact, key farmers were able to relate these local maps with their field water flow regime and cultivation practices.From interactions with farmers, it becomes clear that LAI variability maps at dense time points provide relevant information on crop failure of immediate use for improving their agro-practices.Note that non-rice areas were masked out during the retrieval process using the rice maps based on Sentinel-1A data acquired over each study area.These LAI maps can be used to continuously monitor the on-going cropping season and for early detection of crop growth anomaly connected to potential rice damage and yield loss.Preliminary effective application was conducted during the 2016 cropping season in the framework of the ERMES project.LAI maps were produced in near real time and provided through the ERMES geo-portal to farmers and users.In Italy, this information resulted in being useful to identify and assess the impact of the misuse of an herbicide.The maps, computed after the event and compared with the previous ones, were used together with the farmer to clearly delineate which part of the field was mainly impacted and needed recovery.On the other hand, anomalies in Spanish LAI time series were identified and related to fields affected by a rice disease thanks to the interaction with users.These results revealed that operational production of decametric LAI time series can represent a useful screening tools to support field inspection to identify situations likely affected by a crop problem that if not managed will produce yield shortage.This kind of information is also critical in crop modeling, since the identification of anomalous drops in LAI values can improve the accuracy of yield forecasting.
The availability of LAI ground measurements allowed assessing the high-resolution LAI estimates.The root mean squared error (RMSE), mean error (ME), mean absolute error (MAE) and coefficient of determination (R 2 ) were computed in order to assess the accuracy of the retrievals, bias and goodness-of-fit (see the values in Figure 4 and the metrics in Table S3 of the Supplementary Materials section).Good accuracy and high correlation were found revealing an overall RMSE of 0.69 m 2 /m 2 and a coefficient of determination of 0.95.A slight bias was observed between Sentinel-2A LAI estimates and ground data.More in detail, Table 2 shows the statistical indicators between ground data and Sentinel-2A estimates over each study area.In general, very high correlations were found.Low biases were observed in the case of the Spanish study area, while in Italy and Greece, small biases were observed with respect to the ground data (see Table 2).Biases could be explained due the rice ESU heterogeneity and the slight underestimation of PocketLAI regarding other measurements [32].Relative RMSE to mean (rRMSE m ) and relative RMS to range (rRMSE r ) are also reported in Table 2 and Figure 4 showing good accuracies and revealing that the Sentinel-2A LAI estimates are suitable to be integrated into crop models [48].Table 2. Statistical indicators (coefficient of determination (R 2 ), root mean squared error (RMSE), mean error (ME), mean absolute error (MAE), relative root mean squared error to mean (rRMSE m = RMSE/ ŷi )) and relative root mean squared error to range (rRMSE r = RMSE/(max(y i ) − min(y i ))) between the ground measurements and LAI estimates during the 2016 rice season over the study areas.

Spain
Italy Greece

Spatial-Temporal Consistency between Sentinel-2A and Landsat LAI Estimates
Sentinel-2A LAI estimates were spatially assessed with estimates obtained following the approach described in [18] for Landsat-8 LAI retrieval.Spatial consistency between Sentinel-2A and Landsat-8 LAI estimates was assessed computing the difference map of the pair of the closest acquisition dates between Landsat-8 and Sentinel-2A over each study area (see Figure 5).For this purposes, the Sentinel-2A images were resampled to Landsat-8 resolution, and the comparison was achieved averaging the LAI value computed over a window of 3 × 3 valid pixels in order to reduce coregistration errors between images [49] and inconsistencies associated with differences in the point spread functions [50].Good correlations were found in all three study areas, and a notable agreement was found in Italy with practically no bias and low RMSE (see Figure 5).Leaf area index difference maps show the predominance of the light green color, which demonstrates that the retrievals from the two sensors are coherent, exhibiting LAI differences around zero over the majority of the rice areas.Higher differences are mainly due to differences in surface reflectances from Landsat-8 and Sentinel-2A generated by the different atmospheric correction methods.In the case of Spain (Figure 5, top-left), a surrounding halo effect can be discerned.This effect is due a suboptimal Landsat-8 atmospheric correction performance over the edge region between land and sea and could explain the wider dispersion observed in the scatter plot (see Figure 5, bottom-left).In addition to the good spatial consistency between Sentinel-2A and Landsat-8 estimates, the 2016 LAI temporal behaviors are shown in Figure 6 which exhibits the temporal consistency of LAI estimates from Sentinel-2A (blue dots) with the ones from Landsat-7/8 (orange/green dots).Each dot is the estimated LAI averaged over a rice field; over each study site.The fields shown in Figure 6 were selected to cover LAI variability as large as possible and, where available, cultivated with local traditional varieties (Bomba and Carnaroli in the top panels of Spain and Italy).All of the other panels show modern rice varieties (Olympiada and Mare CL in Greece, Sirio CL and Selenio in the bottom panels of Spain and Italy, respectively).In Italy, there is larger variability of agro-practices, which is highlighted in the example field of Figure 6 by the significant difference in sowing dates (vertical dashed bars).Moreover, according to in situ observations, Field #68 is characterized by the presence of a cover crop preceding rice sowing, which is also confirmed by LAI profiles (high LAI values before the abrupt drop due to field preparation (ploughing and laser leveling) for rice sowing).
Cloud contamination is a common problem, which limits the utility of passive optical multispectral images.The use of multi-sensor approach allows filling time gaps typical of single sensor time series, due primarily to cloud cover, thus improving the usefulness of the satellite estimation.Figure 6 demonstrates the temporal consistency of multi-source LAI estimates as provided by Sentinel-2A and Landsat-7/8 images, since neither bias nor shifts can be clearly recognized.For example, in Field #36, missing Landsat-7/8 observations at the peak of the season due to cloud cover are filled by Sentinel-2A, thus improving retrieval accuracy of the maximum LAI value at the peak of the season.In'detail, the merged time series reaches a maximum LAI value of 5.3 m 2 /m 2 at DOY (day of year) 203, while the Landsat-7/8 series has a maximum LAI value of 4.3 m 2 /m 2 at DOY 215.As in this example, a multi-sensor approach provides more accurate estimates of LAI at key phenological stages related to the shape of the curve given by the temporal profile (e.g., peaks, inflection points).Similar remarks can be done for Field #8267 in Greece, where three Sentinel-2A images are available at the peak of the season, but only one from Landsat-7/8.The Spanish study site does not suffer a lack of Landsat-7/8 images because it is located in the overlapping region between two Landsat paths, which enhances the frequency of observation up to seven and nine days and the likelihood of clear sky conditions at image acquisition time.These preliminary results suggest that a very frequent time series of LAI at high resolution can be obtained from a multi-sensor approach to better outline rice-growing behavior.The combined curves of LAI could be exploited for: (i) the retrieval of phenological stages for crop modeling purposes [51][52][53]; (ii) deriving multitemporal training sets for mapping purposes [54]; or (iii) monitoring vegetation production by computing the area under the curve [55].

On the Sources of Error and Uncertainty
It is worth mentioning that indirect and non-destructive methods, such as PocketLAI, provide estimates associated with several sources of measurement error, including wrong measurements and suboptimal sampling within an ESU, simplification of canopy leaf optical properties, illumination conditions and the optical signal saturation in dense vegetation.Specifically, PocketLAI measurements over rice crops showed an underestimation with respect to classical instrumentation, such as DHP [32].The PocketLAI underestimation can be due to the fact that measurement is based on upward photography, which underestimates LAI in rice crops [56], while DHP techniques over rice typically use downward-looking photography due to the characteristics of rice crops (low height and flooded background).In addition, LAI measurements over rice present a saturation during flowering, fruit development and ripening phenological rice stages.Moreover, the specific background conditions of rice crops, which are flooded most of the time, except on some dates, when water is pumped out for agronomic purposes, may affect the temporal trend of non-destructive LAI ground measurements.
In this study, a specific PROSAIL parametrization for simulating rice cropping areas was used obtaining thus a single model, which was applied over the three countries.The sensitivity of the model to each region can be observed in Table 2, where estimates were more accurate and less skewed in Spain; nevertheless, the results for Italy and Greece were also good, revealing the robustness of the model.On the other hand, The GPR model not only provides a mean estimate for the predictions, but also a prediction uncertainty related to the confidence of the estimates (the lower the uncertainty, the higher the confidence).Although this uncertainty is a numerical value, it must be interpreted as a quality-driven indicator [22,57] related, in the case of rice, to (i) an incorrect simulation of rice background, which has a high impact on the reliability of LAI retrievals in rice crops, and (ii) to spectral data under-represented in the PROSAIL training database, such as water bodies and man-made surfaces, for which its associated LAI estimates presented significantly higher prediction uncertainty (lower confidence).Uncertainty values over the three areas remained constant at a value of about 0.8 during the rice growing period, revealing the robustness of the GPR model.In addition, unexpected estimates (pronounced drops) provided by the GPR were also observed during the season and related to suboptimal atmospheric correction or undetected clouds in Sentinel-2A data, which were labeled as poor-quality estimates and consequently excluded.This feature allows one to create a quality flag indicator useful for users and crop modelers.Figure 7 shows the temporal evolution of Sentinel-2A LAI and uncertainty within a rice field.Note that the unexpected LAI drop reported on 7 September the associated σ peak were related to a cloudy pixel undetected by Sen2cor.

Conclusions
This study presents 2016 seasonal LAI estimation from Sentinel-2A imagery and validation over the major European rice areas.The approach relies on the inversion of the PROSAIL radiative transfer model with Gaussian process regression on rice field automatically detected using Sentinel-1A data.The proposed approach allows us to retrieve multitemporal high-resolution LAI estimates by exploiting the enhanced spatial and spectral resolutions of the ESA Sentinels constellation.
Sentinel-2A estimates are in agreement with field measurements acquired in the three rice areas during the 2016 crop season.In situ measurements were acquired with smartphones using the PocketLAI app, which allowed conducting coordinated field campaigns from the beginning of the leaf development up to the flowering phase.The spatial and multi-temporal consistency of Sentinel-2A LAI estimates was assessed by comparison with LAI maps obtained from a consolidated and published approach applied to Landsat-7/8 images [18].The single date pixel-by-pixel agreement shows the goodness of fit of Sentinel-2A and Landsat-8 LAI estimations over the three rice areas (R 2 ranging from 0.83-0.94).In addition, the uncertainty provided by the GPR allows enhancing the detection of artifacts, such as undetected clouds by the Sen2cor module.
Preliminary tests conducted by combining Sentinel-2A and Landsat-7/8 LAI time series highlighted the feasibility of deriving high temporal resolution and multi-sensor LAI estimates.The increased temporal resolution of multi-sensor LAI allows filling gaps mainly due to atmospheric interference to obtain more reliable time series for precision agriculture applications and rice monitoring.A dense temporal dataset of LAI maps is in fact fundamental to perform expert crop monitoring and/or improve crop model estimations exploiting assimilation techniques.The free
n accounts for the noise standard deviation and σ b is a dedicated parameter controlling the spread of the relations for each particular spectral band b (b = 1, . . ., B). Model hyperparameters are collectively grouped into θ = [ν, σ n , σ 1 , . . ., σ B ], and model weights α can be automatically obtained by maximizing the marginal likelihood in the training set.The obtained weights α after optimization give the relevance of each datum x i , while the inverse of σ b represents the relevance of each band b.Hence, low values of σ b indicate a higher informative content of this certain band b to the training function k.

Figure 2 .
Figure 2. Rice area for Italy (IT) (left), Greece (GR) (middle) and Spain (ES) (right) at the end of the 2016 season and points/field boundaries over rice (red) and not rice (blue) fields.

Figure 3 .
Figure 3. High-resolution LAI map obtained with Sentinel-2A data acquired in the Spanish study area on 9 August 2016.The right panel provides a zoom over a monitored field in the rice season.

Figure 4 .
Figure 4. Scatter plots of Sentinel-2A estimated LAI values versus in situ LAI measurements during the 2016 rice season.For the sake of visualization, the standard deviation of measurements and Sentinel-2A prediction uncertainty (≈±0.8) are not shown.

Figure 5 .
Figure 5. LAI difference maps (Top) and the corresponding scatter plots of the closest Sentinel-2A and Landsat-8 LAI estimates (Bottom) over Spain (Left), Italy (Middle) and Greece (Right).

Figure 6 .
Figure 6.Example of merged series of LAI on three sites.Dotted vertical bars indicate the sowing date.Numbers indicate fields ID.For the sake of visualization, standard deviations are not shown.

Figure 7 .
Figure 7. Temporal evolution of LAI and uncertainty (σ) over a rice field in Spain during the 2016 rice season.

Table 1 .
Distribution of the parameters within the PROSAIL RTM at canopy (leaf area index (LAI), leaf angle inclination (ALA), hot-spot parameter (hotspot) and vegetation cover (vCover)) and leaf (mesophyll structural parameter (N), leaf chlorophyll (C ab ), dry matter (C m ) and relative water (C w ) contents) levels and the soil brightness parameter (β s ).