remote sensing Article Retrieval of High Spatiotemporal Resolution Leaf Area Index with Gaussian Processes, Wireless Sensor Network, and Satellite Data Fusion Gaofei Yin 1, * , Aleixandre Verger 2 , Yonghua Qu 3 , Wei Zhao 4 , Baodong Xu 5 , Yelu Zeng 6 , Ke Liu 7 , Jing Li 8, * and Qinhuo Liu 8 1 Faculty of Geosciences and Environmental Engineering, Southwest Jiaotong University, Chengdu 610031, China 2 CREAF, 08193 Cerdanyola del Vallès, Catalonia, Spain; verger@creaf.uab.cat 3 State Key Laboratory of Remote Sensing Science, Beijing Key Laboratory for Remote Sensing of Environment and Digital Cities, Institute of Remote Sensing Science and Engineering, Faculty of Geography Science, Beijing Normal University, Beijing 100875, China; qyh@bnu.edu.cn 4 Research Center for Digital Mountain and Remote Sensing Application, Institute of Mountain Hazards and Environment, Chinese Academy Sciences, Chengdu 610010, China; zhaow@imde.ac.cn 5 Macro Agriculture Research Institute, College of Resource and Environment, Huazhong Agricultural University, Wuhan 430070, China; xubaodong@mail.hzau.edu.cn 6 Department of Global Ecology, Carnegie Institution for Science, Stanford, CA 94305, USA; yzeng@carnegiescience.edu 7 Institute of Remote Sensing Application, Sichuan Academy of Agricultural Science, Chengdu 610066, China; billc_st@163.com 8 State Key Laboratory of Remote Sensing Science, Institute of Remote Sensing and Digital Earth, Chinese Academy of Sciences, Beijing 100101, China; liuqh@radi.ac.cn * Correspondence: yingf@swjtu.edu.cn (G.Y.); lijing01@radi.ac.cn (J.L.)  Received: 21 December 2018; Accepted: 21 January 2019; Published: 24 January 2019  Abstract: Many applications, including crop growth and yield monitoring, require accurate long-term time series of leaf area index (LAI) at high spatiotemporal resolution with a quantification of the associated uncertainties. We propose an LAI retrieval approach based on a combination of the LAINet observation system, the Consistent Adjustment of the Climatology to Actual Observations (CACAO) method, and Gaussian process regression (GPR). First, the LAINet wireless sensor network provides temporally continuous field measurements of LAI. Then, the CACAO approach generates synchronous reflectance data at high spatiotemporal resolution (30-m and 8-day) from the fusion of multitemporal MODIS and high spatial resolution Landsat satellite imagery. Finally, the GPR machine learning regression algorithm retrieves the LAI maps and their associated uncertainties. A case study in a cropland site in China showed that the accuracy of LAI retrievals is 0.36 (12.7%) in terms of root mean square error and R2 = 0.88 correlation with ground measurements as evaluated over the entire growing season. This paper demonstrates the potential of the joint use of newly developed software and hardware technologies in deriving concomitant LAI and uncertainty maps with high spatiotemporal resolution. It will contribute to precision agriculture, as well as to the retrieval and validation of LAI products. Keywords: leaf area index; uncertainty; Gaussian processes; wireless sensor network; data fusion; Landsat; MODIS; validation Remote Sens. 2019, 11, 244; doi:10.3390/rs11030244 www.mdpi.com/journal/remotesensing Remote Sens. 2019, 11, 244 2 of 18 1. Introduction The leaf area index (LAI) describes the surface area available per unit ground surface for energy and mass exchanges between vegetation and the atmosphere. It plays a key role in several surface processes, including photosynthesis, respiration, and evapotranspiration. Spatiotemporal continuous LAI data is critical for climate model validation [1], regional to global scale carbon budget estimation [2], disturbance detection [3], and many other applications. Remote sensing is the only tool that can effectively generate spatiotemporal continuous LAI products [4,5]. The methods to estimate LAI values from satellite data can be divided into empirical methods based on statistical relationships with ground measurements [6–8], physical methods based on the inversion of radiative transfer (RT) models [9–11], and hybrid methods that combine statistical methods with RT models [12–14]. The physical and hybrid methods are increasingly used in the scientific community because of their strong generalization ability [15,16]. However, the inversion of a canopy RT model is an under-determined problem because of the limited information content of the radiometric signal compared to the high number of unknowns influencing the canopy reflectance [17]. The ill-posed nature of the inversion problem—several sets of input variables can yield very similar spectra [18,19]—hampers the application of the physical and hybrid methods [18], so empirical methods are still popular, especially in regional-scale studies. In situ LAI measurements, required in the calibration of the empirical methods, are usually derived from optical instruments [20–22], including digital hemispheric photographs [23,24], TRAC [25], LAI-2000 Plant Canopy Analyzer [26], or its successor LAI-2200. The manual nature of the optical instruments makes the traditional in situ measurements time-consuming and labor-intensive. The emergence of more portable or even fully automatic observation instruments, such as LED sensors [27] and smartphones [28,29], represents a change in paradigm in near-surface remote sensing. Several operational LAI products, including MODIS [30], GLASS [31], and Copernicus Global Land [32,33] products, are regularly provided by space agencies over long periods at high temporal resolution (between 8- and 15-day frequency) but at coarse spatial resolution (between 300 m and 8 km). On the other hand, several studies have generated LAI maps with high spatial resolution but with a limited temporal sampling [34–36]. In fact, the technical constraints of remote sensing instruments [37] involves a tradeoff between spatial and temporal resolution for remote sensing application. This hinders the utility of LAI products in areas with high spatial heterogeneity or substantial temporal variations [38–40]. LAI time series with high spatiotemporal resolution are urgently needed for characterizing terrestrial ecosystem processes at the parcel level. In a previous work [8], we proposed a framework for the retrieval of high spatiotemporal LAI products based on the combined use of wireless sensor network (WSN) and data blending (DB) technologies. The LAINet observation system [41,42] was used to estimate temporally continuous field measurements at high temporal frequency. The Consistent Adjustment of the Climatology to Actual Observations (CACAO) [43] was used as a DB technique to produce fine spatial resolution satellite imagery. Then, an exponential function relating the reconstructed remote sensing observations to the LAINet measurements was calibrated and used to predict LAI maps with high spatiotemporal resolution. However, this approach has two main limitations: (i) accuracy of LAI retrievals not meeting the user requirements (15%) proposed by the Global Climate Observing System [44] and (ii) non-retrieval of the uncertainty associated with the LAI maps. Characterizing the uncertainties associated with LAI maps is required in many applications including the assimilation into land surface models [45,46]. Some of the existing LAI products provide uncertainty estimates. In the MODIS LAI products, the uncertainty is the standard deviation over all acceptable solutions of the inversion process [4]. In other cases, such as Copernicus Global Land 1 km Version 1 LAI products [32], the uncertainty is estimated using a dedicated algorithm that is independent of the main LAI retrieval algorithm. The rapid development of machine learning regression methods, particularly those based on Bayesian approaches, sheds light on the simultaneous estimation of LAI and uncertainty in a single algorithm [47]. Gaussian process regression (GPR) [48,49] Remote Sens. 2019, 11, 244 3 of 18 has been Remote Sens.reported to improve 2019, 11, x FOR accuracy PEER REVIEW in LAI retrieval and provide uncertainty estimates directly 3 of 18 through Gaussian probabilities [7]. We aim to We aim to develop develop aa method method toto generate generate more more accurate accurate LAI LAI time time series series and and their their associated associated uncertainties at high spatiotemporal resolution. This will contribute to the temporally continuous uncertainties at high spatiotemporal resolution. This will contribute to the temporally continuous validation validation of of coarse coarse spatial spatial resolution resolution LAI LAI products products and and the assimilation of the assimilation LAI estimates of LAI estimates into into land land surface surface models, for which high spatiotemporal resolution maps of both LAI estimations and the models, for which high spatiotemporal resolution maps of both LAI estimations and the corresponding uncertainties are corresponding uncertainties are needed. needed. 2. Material and Methods 2.1. Study Area research was The research wasconducted km× 5×km conductedinina a5 5km 5 km region region in Huailai, in Huailai, northern northern ChinaChina (Figure (Figure 1). 1). The The annual average precipitation ◦ C, respectively. Croplands are annual average precipitation and and temperature temperature are mm are 396 396 mm and and 9.6 respectively. 9.6 °C, Croplands are the the main main landland cover cover type, type, andand thethe growing growing season season typically typically spansfrom spans frommid-May mid-Maytotolate lateSeptember September[50]. [50]. The remaining remaininglandlandcover covertypes typesinclude water include in in water thethe northwest andand northwest impervious surfaces impervious scattered surfaces over scattered the study area. over the study area. Figure 1. Map of thethe study studyarea. area.The Thebackground backgroundimage imageisisthe thefalse false color color composition composition with with bands bands 6- 6-5-4 5-4 of the Landsat 8 Operational Land Imager (OLI) image acquired on August 23, 2013. The points of the Landsat 8 Operational Land Imager (OLI) image acquired on August 23, 2013. The points and numbers and indicate numbers the locations indicate of the twelve the locations below nodes of the twelve belowinnodes the LAINet in the observation system deployed LAINet observation system in the study deployed area.study area. in the 2.2. Satellite 2.2. Satellite Data Data 2.2.1. Landsat 8 OLI 2.2.1 Landsat 8 OLI We used Landsat 8 Operational Land Imager (OLI) data acquired through the Earth Resources We used Landsat 8 Operational Land Imager (OLI) data acquired through the Earth Resources Observation and Science Center Science Processing Architecture (ESPA) (https://espa.cr.usgs.gov). Observation and Science Center Science Processing Architecture (ESPA) (https://espa.cr.usgs.gov). ESPA provides top of canopy reflectances based on the Second Simulation of the Satellite Signal in ESPA provides top of canopy reflectances based on the Second Simulation of the Satellite Signal in the Solar Spectrum Vectorial (6SV) atmospheric correction algorithm [51]. Spectral bands in the blue, green, the Solar Spectrum Vectorial (6SV) atmospheric correction algorithm [51]. Spectral bands in the blue, red, near infrared (NIR), and the two shortwave infrared (SWIR1 and SWIR2) bands were used (Table 1). green, red, near infrared (NIR), and the two shortwave infrared (SWIR1 and SWIR2) bands were used (Table 1). Table 1. Band numbers (BN), central wavelengths (CW), and full width at half maximum (FWHM) of the OLI and MODIS sensors for the selected spectral bands. OLI MODIS Spectral band CW FWHM CW FWHM BN BN (nm) (nm) (nm) (nm) Blue 2 482.0 60.1 3 469.0 11.5 Green 3 561.4 60.1 4 555.0 12.9 Remote Sens. 2019, 11, 244 4 of 18 Table 1. Band numbers (BN), central wavelengths (CW), and full width at half maximum (FWHM) of the OLI and MODIS sensors for the selected spectral bands. Remote Sens. 2019, band 11, x FOR PEER REVIEW OLI MODIS 4 of 18 Spectral BN CW (nm) FWHM (nm) BN CW (nm) FWHM (nm) Red Blue 2 4 654.6 482.0 60.137.5 3 1 645.0 469.0 10.5 11.5 Green Near infrared 3 5 561.4 864.7 60.128.2 4 2 555.0 858.5 12.9 11.6 Red 4 654.6 37.5 1 645.0 10.5 Short wavelength Near infrared infrared 5 6 1608.9 864.7 28.220.4 2 6 1640.0 858.5 14.0 11.6 Short wavelength Short wavelength infrared infrared 6 7 1608.9 2200.7 20.484.7 6 7 1640.0 2130.0 14.0 15.8 Short wavelength infrared 7 2200.7 84.7 7 2130.0 15.8 The dataset includes six cloud-free Landsat 8 OLI scenes (path: 123; row: 32) acquired on the day The of year dataset (DOY) includes 139, six267, 187, 235, cloud-free 299, andLandsat 315 for8the OLI scenes year 2013(path: 123; (Figure 2).row: This32) acquired period on the day corresponds to of year (DOY) 139, 187, 235, 267, 299, and 315 for the year 2013 (Figure 2). This period the growing season of the study area (Figure 4). Note that only two Landsat 8 OLI scenes (acquired corresponds to theDOY on growing 187 season and 235)of are the within study areathe (Figure 4).acquisition period of Note that only two Landsat of ground 8 OLIJuly data (from scenes 1 to(acquired Septemberon DOY 187 and 235) are within the period of acquisition of ground data (from July 1 to September 14, 2013, as detailed in Section 2.3.1), that is, only two ground-based LAI maps could be generated 14, 2013, as detailed from in Section Landsat 2.3.1), data data without that is, only two ground-based LAI maps could be generated from Landsat blending. data without data blending. Figure 2. Figure Temporal distributions 2. Temporal distributions of of the the MODIS MODIS and and Landsat Landsat acquisitions acquisitions during during the the year year 2013 2013 over over the study area. The shaded area represents the period of acquisition of ground LAINet measurements. the study area. The shaded area represents the period of acquisition of ground LAINet measurements. DOY: day of year. DOY: day of year. 2.2.2. MCD43A4 2.2.2. MCD43A4 We also used the nadir directionally normalized reflectance MODIS Nadir Bidirectional We alsoDistribution Reflectance used the nadirFunctiondirectionally normalized (BRDF) Adjusted reflectance Reflectance MODIS product Nadir Bidirectional (MCD43A4) [52]. It is Reflectance Distribution Function (BRDF) Adjusted Reflectance product produced at 500 m spatial resolution every 8 days using a compositing window of 16 days (MCD43A4) [52]. with It is produced data from at 500Aqua both m spatial resolution and Terra everyAll satellites. 8 days using in the scenes a compositing 2013 covering window of 16 the study days area with (tile data number from both Aqua and Terra satellites. All the scenes in 2013 covering the study of h24v04) were downloaded from the Land Processes Distributed Active Archive Center (LP DAAC) area (tile number of h24v04) were downloaded from the Land Processes Distributed Active Archive (https://lpdaac.usgs.gov). The cloud-contaminated scenes were all filtered out, resulting in 36 clean Center (LP DAAC) (https://lpdaac.usgs.gov). scenes with a full coverageThe cloud-contaminated of the growing season ofscenes were the study all(see area filtered out, Figure 2).resulting in 36 in The six bands clean the scenes with a full coverage of the growing season of the study area (see MCD43A4 product analogous to the Landsat 8 OLI were selected (Table 1). The MCD43A4 data wereFigure 2). The six bands in the MCD43A4 resampled product to 30-m analogous spatial resolutionto the Landsat using 8 OLIneighbor the nearest were selected sampling (Table 1). The MCD43A4 data method. were resampled to 30-m spatial resolution using the nearest neighbor sampling method. 2.3. Retrieval Algorithm 2.3. Retrieval Algorithm Figure 3 shows the general outline of the proposed framework to generate concomitant LAI Figure 3 shows and uncertainty maps thewith general highoutline of the proposed spatiotemporal framework resolution. to generate It includes LAINet,concomitant CACAO, and LAIGPR and uncertainty methods: maps with high spatiotemporal resolution. It includes LAINet, CACAO, and GPR methods: (i) (i)LAINet provided LAINet the LAI provided fieldfield the LAI measurements measurements for calibrating and and for calibrating validating GPR;GPR; validating (ii) (ii) CACAO CACAO was used for blending the Landsat 8 OLI and MODIS MCD43A4and was used for blending the Landsat 8 OLI and MODIS MCD43A4 products generating products and reflectance generating data with reflectance datahigh withspatiotemporal resolution; high spatiotemporal and, and, resolution; (iii) (iii) GPRGPR waswas usedused for the for retrieval of LAI the retrieval and and of LAI uncertainty maps. uncertainty maps. uncertainty maps with high spatiotemporal resolution. It includes LAINet, CACAO, and GPR methods: (i) LAINet provided the LAI field measurements for calibrating and validating GPR; (ii) CACAO was used for blending the Landsat 8 OLI and MODIS MCD43A4 products and generating reflectance data with high spatiotemporal resolution; and, Remote Sens. 2019, 11, 244 5 of 18 (iii) GPR was used for the retrieval of LAI and uncertainty maps. Figure 3. Flowchart of the methodology combining LAINet, CACAO, and GPR for the generation of LAI and uncertainty maps with high spatiotemporal resolution. 2.3.1. LAINet Field Measurements Field measurements were collected through the LAINet observation system. LAINet comprises 12 below nodes which receive the transmitted solar radiation below the canopy (EB ) (see Figure 1 for their spatial distribution), and 1 above node which receives downward radiation above the canopy (EA ) [42]. The LAI estimation is based on gap fraction measurements calculated as EB /EA . For details regarding LAINet, please refer to Qu et al. (2014) and our previous work [8]. Ground measurements were collected from July 1 to September 14, 2013 (see Figure 2). To reduce the random errors, the estimated LAI values from each node of the LAINet observation system were averaged over an 8-day interval. The 8-day composite measurements are hereafter indicated by the DOY of the first day of the compositing period. Due to unexpected instrument failures and unsatisfactory weather conditions, some LAINet measurements were invalid. Only four nodes (numbered 3, 5, 8, and 12 in Figure 1) had valid LAI values for every date. A total of 67 measurements were finally available. 2.3.2. CACAO Reconstruction of High Spatiotemporal Satellite Data The CACAO [43] method was used to blend Landsat 8 OLI data of high spatial resolution and MCD43A4 data of high temporal resolution, and to reconstruct high spatiotemporal resolution reflectance data. CACAO is a temporal smoothing and gap filling method with temporal constraints [43]. It fits a phenology model to the actual observations based on a two-step process: (i) A phenology model was first constructed per pixel and per band based on the temporal evolution of MCD43A4 reflectance data. The time series of MCD43A4 reflectance were smoothed using Savitzky–Golay (SG) filtering [53] to generate the phenology model. SG smooths the time series by fitting a low-degree polynomial in the local temporal window. In this study, the degree of the smoothing polynomial was fixed at 2 and the half-width of the smoothing window at 16 days. (ii) The phenology model was then fitted to the actual Landsat observations. It was shifted and scaled in order to minimize its difference with the Landsat 8 OLI data in terms of the root mean square error (RMSE): s 1 n n i∑ arg min ( ROLI (ti ) − scaleR pm (ti + shi f t))2 , (1) scale,shi f t =1 where ROLI represents the Landsat 8 OLI reflectance, Rpm represents the reflectance of the phenology model, n is the number of available Landsat 8 OLI observations (in this study, n = 6), t is the time in days, and scale and shift are the parameters to be adjusted. The minimization of the RMSE was done pixel by pixel, and band by band. After obtaining the values of scale and shift for each pixel in each band, Landsat-like reflectance can be computed at any arbitrary date from the MODIS phenology model. Remote Sens. 2019, 11, 244 6 of 18 The CACAO fitting is illustrated in Figure 4. Generally, the Landsat 8 OLI reflectance and the MODIS phenology model (SG-filtered MCD43A4 reflectance) show a similar seasonality with some differences in terms of magnitude and phase. The reconstructed reflectance retains the seasonality of the phenology model and fits the Landsat 8 OLI reflectance values. The reconstructed reflectances have an 8-day temporal frequency and a 30-m spatial resolution. Their temporal and spatial characteristics inherit from MCD43A4 and Landsat 8 OLI data, respectively. The time Remote Sens.series spans 2019, 11, from x FOR PEERDOY 121 to 321, 2013, covering the entire growing season in our study REVIEW area. 6 of 18 Figure 4. Figure Illustration of 4. Illustration of the the CACAO CACAO fitting fitting of of the the MODIS MODIS phenology phenology model model toto actual actual Landsat Landsat observations for observations for reconstructing reconstructing high high spatiotemporal spatiotemporal reflectance reflectance data. data. Case Case of of reflectance reflectance in in the the red red band for band for the the central central pixel pixel of of the the study study area. area. The The shift shift and and scale scale fitting fitting parameters parameters are are indicated. indicated. 2.3.3. GPR LAI and Uncertainty Retrieval 2.3.3. GPR LAI and Uncertainty Retrieval GPR can learn the complex relationship between input and output variables based on a GPR can learn the complex relationship between input and output variables based on a probabilistic approach [49]. In addition to highly accurate retrievals, GPR provides pixel-wise probabilistic approach [49]. In addition to highly accurate retrievals, GPR provides pixel-wise confidence intervals of the estimations [7]. Details regarding GPR can be found in a series of studies confidence intervals of the estimations [7]. Details regarding GPR can be found in a series of studies from Verrelst et al. [7,46,54,55], and we provide only a brief description of GPR adopting their notation. from Verrelst et al. [7,46,54,55], and we provide only a brief description of GPR adopting their Under the GPR paradigm, the estimated LAI from GPR obey Gaussian distribution with mean notation. and variance given by the following: Under the GPR paradigm, the estimated LAI from GPR obey Gaussian distribution with mean and variance given by the following:   −1 µ∗ = K∗T K + σTn2 I y, (2) μ* = K* ( K + σ n2 I )-1 y , (2) σ 2 = K - K T ( K−+1σ n2 I ) -1 K* , (3) σ∗2 = K∗∗*− K∗T **K + *σn2 I K∗ , (3) where K* is a vector used for characterizing the similarity between the test and training samples, and where K K*prior ** is the is a covariance. vector usedAfor characterizing positive the similarity term is subtracted from the between the test and prior covariance to training estimate samples, the final and K variance.** is the prior covariance. A positive term is subtracted from Therefore, the posterior variance is always smaller than the prior variance, since the the prior covariance to estimate the final variance. observations Therefore, contribute some the additional posterior variance is always information [49].smaller than the prior Theoretically, variance, since the GPR-estimated the observations contribute some additional information [49]. uncertainty depends on the training reflectance. If the test reflectance is similar to the training Theoretically, the GPR-estimated uncertainty the reflectance, depends on thefor uncertainty training the testreflectance. set is relatively If thesmall test reflectance because weisare similar dealing to with the training similar reflectance, input featuresthe[46]. uncertainty for thestandard In this study, test set is deviation relatively (SD), small namely, because we the are squaredealing rootwith similar of the input predictive features [46]. In this study, standard deviation variance, was used as a proxy of the retrieval uncertainty. (SD), namely, the square root of the predictive variance, was used as a proxy In Equations (2)ofand the (3), retrieval K is auncertainty. kernel function that evaluates the similarity between the test In Equations reflectance x and all(2)training and (3), K is a kernel reflectance values. function In thisthat study,evaluates the similarity the automatic relevance between the test determination reflectance kernel was used: x and all training reflectance values. In this study, the automatic relevance determination kernel was used: B ( xi − x j ) b b 2 K ( xi , x j ) = v exp( −( b xi b−=1 x ) b 2 ) +  ij n2 , (4) j 2 b 2 B K ( xi , x j ) = v exp(−∑b=1 2 2 ) + δij σn , (4) where v is a scaling factor, B is the dimensions of the input2σ(in b this study, B = 6), σb controls the spread of the relations for each dimension of the input, and high values of σb mean lower informative content of the corresponding dimension for the regression and vice versa, σn is the noise standard deviation, and δij is Kronecker’s symbol. Model hyper-parameters involved in GPR were automatically optimized by maximizing the marginal likelihood in the training set [49]. In this study, GPR was implemented using the simpleR Remote Sens. 2019, 11, 244 7 of 18 where v is a scaling factor, B is the dimensions of the input (in this study, B = 6), σb controls the spread of the relations for each dimension of the input, and high values of σb mean lower informative content of the corresponding dimension for the regression and vice versa, σn is the noise standard deviation, and δij is Kronecker’s symbol. Model hyper-parameters involved in GPR were automatically optimized by maximizing the marginal likelihood in the training set [49]. In this study, GPR was implemented using the simpleR MATLAB toolbox, freely available at https://github.com/IPL-UV/simpleR (accessed on 18 January 2019). The training dataset is the prerequisite to train a GPR suitable for LAI and uncertainty retrieval. The 67 LAINet measurements and the associated reconstructed Landsat-like reflectances were therefore collected and added to the training dataset. Our study area, beyond the cropping period, is characterized by seasonal inundation and bare soil after harvest. Two dry-soil and two flooded-soil spectral reflectances were therefore added to the training database (with LAI values equal to zero) to improve its representativeness. This treatment was recommended to improve the accuracy of GPR, especially for the input reflectances without associated in situ measurements [34,46]. This resulted in the final training dataset made of 71 LAI-reflectance pairs. The trained GPR was used to generate the concomitant LAI and uncertainty maps with 30-m and 8-day resolutions. 2.4. Validation Approach The proposed framework approach requires accurate reconstruction of reflectance data at high spatial resolution. The accuracy of the reconstructed reflectances was first assessed (Section 3.1). The original reflectances from the six available Landsat 8 OLI images (Figure 2) were used as the reference data. We conducted a leave-one-out cross-validation to evaluate the performance of GPR in LAI retrieval (Section 3.2). The n = 67 number of pairs of LAINet-reflectances were split in n-1 for the training and 1 for the validation. Then, the calibrated GPR was used to predict LAI for the observation left out. This process was repeated for all the samples. Finally, we evaluated the performance comparing the predicted LAI with the ground observations for all the samples. For comparison purposes, we also assessed the performance of the exponential regression (ER) approach proposed in our previous work [8]. Both approaches share the same three-step flowchart described in Figure 3, but with some major updates: (1) the exponential regression (y = 0.21exp(3.44x)) proposed in [8] has now been replaced with GPR, to generate LAI and uncertainty maps from ground measurements and high spatiotemporal reflectance data; (2) spectral reflectance was directly ingested into GPR, instead of Normalized Difference Vegetation Index (NDVI) ingested into ER in the previous work; (3) SG filtering was used to generate the phenology model in CACAO, rather than the Double Logistic function adopted in the previous work. Finally, the temporal evolution of the predicted LAI and associated uncertainties was evaluated (Section 3.3). 3. Results 3.1. Evaluation of the Reconstructed High Spatiotemporal Satellite Data For brevity, Figure 5 only shows the histograms of residuals between the CACAO reconstructed reflectances and OLI observed reflectances on DOY 187 and 235, which are temporally synchronous with the LAINet measurements. The comparison for other dates is shown in Appendix A. In general, the distributions of residuals have peak shapes centered at 0 without significant bias, confirming a good agreement between the reconstructed and the observed reflectances. Remote Sens. 2019, 11, x FOR PEER REVIEW 8 of 18 Remote Sens. 2019, 11, 244 8 of 18 Remote Sens. 2019, 11, x FOR PEER REVIEW 8 of 18 Figure 5. Histograms of residuals between the CACAO reconstructed reflectances and the Landsat 8 Figure5.5.Histograms Figure Histograms of residuals between between the theCACAO CACAOreconstructed reconstructedreflectances reflectancesandandthe the Landsat Landsat 8 8 OLI observed reflectances on DOY 187 (solid lines) and 235 (dashed lines) for the different spectral OLI OLI observed reflectancesonon observed reflectances DOYDOY 187187 (solid (solid lines)lines) and(dashed and 235 235 (dashed lines) forlines) for the spectral the different different bands. The statistics refer to the relative differences (RD). spectral bands. bands. The The statistics statistics refer refer to the to thedifferences relative relative differences (RD). (RD). 3.2. 3.2. Cross-validation 3.2.Cross-Validation Cross-validation The The results of the cross-validation cross-validation forER ER andGPRGPR approachesare are shownininFigure Figure 6. Compared Theresults results of the cross-validationfor for ERand and GPRapproaches approaches areshown shown in Figure6.6. Compared Compared to to ourour previous ER approach [8] 2 (R 2 = 0.70, RMSE = 0.52, relative RMSE (R-RMSE) = 18.4%), GPR to our previous ER approach [8] (R = 0.70, RMSE = 0.52, relative RMSE (R-RMSE) = 18.4%),showed previous ER approach [8] (R = 0.70, 2 RMSE = 0.52, relative RMSE (R-RMSE) = 18.4%), GPR GPR showed showedimproved improved improved performances performances performances (R (R2 = 0.88, 2 = 0.88, RMSE = 0.36, relative RMSE = 12.7%). In addition, ER (R2RMSE = 0.88,=RMSE 0.36, relative RMSE RMSE = 0.36, relative = 12.7%). In addition, = 12.7%). ER clearly In addition, ER clearly clearlyunderestimated underestimated LAI LAIobservations LAI observations underestimated for LAIfor observations LAI values for values higher higherhigher LAI values than than 4.0 than 4.0because because 4.0 because ofthe thesaturation saturation of theofsaturation ofof of NDVI. NDVI. This NDVI. This Thisunderestimation underestimation phenomenon phenomenon underestimation for the for phenomenon high for the the high LAI LAI values values high LAI values was wasclearly clearly was clearly mitigatedmitigated ininGPR. in GPR. mitigated GPR. Figure Figure 6.6.Cross-validation Figure Cross-validation of(a) of (a)exponential 6. Cross-validation exponential of regression regression (a) exponential (ER) (ER) predicted regressionpredicted LAI, LAI,and (ER) predicted and(b) (b) LAI, Gaussian Gaussian and process process (b) Gaussian regression regression process regression (GPR)with (GPR) withfield (GPR) fieldLAI with field LAIobservations. LAI observations.The observations. The dashed dashed line The dashed line represents line represents the the 1:1represents 1:1line. line. the 1:1 line. 3.3. Temporal 3.3.Temporal Evolution TemporalEvolution Evolutionof of LAI ofLAI and LAIand Uncertainty and Uncertainty Retrievals Uncertainty Retrievals Retrievals 3.3. The LAI TheLAI maps LAImaps (Figure maps(Figure (Figure7)7) provide detailed 7)provide provide detailed informationofofthe the expectedspatiotemporal spatiotemporal dynamics The detailed information information of theexpected expected spatiotemporaldynamics dynamics ofofvegetation in the study area. Before green-up (DOY 161), consistent low LAIvalues values close to zero of vegetation in the study area. Before green-up (DOY 161), consistent low LAI values closetotozero vegetation in the study area. Before green-up (DOY 161), consistent low LAI close zero are observed. areobserved. After observed.After green-up, Aftergreen-up, the green-up, the vegetation the vegetation grew vegetation grew rapidly, grew rapidly, and reached the growing peak on DOY 217. are rapidly, and and reached reachedthe thegrowing growingpeak peakononDOY DOY During this period 217. During (from(from this period DOYDOY201 to 201217), areas to 217), around areas the the around water (in (in water thethe northwestern northwestern part ofof part the 217. During this period (from DOY 201 to 217), areas around the water (in the northwestern part of study area)area) the study werewerecharacterized by low characterized LAI LAI by low values (<2.0) values duedue (<2.0) to later sowing. to later After sowing. the the After growing peak growing the study area) were characterized by low LAI values (<2.0) due to later sowing. After the growing on DOY peak on 217, DOYthe 217,vegetation beganbegan the vegetation to wither, and and to wither, the LAI decreased the LAI decreasedgradually to zero. gradually During to zero. Duringthis peak on DOY 217, the vegetation began to wither, and the LAI decreased gradually to zero. During Remote Sens. Remote 2019,11, 2019, Sens. 11,244 x FOR PEER REVIEW 9 of1818 9 of this withering period, areas around the water were also characterized by rather low LAI values, withering period, because of areas seasonal around the water were also characterized by rather low LAI values, because of inundation. seasonal inundation. Figure7.7.Maps Figure Mapsofofthe theestimated estimatedLAI LAIatathigh highspatiotemporal spatiotemporal resolution resolution (30-m (30-mand and8-day). 8-day).The Thewhite white areasrepresent areas representnon-vegetated non-vegetatedcover covertypes. types. The Theretrieved retrieveduncertainty uncertaintymaps maps (Figure (Figure8) 8) provide insight provide into insight thethe into confidence confidence of the estimated of the LAI. estimated LAI. Before green-up, the middle area covering the LAINet nodes was characterized Before green-up, the middle area covering the LAINet nodes was characterized by high uncertainty by high uncertainty (SD (SDmaps > 0.8). LAI >0.8).near LAI and maps near and during the during period the period of implementation of implementation of the(from of the LAINet LAINet (from DOY 185 toDOY 233) 185 weretomore 233) were more credible credible (SD < 0.7). (SD <0.7). During During the the dormancy dormancy period period (e.g., DOY (e.g., DOY 289-321), 289-321), LAI LAI retrievals inretrievals in the region the inundated inundated wereregion were(SD uncertain uncertain > 0.9).(SD >0.9). Remote Sens. 2019, 11, 244 10 of 18 Remote Sens. 2019, 11, x FOR PEER REVIEW 10 of 18 Figure Figure 8. Maps 8. Maps ofstandard of the the standard deviation deviation (SD) associated (SD) associated with with the theretrievals LAI LAI retrievals (Figure (Figure 7). The7). The white white areas areas represent represent non-vegetated non-vegetated cover cover types. types. TheThe temporal temporal evolution evolution of GPR of GPR andand ERER retrievals retrievals was was further further evaluated evaluated over over thethe nodes nodes 3, 3, 5, 5, 8 8 and and1212ofof LAINet, LAINet, which whichprovided provided full temporal full temporal coverage coverage forfor field measurements field measurements (Figure (Figure 9).9). In In general, both the general, bothGPRthe and GPRER andestimates had similar ER estimates temporal had similar patterns temporal and followed patterns and followedthe temporal the temporaldynamics of the observed dynamics LAI. However, of the observed three main LAI. However, three differences can be main differences clearly can observed: be clearly (1) the observed: (1) local temporal the local temporal variations variations in groundin ground LAI measurements LAI measurements were were well captured well captured by GPR, by GPR, whereas whereas the ER theover-smoothed ER over- smoothed the LAIthe LAIseries; time time series; (2) the(2) the magnitude magnitude of LAIof values LAI values beforebefore green-up green-up (DOY (DOY <161) < 161) was was closeto zero close to zero for GPR for GPR as expected, as expected, whereas whereas ER ER slightly slightly over-estimated over-estimated thethe LAI; LAI; (3)(3) during during thethe witheringperiod withering period between between DOYDOY 257 and 257 and 297, 297, the GPR-estimated the GPR-estimated LAIsLAIs werewere largerlarger thanthan the ER-estimated the ER-estimated ones.ones. Note that Note that nearly nearly all fieldallobservations field observations (except(except the one the forone for8plot plot 8 on DOY on DOY 185) 185) lie lie within within the confidence the confidence interval of interval GPR of GPR predictions. predictions. Remote Sens. 2019, 11, 244 11 of 18 Remote Sens. 2019, 11, x FOR PEER REVIEW 11 of 18 Figure Figure 9. The temporal 9. The temporal evolution evolution of of the the LAI derived from Gaussian process regression (GPR), exponential (ER), and exponential regression (ER), andLAINet LAINetfield fieldobservations observationsforfor (a)(a) node node 3, (b) 3, (b) node node 5, node 5, (c) (c) node 8, 8, and and (d) node (d) node 12 (see 12 (see Figure Figure 1 for1node for node location). location). The The shaded shaded areas areas represent represent the GPR the GPR uncertainty. uncertainty. 4. Discussion 4. Discussion 4.1. Improvements Over Previous Approach 4.1. Improvements Over Previous Approach This study is based on our previous approach [8] with three major updates: (1) the exponential This study is based on our previous approach [8] with three major updates: (1) the exponential regression was replaced with GPR, to generate LAI and uncertainty maps from ground measurements regression was replaced with GPR, to generate LAI and uncertainty maps from ground and high spatiotemporal reflectance data; (2) spectral reflectance was directly ingested into the transfer measurements and high spatiotemporal reflectance data; (2) spectral reflectance was directly ingested function, instead of the vegetation index in the previous work; (3) SG filtering was used to generate into the transfer function, instead of the vegetation index in the previous work; (3) SG filtering was the phenology model in CACAO, rather than the Double Logistic function adopted previously. used to generate the phenology model in CACAO, rather than the Double Logistic function adopted These updates bring several improvements. First, the accuracy of the retrievals has improved. previously. These updates bring several improvements. First, the accuracy of the retrievals has The R2 increased from 0.70 to 0.88, and the RMSE decreased from 0.52 LAI (18.4%) to 0.36 LAI improved. The R2 increased from 0.70 to 0.88, and the RMSE decreased from 0.52 LAI (18.4%) to 0.36 (12.7%), meeting the accuracy requirements proposed by the Global Climate Observing System [44]. LAI (12.7%), meeting the accuracy requirements proposed by the Global Climate Observing System Second, the temporal evolution of LAI is better captured by GPR, whereas ER over-smooths the [44]. Second, the temporal evolution of LAI is better captured by GPR, whereas ER over-smooths the time series. Further, GPR shows a better temporal extrapolation capacity than ER, particularly for the time series. Further, GPR shows a better temporal extrapolation capacity than ER, particularly for the dormancy period with LAI retrievals close to the expected zero LAI value. This allows the generation of dormancy period with LAI retrievals close to the expected zero LAI value. This allows the generation temporally continuous LAI maps over the entire growing season. Finally, GPR provides the uncertainty of temporally continuous LAI maps over the entire growing season. Finally, GPR provides the associated with the LAI retrievals [46,54]. uncertainty associated with the LAI retrievals [46,54]. Several reasons can explain this improvement. (1) GPR is an advanced regression method and Several reasons can explain this improvement. (1) GPR is an advanced regression method and can learn the complex nonlinear relationship between remote sensing observations and biophysical can learn the complex nonlinear relationship between remote sensing observations and biophysical variables [48,49]. Many studies have shown that GPR outperforms parametric regression methods variables [48,49]. Many studies have shown that GPR outperforms parametric regression methods in in most cases [5,7,34]. Our study further confirms this finding. (2) GPR has the capability to directly most cases [5,7,34]. Our study further confirms this finding. (2) GPR has the capability to directly exploit the full spectrum [55] and, for this reason, individual bands rather than the derived vegetation exploit the full spectrum [55] and, for this reason, individual bands rather than the derived vegetation index were used in our GPR. The procedure for calculating the vegetation index (i.e., the selection index were used in our GPR. The procedure for calculating the vegetation index (i.e., the selection of of bands and the formulation of the selected bands) would cause loss of information useful for LAI bands and the formulation of the selected bands) would cause loss of information useful for LAI retrieval [56]. (3) SG filtering can reserve the relatively high frequency variation in the original retrieval [56]. (3) SG filtering can reserve the relatively high frequency variation in the original reflectance time series [53], which is therefore suitable for capturing the local scale changes in the time reflectance time series [53], which is therefore suitable for capturing the local scale changes in the series of LAI and preventing the over-smoothing observed in Double Logistic. (4) Two dry-soil and time series of LAI and preventing the over-smoothing observed in Double Logistic. (4) Two dry-soil two flooded-soil spectral reflectances (with LAI values equal to zero) were introduced to establish an and two flooded-soil spectral reflectances (with LAI values equal to zero) were introduced to extended training dataset. The original training dataset consisting of only the 67 LAINet measurements establish an extended training dataset. The original training dataset consisting of only the 67 LAINet measurements under-sampled the low LAI values. The introduction of dry- and flooded-soil spectral Remote Sens. 2019, 11, 244 12 of 18 under-sampled the low LAI values. The introduction of dry- and flooded-soil spectral reflectances improved the representativeness of the trained GPR for low LAI values, and therefore improved its temporal extrapolation capacity especially for dormancy periods. The GPR-estimated LAI values are higher than the ER retrievals during the withering period (see Figure 9). One possible explanation is that the NDVI represents the greenness of the vegetation [57], and the gradual decrease of leaf chlorophyll content during the withering period [58] can therefore decrease the value of NDVI and further cause an underestimation of the LAI value. The temporal variation of leaf chlorophyll content (which further causes the variation of leaf spectra) underlines the importance of the direct use of multi- or even hyper-spectral reflectances, especially those containing information on leaf spectra. 4.2. Potential Applications The LAI time series generated from the proposed framework are characterized by high spatiotemporal resolution and are accompanied with concomitant uncertainty maps. These datasets may be useful in many applications. The generated high spatiotemporal resolution LAI maps can support the validation activities for the coarse resolution LAI products. Our framework complies with the hierarchical four-stage validation approach proposed by the Land Product Validation (LPV) subgroup of the Committee Earth Observing Satellites’ Working Group on Calibration and Validation (CEOS WGCV) [59]. Therefore, the generated LAI maps can directly serve as reference maps in validation activities. The obvious advantage of our framework over existing studies is that it can derive temporally continuous LAI reference maps. Therefore, it can be used to assess the temporal performance of the coarse resolution LAI products, which is the prerequisite for long-term global change research [60,61]. The derived high spatiotemporal resolution LAI maps can also be used to quantify the spatial and temporal variability of crop distribution and status at field scale. This kind of information is of paramount importance for precision agriculture [62]. The high spatiotemporal resolution LAI maps can provide guidance on agricultural activities and aid farmers to decide, for example, which is the appropriate time to fertilize or irrigate their fields. The retrieved uncertainty information would allow a proper use of the LAI retrieval. For example, the high reliable retrieval with low SD can be easily extracted for subsequent applications. The inconsistencies among existing coarse resolution LAI products call for quality-assured long-term LAI products [60,63]. Our framework may provide an effective method to fulfill this target. In addition, the concomitant uncertainty maps facilitate the integration of remote sensing and process models at field scale [64]. For example, the LAI retrievals can be easily assimilated by crop models, because their contribution in the assimilation process can be automatically weighted according to the concomitant uncertainty maps. 4.3. Future Prospects of Research This study proposed a framework to generate concomitant LAI and uncertainty maps with high spatiotemporal resolution. Several issues are worth being noted to further improve the framework. High spatial resolution remote sensing observations that are temporally synchronous with field measurements are needed to reduce the uncertainty during LAI estimation caused by the dynamic change of vegetation [65]. The tradeoff between temporal and spatial resolutions of a remote sensing sensor makes it difficult to acquire remote sensing observation with both high spatial resolution and frequent coverage. Data blending techniques were used here to obtain high spatial resolution remote sensing observations that were temporally synchronous with LAINet observations. In our study area, one year of MODIS data was enough to provide temporally continuous coverage although multi-year data may be required to build the phenology model in regions with high cloud coverage. Note that a spectral normalization between MCD43A4 and Landsat 8 OLI data was not applied. We assumed that the difference in their spectral response functions could be implicitly accounted for during the Remote Sens. 2019, 11, 244 13 of 18 shifting and scaling procedure, considering the significant linear correlation between MODIS and OLI reflectance products [51]. An alternative to generate high spatiotemporal resolution remote sensing data is the combined use of multi-sensor data. For example, Zhao et al. [16] used HJ1/CCD and Landsat 8/OLI simultaneously to fill temporal gaps which existed in any one single data source. In addition, some newly launched satellites (e.g., GF-6 or Venµs) can provide remote sensing observations with both high spatial and temporal resolutions. The performance of our proposed method for multi-sensor or high spatiotemporal data will be evaluated in future studies. The analysis of the uncertainty maps suggested that under-sampled observations (e.g., middle areas before green-up and inundated areas during the dormancy period, see Figure 9) are often associated with high uncertainty. Therefore, a longer temporal span and more below nodes should be added to the LAINet system in our study area to enhance the representativeness of the training dataset [50,66]. To partly cope with the training dataset under-sample for low LAI, we deliberately extended the original training dataset with four new spectra that covered two main kinds of non-vegetated surfaces (i.e., dry and flooded soils) with an LAI value of zero. This extended training dataset improved the extrapolation capacity of GPR, especially for low values which were not sufficiently sampled by LAINet. Surprisingly, the estimates with low LAI also represent high uncertainty after introducing the non-vegetated spectra. This may be because the non-vegetated and low-LAI surfaces have significant heterogeneity [67] and four spectra cannot represent them very well. However, we did not try to add more non-vegetated spectra, because the number of LAINet measurements was limited (67), and more non-vegetated spectra would have diluted the information involved in LAINet measurements, and further reduced the retrieval accuracy. Another solution to improve the representativeness of the training dataset is to add reflectance simulated by RT models, which can cover a broader range of scenarios. As proposed by [12], the simulated reflectance can be imported through the joint Gaussian process to extrapolate the observed reflectance. The predictive variance from GPR is an indicator of the real uncertainty; however, the physical meaning of this predictive variance is still not clear [46]. Intuitively, it is an indicator of the similarity between the test and the training samples in the LAI space transferred from the reflectance space through a kernel function (see Equation (4)). The kernel used assumed that the measurement noise was independent from the signal (see Equation (4)), which does not hold in most of the problems. The heteroscedastic Gaussian processes [68] can be employed to further improve the estimation results in future studies. In addition, several other methods exist in the literature to quantify the uncertainty associated with retrievals including partial least squares regression [69], Markov chain Monte Carlo [70], random forest [71], and Kriging [72]. These methods will be assessed in future studies. 5. Conclusions We proposed a framework to generate long-term time series of LAI and their associated uncertainty maps with high spatiotemporal resolution. It was based on the combination of LAINet, CACAO, and GPR. The performance of the proposed method was evaluated over the entire growing season in a crop site in northern China. The accuracy of the retrieved LAI maps (0.36 LAI (12.7%) in terms of RMSE with R2 = 0.88 correlation) met the up-to-date uncertainty threshold (15%) proposed by the Global Climate Observing System. To our knowledge, this is the first work combining wireless sensor network, data blending, and machine learning technologies for retrieving LAI and its uncertainty at high spatiotemporal resolution. Our framework will contribute to precision agriculture, as well as to the retrieval and validation of LAI products. Author Contributions: G.Y. conceived the experiments and wrote the manuscript; Y.Q. implemented the field experiments; A.V., W.Z., B.X., Y.Z., K.L., J.L. and Q.L. reviewed and edited the manuscript. Funding: This work was funded by the National Natural Science Foundation of China (41531174), the GF6 Project (30-Y20A03-9003-17/18), the National Natural Science Foundation of China (41601403 and 41631180), the Youth Innovation Promotion Association CAS (2016333), the China Postdoctoral Science Foundation (2018T110996), Remote Sens. 2019, 11, 244 14 of 18 the “Innovation Remote Sens. 2019, Ability 11, x FORPromotion” Program of the Sichuan Provincial Department of Finance (2016QNJJ-023), PEER REVIEW 14 of 18 and the EC Copernicus Global Land Service (CGLOPS-1, 199494-JRC). Acknowledgments: The authors would Acknowledgments: would like like to to acknowledge acknowledgethe theImage ImageProcessing ProcessingLaboratory, Laboratory,Universitat Universitatdede València, València,for forthe thesimpleR simpleR toolbox. toolbox. Conflicts of Interest: The authors declare no conflict of interest. Conflicts of Interest: The authors declare no conflict of interest. Appendix A Appendix A Figure A1. Histograms Figure A1. Histograms ofof residuals residuals between between the the CACAO CACAO reconstructed reconstructed reflectances reflectances and and the the Landsat Landsat 8 OLI observed reflectances on DOY 139 (red solid lines), 235 (blue dashed lines), 299 (orange dotted 8 OLI observed reflectances on DOY 139 (red solid lines), 235 (blue dashed lines), 299 (orange dotted lines), and 315 (green dash-dotted lines) for the different spectral bands. The statistics refer to the lines), and 315 (green dash-dotted lines) for the different spectral bands. The statistics refer to the relative differences (RD). (RD). relative differences References References 1. Wang, G.L.; Yu, M.; Pal, J.S.; Mei, R.; Bonan, G.B.; Levis, S.; Thornton, P.E. On the development of a coupled 1. Wang, G.L.; Yu, M.; Pal, J.S.; Mei, R.; Bonan, G.B.; Levis, S.; Thornton, P.E. On the development of a coupled regional climate-vegetation model RCM-CLM-CN-DV and its validation in Tropical Africa. Clim. Dyn. 2016, regional climate-vegetation model RCM-CLM-CN-DV and its validation in Tropical Africa. Clim. Dyn. 46, 515–539. 2016, [CrossRef] 46, 515–539, doi:10.1007/s00382-015-2596-z. 2. 2. Liu, Liu, J.; Chen, J.M.; J.; Chen, Cihlar, J.; J.M.; Cihlar, J.; Park, Park, W.M. W.M. A A process-based process-based boreal borealecosystem ecosystemproductivity productivitysimulator simulatorusing using remote sensing inputs. Remote Sens. Environ. 1997, 62, 158–175. [CrossRef] remote sensing inputs. Remote Sens. Environ. 1997, 62, 158–175, doi:10.1016/s0034-4257(97)00089-8. 3. 3. Serbin, Serbin, S.P.; S.P.; Ahl, Ahl, D.E.; D.E.; Gower, Gower, S.T.S.T. Spatial and temporal Spatial and temporalvalidation validationofofthetheMODIS MODISLAI LAIandandFPAR FPARproducts products across a boreal forest wildfire chronosequence. Remote Sens. Environ. across a boreal forest wildfire chronosequence. Remote Sens. Environ. 2013, 133, 71–84, 2013, 133, 71–84. [CrossRef] 4. Myneni, R.B.; Hoffman, S.; Knyazikhin, Y.; Privette, J.L.; Glassy, J.; Tian, Y.; Wang, Y.; Song, X.; Zhang, Y.; doi:10.1016/j.rse.2013.01.022. 4. Smith, Myneni, R.B.;etHoffman, G.R.; al. Global S.; products Knyazikhin, of vegetation Y.; Privette,leaf J.L.;area andJ.;fraction Glassy, Tian, Y.;absorbed Wang, Y.;PAR from Song, year one X.; Zhang, Y.;of MODIS data. Remote Sens. Environ. 2002, 83, 214–231. [CrossRef] Smith, G.R.; et al. Global products of vegetation leaf area and fraction absorbed PAR from year one of 5. García-Haro, MODIS data. Remote F.J.; Campos-Taberner, Sens. Environ. 2002,M.; 83, Muñoz-Marí, J.; Laparra, V.; Camacho, F.; Sánchez-Zapero, J.; 214–231, doi:10.1016/s0034-4257(02)00074-3. 5. Camps-Valls, G. Derivation of global vegetation García-Haro, F.J.; Campos-Taberner, M.; Muñoz-Marí, J.; Laparra, biophysical parameters from EUMETSAT V.; Camacho, Polar System. F.; Sánchez-Zapero, J.; ISPRS J. Photogramm. Remote Sens. 2018, 139, 57–74. [CrossRef] Camps-Valls, G. Derivation of global vegetation biophysical parameters from EUMETSAT Polar System. 6. Heiskanen, J.; Rautiainen, ISPRS J. Photogramm. M.; Stenberg, Remote P.; 139, Sens. 2018, Mottus, M.;doi:10.1016/j.isprsjprs.2018.03.005. 57–74, Vesanto, V.H. Sensitivity of narrowband vegetation 6. indices Heiskanen,to boreal forest LAI, J.; Rautiainen, M.;reflectance Stenberg, seasonality P.; Mottus, M.;andVesanto, species composition. V.H. Sensitivity ISPRS J. Photogramm. of narrowband 2013, 78, vegetation 1–14. [CrossRef] indices to boreal forest LAI, reflectance seasonality and species composition. ISPRS J. Photogramm. 2013, 7. Verrelst, 78, 1–14, J.; Munoz, J.; Alonso, L.; Delegido, J.; Rivera, J.P.; Camps-Valls, G.; Moreno, J. Machine learning doi:10.1016/j.isprsjprs.2013.01.001. 7. regression algorithms Verrelst, J.; Munoz, for biophysical J.; Alonso, L.; Delegido,parameter retrieval: J.; Rivera, J.P.; Opportunities Camps-Valls, G.; Moreno, forJ. Sentinel-2 and -3. Machine learning Remote regressionSens.algorithms Environ. 2012, 118, 127–139. for biophysical [CrossRef] parameter retrieval: Opportunities for Sentinel-2 and -3. Remote Sens. Environ. 2012, 118, 127–139, doi:10.1016/j.rse.2011.11.002. 8. Yin, G.F.; Li, A.N.; Jin, H.A.; Zhao, W.; Bian, J.H.; Qu, Y.H.; Zeng, Y.L.; Xu, B.D. Derivation of temporally continuous LAI reference maps through combining the LAINet observation system with CACAO. Agric. For. Meteorol. 2017, 233, 209–221, doi:10.1016/j.agrformet.2016.11.267. Remote Sens. 2019, 11, 244 15 of 18 8. Yin, G.F.; Li, A.N.; Jin, H.A.; Zhao, W.; Bian, J.H.; Qu, Y.H.; Zeng, Y.L.; Xu, B.D. Derivation of temporally continuous LAI reference maps through combining the LAINet observation system with CACAO. Agric. For. Meteorol. 2017, 233, 209–221. [CrossRef] 9. Houborg, R.; McCabe, M.; Cescatti, A.; Gao, F.; Schull, M.; Gitelson, A. Joint leaf chlorophyll content and leaf area index retrieval from Landsat data using a regularized model inversion system (REGFLEC). Remote Sens. Environ. 2015, 159, 203–221. [CrossRef] 10. Laurent, V.C.E.; Schaepman, M.E.; Verhoef, W.; Weyermann, J.; Chavez, R.O. Bayesian object-based estimation of LAI and chlorophyll from a simulated Sentinel-2 top-of-atmosphere radiance image. Remote Sens. Environ. 2014, 140, 318–329. [CrossRef] 11. Zeng, Y.; Li, J.; Liu, Q.; Huete, A.R.; Yin, G.; Xu, B.; Fan, W.; Zhao, J.; Yan, K.; Mu, X. A Radiative Transfer Model for Heterogeneous Agro-Forestry Scenarios. IEEE Trans. Geosci. Remote Sens. 2016, 54, 4613–4628. [CrossRef] 12. Svendsen, D.H.; Martino, L.; Campos-Taberner, M.; Garcia-Haro, F.J.; Camps-Valls, G. Joint Gaussian Processes for Biophysical Parameter Retrieval. IEEE Trans. Geosci. Remote Sens. 2017, 1–10. [CrossRef] 13. Teubner, I.E.; Forkel, M.; Jung, M.; Liu, Y.Y.; Miralles, D.G.; Parinussa, R.; van der Schalie, R.; Vreugdenhil, M.; Schwalm, C.R.; Tramontana, G.; et al. Assessing the relationship between microwave vegetation optical depth and gross primary production. Int. J. Appl. Earth Obs. Geoinf. 2018, 65, 79–91. [CrossRef] 14. Verrelst, J.; Camps-Valls, G.; Munoz-Mari, J.; Rivera, J.P.; Veroustraete, F.; Clevers, J.; Moreno, J. Optical remote sensing and the retrieval of terrestrial vegetation bio-geophysical properties—A review. ISPRS J. Photogramm. 2015, 108, 273–290. [CrossRef] 15. Zeng, Y.; Xu, B.; Yin, G.; Wu, S.; Hu, G.; Yan, K.; Yang, B.; Song, W.; Li, J. Spectral Invariant Provides a Practical Modeling Approach for Future Biophysical Variable Estimations. Remote Sens. 2018, 10, 1508. [CrossRef] 16. Zhao, J.; Li, J.; Liu, Q.; Fan, W.; Zhong, B.; Wu, S.; Yang, L.; Zeng, Y.; Xu, B.; Yin, G. Leaf Area Index Retrieval Combining HJ1/CCD and Landsat8/OLI Data in the Heihe River Basin, China. Remote Sens. 2015, 7, 6862–6885. [CrossRef] 17. Yin, G.F.; Li, J.; Liu, Q.H.; Fan, W.L.; Xu, B.D.; Zeng, Y.L.; Zhao, J. Regional leaf area index retrieval based on remote sensing: The role of radiative transfer model selection. Remote Sens. 2015, 7, 4604–4625. [CrossRef] 18. Quan, X.W.; He, B.B.; Li, X. A Bayesian network-based method to alleviate the ill-posed inverse problem: A case study on leaf area index and canopy water content retrieval. IEEE Trans. Geosci. Remote Sens. 2015, 53, 6507–6517. [CrossRef] 19. Zeng, Y.; Li, J.; Liu, Q.; Huete, A.R.; Xu, B.; Yin, G.; Zhao, J.; Yang, L.; Fan, W.; Wu, S.; et al. An Iterative BRDF/NDVI Inversion Algorithm Based on A Posteriori Variance Estimation of Observation Errors. IEEE Trans. Geosci. Remote Sens. 2016, 54, 6481–6496. [CrossRef] 20. Breda, N.J.J. Ground-based measurements of leaf area index: A review of methods, instruments and current controversies. J. Exp. Bot. 2003, 54, 2403–2417. [CrossRef] 21. Jonckheere, I.; Fleck, S.; Nackaerts, K.; Muys, B.; Coppin, P.; Weiss, M.; Baret, F. Review of methods for in situ leaf area index determination—Part I. Theories, sensors and hemispherical photography. Agric. For. Meteorol. 2004, 121, 19–35. [CrossRef] 22. Weiss, M.; Baret, F.; Smith, G.J.; 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. [CrossRef] 23. Zeng, Y.; Li, J.; Liu, Q.; Hu, R.; Mu, X.; Fan, W.; Xu, B.; Yin, G.; Wu, S. Extracting Leaf Area Index by Sunlit Foliage Component from Downward-Looking Digital Photography under Clear-Sky Conditions. Remote Sens. 2015, 7, 13410–13435. [CrossRef] 24. Demarez, V.; Duthoit, S.; Baret, F.; Weiss, M.; Dedieu, G. Estimation of leaf area and clumping indexes of crops with hemispherical photographs. Agric. For. Meteorol. 2008, 148, 644–655. [CrossRef] 25. Leblanc, S.G. Correction to the plant canopy gap-size analysis theory used by the Tracing Radiation and Architecture of Canopies instrument. Appl. Opt. 2002, 41, 7667–7670. [CrossRef] [PubMed] 26. LI-COR. LAI-2000 Plant Canopy Analyzer Operating Manual. 1991. Available online: http://www.ecotek. com.cn/download/Manual-LAI-2000-EN.pdf (accessed on 23 January 2019). Remote Sens. 2019, 11, 244 16 of 18 27. Ryu, Y.; Lee, G.; Jeon, S.; Song, Y.; Kimm, H. Monitoring multi-layer canopy spring phenology of temperate deciduous and evergreen forests using low-cost spectral sensors. Remote Sens. Environ. 2014, 149, 227–238. [CrossRef] 28. Campos-Taberner, M.; Garcia-Haro, F.J.; Confalonieri, R.; Martinez, B.; Moreno, A.; Sanchez-Ruiz, S.; Gilabert, M.A.; Camacho, F.; Boschetti, M.; Busetto, L. Multitemporal monitoring of plant area index in the Valencia rice district with PocketLAI. Remote Sens. 2016, 8, 202. [CrossRef] 29. Qu, Y.H.; Wang, J.; Song, J.L.; Wang, J.D. Potential and limits of retrieving conifer leaf area index using smartphone-based method. Forests 2017, 8, 217. [CrossRef] 30. Yan, K.; Park, T.; Yan, G.J.; Chen, C.; Yang, B.; Liu, Z.; Nemani, R.R.; Knyazikhin, Y.; Myneni, R.B. Evaluation of MODIS LAI/FPAR product collection 6. Part 1: Consistency and improvements. Remote Sens. 2016, 8, 359. [CrossRef] 31. Liang, S.L.; Zhao, X.; Liu, S.H.; Yuan, W.P.; Cheng, X.; Xiao, Z.Q.; Zhang, X.T.; Liu, Q.; Cheng, J.; Tang, H.R.; et al. A long-term Global LAnd Surface Satellite (GLASS) data-set for environmental studies. Int. J. Digit. Earth 2013, 6, 5–33. [CrossRef] 32. Baret, F.; Weiss, M.; Lacaze, R.; Camacho, F.; Makhmara, H.; Pacholcyzk, P.; Smets, B. GEOV1: LAI and FAPAR essential climate variables and FCOVER global time series capitalizing over existing products. Part1: Principles of development and production. Remote Sens. Environ. 2013, 137, 299–309. [CrossRef] 33. Verger, A.; Baret, F.; Weiss, M. Near real-time vegetation monitoring at global scale. IEEE J. Sel. Top. Appl. Earth Obs. Remote. Sens. 2014, 7, 3473–3481. [CrossRef] 34. Campos-Taberner, M.; García-Haro, F.J.; Camps-Valls, G.; Grau-Muedra, G.; Nutini, F.; Crema, A.; Boschetti, M. Multitemporal and multiresolution leaf area index retrieval for operational local rice crop monitoring. Remote Sens. Environ. 2016, 187, 102–118. [CrossRef] 35. Ganguly, S.; Nemani, R.R.; Zhang, G.; Hashimoto, H.; Milesi, C.; Michaelis, A.; Wang, W.; Votava, P.; Samanta, A.; Melton, F. Generating global leaf area index from Landsat: Algorithm formulation and demonstration. Remote Sens. Environ. 2012, 122, 185–202. [CrossRef] 36. Houborg, R.; McCabe, M.F. A hybrid training approach for leaf area index estimation via Cubist and random forests machine-learning. ISPRS J. Photogramm. Remote Sens. 2018, 135, 173–188. [CrossRef] 37. Zhu, X.; Helmer, E.H.; Gao, F.; Liu, D.; Chen, J.; Lefsky, M.A. A flexible spatiotemporal method for fusing satellite images with different resolutions. Remote Sens. Environ. 2016, 172, 165–177. [CrossRef] 38. Xu, B.; Li, J.; Liu, Q.; Huete, A.R.; Yu, Q.; Zeng, Y.; Yin, G.; Zhao, J.; Yang, L. Evaluating Spatial Representativeness of Station Observations for Remotely Sensed Leaf Area Index Products. IEEE J. Sel. Top. Appl. Earth Obs. Remote. Sens. 2016, 9, 3267–3282. [CrossRef] 39. Yin, G.F.; Li, J.; Liu, Q.H.; Li, L.H.; Zeng, Y.L.; Xu, B.D.; Yang, L.; Zhao, J. Improving leaf area index retrieval over heterogeneous surface by integrating textural and contextual information: A case study in the Heihe River Basin. IEEE Geosci. Remote Sens. Lett. 2015, 12, 359–363. [CrossRef] 40. Zeng, Y.L.; Li, J.; Liu, Q.H.; Li, L.H.; Xu, B.D.; Yin, G.F.; Peng, J.J. A Sampling Strategy for Remotely Sensed LAI Product Validation Over Heterogeneous Land Surfaces. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2014, 7, 3128–3142. [CrossRef] 41. Dou, B.; Wen, J.; Li, X.; Liu, Q.; Peng, J.; Xiao, Q.; Zhang, Z.; Tang, Y.; Wu, X.; Lin, X.; et al. Wireless Sensor Network of Typical Land Surface Parameters and Its Preliminary Applications for Coarse-Resolution Remote Sensing Pixel. Int. J. Distrib. Sens. Netw. 2016. [CrossRef] 42. Qu, Y.H.; Han, W.C.; Fu, L.Z.; Li, C.R.; Song, J.L.; Zhou, H.M.; Bo, Y.C.; Wang, J.D. LAINet—A wireless sensor network for coniferous forest leaf area index measurement: Design, algorithm and validation. Comput. Electron. Agric. 2014, 108, 200–208. [CrossRef] 43. Verger, A.; Baret, F.; Weiss, M.; Kandasamy, S.; Vermote, E. The CACAO method for smoothing, gap filling, and characterizing seasonal anomalies in satellite time series. IEEE Trans. Geosci. Remote Sens. 2013, 51, 1963–1972. [CrossRef] 44. GCOS. The Global Observing System for Climate: Implementation Needs. 2016. Available online: https://unfccc. int/sites/default/files/gcos_ip_10oct2016.pdf (accessed on 23 January 2019). 45. Fang, H.; Wei, S.; Jiang, C.; Scipal, K. Theoretical uncertainty analysis of global MODIS, CYCLOPES, and GLOBCARBON LAI products using a triple collocation method. Remote Sens. Environ. 2012, 124, 610–621. [CrossRef] Remote Sens. 2019, 11, 244 17 of 18 46. Verrelst, J.; Rivera, J.P.; 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. [CrossRef] 47. Han, W.C.; Qu, Y.H. Data uncertainty in an improved Bayesian network and evaluations of the credibility of the retrieved multitemporal high-spatial-resolution leaf area index. IEEE J. Sel. Top. Appl. Earth Obs. Remote. Sens. 2016, 9, 3553–3563. [CrossRef] 48. Camps-Valls, G.; Gómez-Chova, L.; Muñoz-Marí, J.; Vila-Francés, J.; Amorós, J.; Valle-Tascon, S.D.; Calpe-Maravilla, J. Biophysical parameter estimation with adaptive Gaussian processes. In Proceedings of the 2009 IEEE International Geoscience and Remote Sensing Symposium, Cape Town, South Africa, 12–17 July 2009; pp. IV-69–IV-72. 49. Rasmussen, C.E.; Williams, C.K. Gaussian Processes for Machine Learning; MIT Press: Cambridge, MA, USA, 2006; Volume 2, p. 4. 50. Zeng, Y.L.; Li, J.; Liu, Q.H.; Qu, Y.H.; Huete, A.R.; Xu, B.D.; Yin, G.F.; Zhao, J. An Optimal Sampling Design for Observing and Validating Long-Term Leaf Area Index with Temporal Variations in Spatial Heterogeneities. Remote Sens. 2015, 7, 1300–1319. [CrossRef] 51. Vermote, E.; Justice, C.; Claverie, M.; Franch, B. Preliminary analysis of the performance of the Landsat 8/OLI land surface reflectance product. Remote Sens. Environ. 2016, 185, 46–56. [CrossRef] 52. Schaaf, C.B.; Gao, F.; Strahler, A.H.; Lucht, W.; Li, X.W.; Tsang, T.; Strugnell, N.C.; Zhang, X.Y.; Jin, Y.F.; Muller, J.P.; et al. First operational BRDF, albedo nadir reflectance products from MODIS. Remote Sens. Environ. 2002, 83, 135–148. [CrossRef] 53. Chen, J.; Jonsson, P.; Tamura, M.; Gu, Z.H.; Matsushita, B.; Eklundh, L. A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky-Golay filter. Remote Sens. Environ. 2004, 91, 332–344. [CrossRef] 54. 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. [CrossRef] 55. Verrelst, J.; Rivera, J.P.; Veroustraete, F.; Munoz-Mari, 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. 2015, 108, 260–272. [CrossRef] 56. Xie, Q.; Huang, W.; Zhang, B.; Chen, P.; Song, X.; Pascucci, S.; Pignatti, S.; Laneve, G.; Dong, Y. Estimating winter wheat leaf area index from ground and hyperspectral observations using vegetation indices. IEEE J. Sel. Top. Appl. Earth Obs. Remote. Sens. 2016, 9, 771–780. [CrossRef] 57. Huete, A.; Didan, K.; Miura, T.; Rodriguez, E.P.; Gao, X.; Ferreira, L.G. Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote Sens. Environ. 2002, 83, 195–213. [CrossRef] 58. Croft, H.; Chen, J.M.; Luo, X.Z.; Bartlett, P.; Chen, B.; Staebler, R.M. Leaf chlorophyll content as a proxy for leaf photosynthetic capacity. Glob. Chang. Biol. 2017, 23, 3513–3524. [CrossRef] [PubMed] 59. Morisette, J.T.; Baret, F.; Privette, J.L.; Myneni, R.B.; Nickeson, J.E.; Garrigues, S.; Shabanov, N.V.; Weiss, M.; Fernandes, R.A.; Leblanc, S.G.; et al. Validation of global moderate-resolution LAI products: A framework proposed within the CEOS Land Product Validation subgroup. IEEE Trans. Geosci. Remote Sens. 2006, 44, 1804–1817. [CrossRef] 60. Jiang, C.Y.; Ryu, Y.; Fang, H.L.; Myneni, R.; Claverie, M.; Zhu, Z.C. Inconsistencies of interannual variability and trends in long-term satellite leaf area index products. Glob. Chang. Biol. 2017, 23, 4133–4146. [CrossRef] [PubMed] 61. Smettem, K.R.J.; Waring, R.H.; Callow, J.N.; Wilson, M.; Mu, Q.Z. Satellite-derived estimates of forest leaf area index in southwest Western Australia are not tightly coupled to interannual variations in rainfall: Implications for groundwater decline in a drying climate. Glob. Chang. Biol. 2013, 19, 2401–2412. [CrossRef] [PubMed] 62. Busetto, L.; Casteleyn, S.; Granell, C.; Pepe, M.; Barbieri, M.; Campos-Taberner, M.; Casa, R.; Collivignarelli, F.; Confalonieri, R.; Crema, A.; et al. Downstream services for rice crop monitoring in Europe: From regional to local scale. IEEE J. Sel. Top. Appl Earth Obs. Remote. Sens. 2017, 10, 5423–5441. [CrossRef] Remote Sens. 2019, 11, 244 18 of 18 63. Peng, J.; Blessing, S.; Giering, R.; Muller, B.; Gobron, N.; Nightingale, J.; Boersma, F.; Muller, J.P. Quality-assured long-term satellite-based leaf area index product. Glob. Chang. Biol. 2017, 23, 5027–5028. [CrossRef] 64. Sabater, J.M.; Rudiger, C.; Calvet, J.C.; Fritz, N.; Jarlan, L.; Kerr, Y. Joint assimilation of surface soil moisture and LAI observations into a land surface model. Agric. For. Meteorol. 2008, 148, 1362–1373. [CrossRef] 65. Tillack, A.; Clasen, A.; Kleinschmit, B.; Foerster, M. Estimation of the seasonal leaf area index in an alluvial forest using high-resolution satellite-based vegetation indices. Remote Sens. Environ. 2014, 141, 52–63. [CrossRef] 66. Yin, G.F.; Li, A.N.; Zeng, Y.L.; Xu, B.D.; Zhao, W.; Nan, X.; Jin, H.A.; Bian, J.H. A cost-constrained sampling strategy in support of LAI product validation in mountainous areas. Remote Sens. 2016, 8, 704. [CrossRef] 67. Xu, B.D.; Li, J.; Park, T.; Liu, Q.H.; Zeng, Y.L.; Yin, G.F.; Zhao, J.; Fan, W.L.; Yang, L.; Knyazikhin, Y.; et al. An integrated method for validating long-term leaf area index products using global networks of site-based measurements. Remote Sens. Environ. 2018, 209, 134–151. [CrossRef] 68. Lazaro-Gredilla, M.; Titsias, M.K.; Verrelst, J.; Camps-Valls, G. Retrieval of biophysical parameters with heteroscedastic Gaussian processes. IEEE Geosci. Remote Sens. Lett. 2014, 11, 838–842. [CrossRef] 69. Wang, Z. Mapping foliar functional traits and their uncertainties across three years in a grassland experiment. Remote Sens. Environ. 2018. [CrossRef] 70. Somarathna, P.; Minasny, B.; Malone, B.P.; Stockmann, U.; McBratney, A.B. Accounting for the measurement error of spectroscopically inferred soil carbon data for improved precision of spatial predictions. Sci. Total Environ. 2018, 631–632, 377–389. [CrossRef] 71. Greaves, H.E. High-resolution mapping of aboveground shrub biomass in Arctic tundra using airborne lidar and imagery. Remote Sens. Environ. 2016. [CrossRef] 72. Viscarra Rossel, R.A.; Webster, R.; Bui, E.N.; Baldock, J.A. Baseline map of organic carbon in Australian soil to support national carbon accounting and monitoring under climate change. Glob. Chang. Biol 2014, 20, 2953–2970. [CrossRef] © 2019 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 (http://creativecommons.org/licenses/by/4.0/).