Retrieval of Leaf Area Index (LAI) and Fraction of Absorbed Photosynthetically Active Radiation (FAPAR) from VIIRS Time-Series Data

Long-term high-quality global leaf area index (LAI) and fraction of absorbed photosynthetically active radiation (FAPAR) products are urgently needed for the study of global change, climate modeling, and many other problems. As the successor of the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor, the Visible Infrared Imaging Radiometer Suite (VIIRS) will continue to provide global environmental measurements. This paper aims to generate longer time series Global LAnd Surface Satellite (GLASS) LAI and FAPAR products after the era of the MODIS sensor. To ensure spatial and temporal consistencies between GLASS LAI/FAPAR values retrieved from different satellite observations, the GLASS LAI/FAPAR retrieval algorithms were adapted in this study to retrieve LAI and FAPAR values from VIIRS surface reflectance time-series data. After reprocessing of the VIIRS surface reflectance to remove remaining effects of cloud contamination and other factors, a database generated from the GLASS LAI product and the reprocessed VIIRS surface reflectance for all Benchmark Land Multisite Analysis and Intercomparison of Products (BELMANIP) sites was used to train general regression neural networks (GRNNs). The reprocessed VIIRS surface reflectance data from an entire year were entered into the trained GRNNs to estimate the one-year LAI values, which were then used to calculate FAPAR values. A cross-comparison indicates that the LAI and FAPAR values retrieved from VIIRS surface reflectance were generally consistent with the GLASS, MODIS and Geoland2/BioPar version 1 (GEOV1) LAI/FAPAR values in their spatial patterns. The LAI/FAPAR values retrieved from VIIRS surface reflectance achieved good agreement with the GLASS LAI/FAPAR values (R2 = 0.8972 and RMSE = 0.3054; and R2 = 0.9067 and RMSE = 0.0529, respectively). However, validation of the LAI and FAPAR values derived from VIIRS reflectance data is now limited by the scarcity of LAI/FAPAR ground measurements.


Introduction
Leaf area index (LAI) is defined as the one-sided green leaf area per unit ground area in broadleaf canopies and as one-half the total needle surface area per unit ground area in coniferous canopies [1]. The fraction of absorbed photosynthetically active radiation (FAPAR) is defined as the fraction of incident photosynthetically active radiation (400-700 nm) absorbed by the green elements of a vegetation canopy [2]. LAI and FAPAR are widely used in agriculture and ecology and are two of the essential variables identified by the Global Climate Observing System (GCOS) [3]. Estimating LAI and FAPAR from satellite observations is the most useful way to generate LAI/FAPAR products at regional and global scales.
Many algorithms have been developed to retrieve LAI or FAPAR from satellite remote-sensing data [4]. In general, two types of methods can be distinguished: empirical and physical. The empirical methods are based on statistical relationships between LAI/FAPAR and spectral vegetation indices (VIs), which are calibrated using field measurements of LAI/FAPAR and reflectance data recorded by a remote sensor or simulated using canopy radiation models (e.g., Wang et al. [5]). These statistical relationships depend not only on vegetation types but also on the points in time and the geographical area the field measurements are collected at. The physical methods are based on inversion of canopy radiative transfer models [6][7][8]. Inversion techniques based on iterative minimization of a cost function require hundreds of runs of the canopy radiative-transfer model for each pixel and are therefore computationally too demanding. For practical applications, lookup table (LUT) methods [9] and machine learning methods [10] are two popular inversion techniques that are based on a pre-computed reflectance database.
Currently, multiple global LAI and FAPAR products have been produced from data acquired by the Advanced Very High Resolution Radiometer (AVHRR) [11], the Multiangle Imaging SpectroRadiometer (MISR) [12], the Moderate Resolution Imaging Spectroradiometer (MODIS) [9,13], VEGETATION [10], the Sea-Viewing Wide Field-of-View Sensor (SeaWiFS) [14], and the Medium-Resolution Imaging Spectrometer (MERIS) [15]. However, these LAI/FAPAR products are generated using different methods developed independently by different investigators under different assumptions, resulting in inconsistencies between LAI/FAPAR products [16,17]. Furthermore, currently available LAI products have limited length of the time series. In fact, a consistent long-term data set is urgently needed for the study of global change, climate modeling, and many other problems. As part of the efforts for generating the Global LAnd Surface Satellite (GLASS) products for the long-term environmental change studies [18,19], Xiao et al. [20] developed a method to retrieve LAI values from satellite observations based on general regression neural networks (GRNNs) and Xiao et al. [21] proposed a method to calculate FAPAR values from the LAI values. These methods were used to generate a long time series of GLASS LAI/FAPAR products (1981-2014) from AVHRR and MODIS reflectance data.
It is well known that the Terra and Aqua satellites equipped with the MODIS sensor had a life expectancy of six years and have far exceeded their design life, and the Visible Infrared Imaging Radiometer Suite (VIIRS) sensor on board the Suomi National Polar-Orbiting Partnership (SNPP) satellite will continue to provide global environmental measurements [22,23]. This paper aims to generate longer time series GLASS LAI and FAPAR products after the era of the MODIS sensor. To ensure spatial and temporal consistencies between GLASS LAI/FAPAR values retrieved from different satellite observations, the GLASS LAI/FAPAR retrieval algorithms were adapted in this study to retrieve LAI and FAPAR values from the VIIRS surface reflectance data. The VIIRS surface reflectance was reprocessed to remove the remaining effects of cloud contamination and other factors, and a database generated from the GLASS LAI product and the reprocessed VIIRS surface reflectance from the Benchmark Land Multisite Analysis and Intercomparison of Products (BELMANIP) sites for 2013 and 2014 were used to train GRNNs. Reprocessed VIIRS surface reflectance data for an entire year were entered into the trained GRNNs to estimate the one-year LAI values, which were then used to calculate FAPAR values. The LAI and FAPAR values retrieved from VIIRS reflectance data in this study were compared with the GLASS, MODIS, and Geoland2/BioPar version 1 (GEOV1) LAI/FAPAR products. This paper is organized as follows. Section 2 describes the VIIRS reflectance product and the existing global LAI and FAPAR products, including GLASS, MODIS, and GEOV1. Several details of the algorithm implementation, including reprocessing of the VIIRS surface reflectance, retrieval of LAI values using GRNNs, and calculation of FAPAR values from the retrieved LAI values, are described in Section 3. Comparisons of the LAI/FAPAR values retrieved from time-series VIIRS surface reflectance data with the MODIS and GEOV1 LAI/FAPAR products are presented in Section 4. The paper concludes in Section 5.

Visible Infrared Imaging Radiometer Suite (VIIRS) Surface Reflectance
As with MODIS, surface reflectance is one of the key products from VIIRS. The VIIRS surface reflectance intermediate product, including three bands (Table 1) with a 500 m spatial resolution and an eight-day temporal sampling period, is based on the heritage MODIS Collection 5 product [24]. The VIIRS surface reflectance intermediate product, in a sinusoidal projection system, is available from the LAADS (Level 1 and Atmosphere Archive and Distribution System) Web [25]. The GLASS LAI and FAPAR products are one of the longest-duration LAI and FAPAR products in the world and were generated and released by the Center for Global Change Data Processing and Analysis of Beijing Normal University [26]. They are also available from the Global Land Cover Facility [27]. The GLASS LAI and FAPAR products have a temporal resolution of eight days, and spans 1981-2014. For 2000-2014, the GLASS LAI product was provided in a sinusoidal projection at a spatial resolution of 1 km and derived from MODIS collection 5 surface reflectance data using GRNNs which were trained using fused time-series LAI values from the CYCLOPES and MODIS collection 5 LAI products and reprocessed time-series MODIS surface reflectance for 2001-2003 over the BELMANIP sites [20]. For 1981-1999, the GLASS LAI product was provided in a geographic latitude/longitude projection at a spatial resolution of 0.05˝(~5 km at the Equator) and derived from AVHRR surface reflectance data from NASA's Land Long-Term Data Record (LTDR) project [28] using GRNNs. The fused LAI values from the CYCLOPES and MODIS collection 5 LAI products were aggregated to 0.05˝spatial resolution using a spatial-average method. Then, the aggregated LAI time-series values and the corresponding reprocessed AVHRR reflectance values over the BELMANIP sites for 2003 and 2004 were used to train GRNNs [20]. Unlike existing neural network methods that use remote-sensing data acquired only at a specific time to retrieve LAI [10,29], the reprocessed MODIS/AVHRR reflectance values from an entire year were inputted to the GRNNs to estimate the one-year LAI profiles. To ensure physical consistency between LAI and FAPAR retrievals, the GLASS FAPAR product was generated from the GLASS LAI product and has the same properties as the GLASS LAI product. The GLASS FAPAR values correspond to total FAPAR at 10:30 A.M. local time, which is a close approximation of daily average FAPAR [21]. The GLASS LAI and FAPAR products are spatially complete and have continuous trajectories. Direct comparison with ground-based estimates from the Validation of Land European Remote sensing Instrument (VALERI) project [30] and the BigFoot project [31] demonstrated that the GLASS LAI values were closer to the ground-based LAI estimates (RMSE = 0.7848 and R 2 = 0.8095) than the GEOV1 LAI values (RMSE = 0.9084 and R 2 = 0.7939) and the MODIS LAI values (RMSE = 1.1173 and R 2 = 0.6705) and the GLASS FAPAR values were also more accurate (RMSE = 0.0716 and R 2 = 0.9292) than the GEOV1 FAPAR values (RMSE = 0.1085 and R 2 = 0.8681) and the MODIS FAPAR values (RMSE = 0.1276 and R 2 = 0.8048) [20,21].
The MODIS LAI and FAPAR products have been available since 2000 and are provided in a sinusoidal projection. The MODIS LAI/FAPAR retrieval algorithm includes a main algorithm and a backup algorithm. The main algorithm is based on LUTs simulated using a three-dimensional radiative-transfer model. When the main algorithm fails, the backup algorithm is used to estimate LAI and FAPAR from biome-specific relationships between LAI/FAPAR and normalized difference vegetation index (NDVI) [13]. Generally, LAI and FAPAR estimates using the backup algorithm are of lower quality, mainly because of residual clouds and poor atmospheric correction [32]. The MODIS FAPAR product is defined as the instantaneous black sky FAPAR (i.e., under direct illumination) at the time of the Terra overpass (10:30 A.M.). Until now, MODIS LAI and FAPAR products have been released in a total of six versions. In collection 5, parameters of the main and backup algorithms are defined for eight main biome classes according to MODIS land cover, and a new stochastic radiative-transfer model was used to provide a better representation of canopy structure and the spatial heterogeneity intrinsic to woody biomes. Collection 5 MODIS LAI/FAPAR products are generated at 1 km spatial resolution and an eight-day time step. Unlike MODAGAGG used in Collection 5, the collection 6 retrieval algorithm uses daily L2G´lite surface reflectance (MOD09GA) as input and uses an improved multi-year land-cover product [33]. Collection 6 MODIS LAI/FAPAR products are generated at a spatial resolution of 500 m. The temporal compositing periods are eight and four days [33]. In this study, both collection 5 and 6 LAI/FAPAR products are compared with the LAI and FAPAR retrieved from VIIRS surface reflectance data. For clarification, the collection 5 LAI and FAPAR products are denoted by MODC5, and the collection 6 LAI and FAPAR products are denoted by MODC6.
The GEOV1 LAI and FAPAR products have been available since 1998 [34]. The products are provided in a Plate Carrée projection at 1/112˝spatial resolution and 10-day frequency. The GEOV1 LAI and FAPAR products were derived from SPOT/VEGETATION sensor data using back-propagation neural networks. The MODIS and CYCLOPES LAI/FAPAR products were fused and scaled to generate best estimates of LAI and FAPAR, which were used to train the back-propagation neural networks using the SPOT/VEGETATION top-of-canopy nadir reflectance values over the BELMANIP sites [29]. The calibrated neural networks were used to generate the GEOV1 LAI and FAPAR products from SPOT/VEGETATION top-of-canopy nadir reflectance data. The GEOV1 FAPAR product corresponds to the instantaneous black sky FAPAR by green parts at 10:15 A.M. local time.

Experimental Area
Tile h09v05 was selected to investigate spatial patterns specific to the LAI/FAPAR values retrieved from the VIIRS surface reflectance data. It covers Texas and New Mexico, USA. The main biome types of this tile include grasses/cereal crops, shrubs, savannah, needleleaf forests, and broadleaf forests according to the MODIS land-cover map (MCD12Q1).
In addition, six sites with different biome types were selected to compare the time series of the reconstructed VIIRS and MODIS surface reflectance and the time series of LAI and FAPAR values. The six sites we chose were: Konza, Argo, Counami, Larose, Laprida, Donga, and their attributes are shown in Table 2.

Reprocessing of the VIIRS Surface Reflectance Data
VIIRS surface reflectance was estimated from satellite observations using atmospheric correction algorithm which consists of aerosol LUTs that are populated by the Second Simulation of the Satellite Signal in the Solar Spectrum (6S) radiative transfer model [35]. Although great efforts have been made, the VIIRS surface reflectance product still contains considerable noise caused, for instance, by cloud or mixed-cloud pixels, resulting in temporal and spatial inconsistencies in subsequent downstream products.
To remove residual clouds from the land-surface reflectance product, Xiao et al. [36] developed a method (referred to hereafter as VIRR algorithm for clarification) which incorporates upper envelopes of time series VIs as constraint conditions to reconstruct time series of surface reflectance for the red, near-infrared (NIR), and shortwave infrared (SWIR) bands. Satellite-retrieved surface reflectance data are used to calculate VIs, and a penalized least square regression based on three-dimensional discrete cosine transform (DCT-PLS) [37] is used to calculate continuous and smooth upper envelopes of VIs. Cloud-contaminated surface reflectance values were detected using the time-series VIs and the upper envelopes of the time-series VIs. Surface reflectance data with good quality in a given time window, along with the continuous and smooth upper envelopes of VIs, were used to estimate the optimal values of coefficients of quadratic polynomial functions fitted to the surface reflectance data. Then, time series of surface reflectance can be reconstructed using the quadratic polynomial functions according to the optimal values of the coefficients. If there are large gaps because of either cloud cover or missing observations, the reflectance values were filled using those in preceding or succeeding years according to better quality. Detailed information about the VIRR method can be found in Xiao et al. [36]. The method was used to reconstruct time series of surface reflectance from the MODIS/TERRA surface reflectance product (MOD09A1). The cloud-free MODIS/AQUA surface reflectance product (MYD09A1) was used as reference data for evaluation of the reconstructed time series of surface reflectance. Comparisons of the collocated surface reflectance values from MYD09A1 and the reconstructed surface reflectance over the BELMANIP sites in 2003 demonstrate that the VIRR method provides good performance for the red band (R 2 = 0.8606 and RMSE = 0.0366) and NIR band (R 2 = 0.6934 and RMSE = 0.0519) [36].
As the successor of the MODIS sensor, the VIIRS sensor has similar spectral bands with the MODIS sensor. In this study, the VIRR method was used to identify cloud-contaminated pixels in the VIIRS surface reflectance data and to generate continuous surface reflectances from VIIRS. The reconstructed surface reflectance values were then used to retrieve LAI and FAPAR.

LAI Retrieval Using General Regression Neural Networks (GRNNs)
In this study, the GLASS LAI retrieval algorithm was adapted to retrieve LAI values consistent with the GLASS LAI product from VIIRS reflectance time-series data using GRNNs. In order to retrieve global LAI from VIIRS surface reflectance using the GRNNs, the training database should be broadly representative. The BELMANIP sites provide a good sampling of biome types and conditions throughout the world [38]. There are a total of 402 BELMANIP sites [39]. A database for the GRNNs was generated over all BELMANIP sites. For each BELMANIP site, 7ˆ7 pixels of the GLASS LAI product and the reprocessed VIIRS reflectance data from 2013-2014 were extracted. The database was randomly split into a training dataset (85% of the data) to train the GRNNs and a testing dataset (15% of the data) to evaluate their performance. To achieve better performance and convergence, the range of variation of input and output values were normalized between´1 and 1 according to the minimum and maximum values of the input and output values respectively. The GRNN, developed by Specht [40], is a generalization of radial basis function networks and probabilistic neural networks. The advantage of this type of neural network is that it can approximate the map inherent in any sample data set. In addition, the GRNN does not require iterative training; the functional estimate is computed directly from the training data. If the kernel function of a GRNN is Gaussian, the fundamental formulation of the GRNN can be expressed as follows: where D 2 i "`X´X i˘T`X´Xi˘r epresents the squared Euclidean distance between the input vector X and the ith training input vector X i , Y i is the output vector corresponding to the vector X i , Y 1 pXq is the estimate corresponding to the vector X, n is the number of samples, and σ is a smoothing parameter that controls the size of the receptor region. Equation (1) shows that the estimate Y 1 pXq, given an input vector X, is the weighted average of all the sample observations Y i , where the weight for each observation is proportional to the Euclidean distance between the vector X and the training input vector X i . When the input to the GRNN is given, the architecture and weights of the GRNN are determined.
To retrieve LAI profiles from VIIRS surface reflectance time-series data using a GRNN, the input vector X of the GRNN is the reprocessed VIIRS time-series reflectance values (for a one-year period); and contains 46ˆm components, where m is the number of bands. The output vector Y 1 " pLAI 1 , LAI 2 ,¨¨¨, LAI 46 q T is the corresponding LAI time series for the year and contains 46 components.
The smoothing parameter σ in Equation (1) is the only free parameter in the GRNN formulation. Therefore, GRNN training essentially involves optimizing the smoothing parameter. In this study, the holdout method was used to find the optimal smoothing parameter at which the following cost function reaches its minimum: whereŶ i pX i q is the estimate corresponding to X i using the GRNN trained over all the training samples except the ith sample. The performance of the trained GRNNs was evaluated over the testing dataset. Figure 1 shows the density scatterplots between the testing LAI values and the LAI values retrieved from VIIRS surface reflectance in the red, NIR, and SWIR bands using the GRNNs. It shows that the training was very efficient for the LAI retrieval. Most points in Figure 1 are around the 1:1 line, which indicates that the retrieved LAI values using the GRNNs achieve excellent agreement with the testing LAI values across the LAI range (R 2 = 0.9478, RMSE = 0.3891 and Bias =´0.0184).
The trained GRNNs were then used to retrieve LAI from the reprocessed VIIRS reflectance data. The reprocessed VIIRS reflectance data for a one-year period were entered into the trained GRNNs, and the output of the trained GRNNs represented the one-year LAI profile for each pixel. The performance of the trained GRNNs was evaluated over the testing dataset. Figure 1 shows the density scatterplots between the testing LAI values and the LAI values retrieved from VIIRS surface reflectance in the red, NIR, and SWIR bands using the GRNNs. It shows that the training was very efficient for the LAI retrieval. Most points in Figure 1 are around the 1:1 line, which indicates that the retrieved LAI values using the GRNNs achieve excellent agreement with the testing LAI values across the LAI range (R 2 = 0.9478, RMSE = 0.3891 and Bias = -0.0184).

Calculation of FAPAR
In this study, the GLASS FAPAR retrieval algorithm was adapted to calculate FAPAR values consistent with the GLASS FAPAR product from the LAI values retrieved from VIIRS surface reflectance. Let τ PAR be the transmittance of PAR down to the soil; FAPAR can then be approximately calculated as follows [34]: FAPAR " 1´τ PAR .
The radiation into the vegetation canopy includes direct and diffuse PAR. Therefore, the transmittance of PAR down to the soil can also be expressed as: where τ dir PAR and τ dif PAR are the fractions of the radiative flux originating from the direct illumination source and the transmitted fraction of the incident diffuse illumination source respectively, and f skyl is the fraction of diffuse skylight. If the leaf area index of a canopy is lai and the absorptivity of leaves for radiation is a, the fraction of the total beam radiation transmitted through a canopy can be approximated using an exponential model [41]: where Ω is the clumping index, ϕ is the solar zenith angle, and k c pϕq is the canopy extinction coefficient for PAR. For an ellipsoidal leaf angle distribution, k c pϕq is calculated as follows [41]. where x is the ratio of average projected areas of canopy elements on horizontal and vertical surfaces. Diffuse radiation comes from all directions. Therefore, the diffuse transmission coefficient can be calculated by integrating the direct transmission coefficient over all illumination directions: In the method just described, LAI is an important input parameter to estimate FAPAR. The LAI values retrieved from VIIRS surface reflectance were used to calculate FAPAR values at 10:30 A.M. local time, which are close approximations of daily average FAPAR. The clumping index is another input parameter. He et al. [42] used the MODIS bidirectional reflectance distribution function (BRDF) parameter to derive a global clumping index map at 500 m resolution. In the present study, the MODIS-derived clumping index map was used to calculate canopy transmittance.

Comparison and Analysis
The reconstructed VIIRS surface reflectance data at the AGRO site are presented to illustrate the performance of the VIRR algorithm. The similarities between the time series of the reconstructed VIIRS and MODIS surface reflectance over the selected sites were quantitatively evaluated. In this study, the following agreement index (S AI ) [43,44] was adopted to measure similarity.
where x i and y i are the time series values at time i, respectively, x is the mean value of x i , M is the length of the time series. The similarity ranges from 0.0% to 100.0%, where 0.0% represents complete disagreement between the time series x i and y i while 100.0% indicates that the time series x i and y i are identical. The VIIRS surface reflectance data, including three bands (red, NIR, and SWIR), were used to retrieve LAI values. For clarification, the LAI and FAPAR products derived in this study are denoted as VIIRS products. The VIIRS LAI and FAPAR products have eight-day temporal and 1 km spatial resolution and are provided in a sinusoidal projection. A cross-comparison was performed to evaluate the spatial and temporal consistencies between the VIIRS LAI/FAPAR products and other existing global LAI/FAPAR products, including GLASS, MODC5, MODC6, and GEOV1 LAI/FAPAR products. For comparisons of spatial consistency, the maps of these LAI and FAPAR products over MODIS tile h09v05 were compared to investigate spatial patterns as well to check the distribution in space of the missing data, and density scatterplots among these LAI/FAPAR products over the tile were also analyzed. To compare temporal consistency, temporal profiles of the differences between GLASS, MODC5, MODC6 LAI/FAPAR products and VIIRS LAI/FAPAR products were checked over the selected sites. The similarities between the series of VIIRS LAI/FAPAR values and GLASS, MODC5, and MODC6 LAI/FAPAR values were calculated to evaluate the agreement between these LAI/FAPAR time series.
For the evaluations of spatial and temporal consistencies, the GEOV1 LAI and FAPAR products were re-projected onto the sinusoidal projection used in the MODIS and GLASS LAI/FAPAR products using the General Cartographic Transformation Package (GCTP) map projection library [45] and resampled to 1 km spatial resolution using the nearest-neighbor resampling method in order not to alter the product values.
In this study, only valid LAI and FAPAR values of these products were used for comparison. For the VIIRS, GLASS and GEOV1 LAI/FAPAR products, all LAI and FAPAR values were considered to be valid. For the MODC5 and MODC6 LAI/FAPAR products, the LAI and FAPAR values retrieved from the backup algorithm were not used for comparison of the LAI and FAPAR products because of their overall lower quality originating from residual clouds and poor atmospheric correction [32]. In other words, only the LAI and FAPAR values retrieved from the main algorithm (QC < 64) were considered to be valid.  Figure 2. The VIIRS surface reflectance, especially in the red and NIR bands, has several abnormal values at this site. The VIRR algorithm identified those observations with high reflectance values as clouds, and the residual clouds were removed in the reconstructed surface reflectance. It is apparent that the reconstructed surface reflectance from the VIRR method is in good agreement with the cloud-free VIIRS surface reflectance. To further evaluate the reconstructed VIIRS surface reflectance, the similarities between the series of the reconstructed VIIRS and MODIS surface reflectance data were calculated ( Table 3). The time-series similarities are high, especially at the Argo and Konza sites where the similarities for the red, NIR and SWIR reflectance are larger than 88%. However, the similarities between the series of the reconstructed VIIRS and MODIS surface reflectance especially in the red and SWIR bands are relatively poor at the Counami site. The similarity for the red reflectance is less than 30%. The low similarity for the Counami site is due to the fact that the number of high-quality satellite observations was very limited in tropical rain forests.   To further evaluate the reconstructed VIIRS surface reflectance, the similarities between the series of the reconstructed VIIRS and MODIS surface reflectance data were calculated ( Table 3).

Evaluation of the Reconstructed VIIRS Surface Reflectance
The time-series similarities are high, especially at the Argo and Konza sites where the similarities for the red, NIR and SWIR reflectance are larger than 88%. However, the similarities between the series of the reconstructed VIIRS and MODIS surface reflectance especially in the red and SWIR bands are relatively poor at the Counami site. The similarity for the red reflectance is less than 30%. The low similarity for the Counami site is due to the fact that the number of high-quality satellite observations was very limited in tropical rain forests.

Spatial Comparison
An example of VIIRS, MODC5, MODC6, GLASS, and GEAV1 LAI products for tile h09v05 on day 185, 2013 is shown in Figure 3. Areas masked in grey correspond to pixels where these LAI products did not provide valid LAI values. As in the GLASS LAI product, no missing values exist in the VIIRS LAI product. It is apparent that the VIIRS, MODC5, MODC6, GLASS, and GEAV1 LAI products are generally consistent in their spatial patterns, depicting a clear correlation with the MODIS biome map (Figure 3f). These LAI products took on their highest values over broadleaf forest regions. LAI values were intermediate over grass/cereal crop, needleleaf forest and savannah regions and were very low over shrub regions. However, discrepancies existed in the magnitudes of these LAI products. Over the broadleaf forest regions, the MODC5 LAI values could reach 6.9 LAI units and were larger than the corresponding VIIRS, GLASS and GEOV1 LAI values, which was partly due to overestimation of the MODC5 LAI values associated with broadleaf forests [16]. Obviously, the MODC6 LAI values have significant improvements. Over needleleaf forest regions, the VIIRS, GLASS and GEOV1 LAI products ranged between 3.0 and 5.0 LAI units and were clearly higher than the MODC5 and MODC6 LAI values (around 1.8 LAI units). Figure 4 shows density scatterplots among the VIIRS, GLASS, MODC5, MODC6, and GEOV1 LAI values for tile h09v05 on day 185, 2013. Only the collocated LAI values among these LAI products are included in Figure 4. Figure 4c,d illustrates that the VIIRS LAI values were slightly higher than those of MODC5 and MODC6 for lower LAI values, but clearly lower than those of MODC5, MODC6, and GEOV1 when the LAI values were greater than 3.5 LAI units. Compared with MODC5, MODC6, and GEOV1 LAI values, those of GLASS are distributed more closely around the 1:1 line against the VIIRS LAI values, which showed that the VIIRS LAI values achieved better agreement with the GLASS LAI values (R 2 = 0.8972 and RMSE = 0.3054) across the LAI range than did the other products. This can be attributed to the same GRNN structure for the GLASS and VIIRS LAI retrieval, and the LAI training data from the GLASS LAI product.       Figure 6. Like the VIIRS LAI product, the VIIRS FAPAR values achieved good agreement with the GLASS FAPAR values (R 2 = 0.9067 and RMSE = 0.0529) across the FAPAR range. However, most points in Figure 6b-d is above the 1:1 line, which indicates that the VIIRS FAPAR values are lower than those from MODC5, MODC6, and GEOV1. Previous observations indicate that the MODC5 and GEOV1 FAPAR values are also slightly larger than those of GLASS [21]. These discrepancies can be partly explained by overestimation in the MODC5, MODC6, and GEOV1 values [17,21].  0.0529) across the FAPAR range. However, most points in Figure 6b-d is above the 1:1 line, which indicates that the VIIRS FAPAR values are lower than those from MODC5, MODC6, and GEOV1. Previous observations indicate that the MODC5 and GEOV1 FAPAR values are also slightly larger than those of GLASS [21]. These discrepancies can be partly explained by overestimation in the MODC5, MODC6, and GEOV1 values [17,21].

Temporal Analysis
Temporal consistency among the VIIRS, GLASS, MODC5, and MODC6 LAI and FAPAR products was evaluated over the selected sites. Figure 7 shows temporal profiles of the differences between GLASS, MODC5, MODC6 LAI/FAPAR products and VIIRS LAI/FAPAR products, and the similarities between the series of VIIRS LAI/FAPAR values and GLASS, MODC5, and MODC6 LAI/FAPAR values were shown in Table 3. Figure 7a shows the temporal trajectories of LAI/FAPAR differences for the Konza site with grass and cereal crop biome type. There is an excellent agreement between VIIRS and GLASS LAI values during the whole year, with maximum differences of only 0.26 LAI units. The MODC6 LAI values were larger than VIIRS LAI values only from Julian days 153 to 185, and the differences

Temporal Analysis
Temporal consistency among the VIIRS, GLASS, MODC5, and MODC6 LAI and FAPAR products was evaluated over the selected sites. Figure 7 shows temporal profiles of the differences between GLASS, MODC5, MODC6 LAI/FAPAR products and VIIRS LAI/FAPAR products, and the similarities between the series of VIIRS LAI/FAPAR values and GLASS, MODC5, and MODC6 LAI/FAPAR values were shown in Table 3. Figure 7a shows the temporal trajectories of LAI/FAPAR differences for the Konza site with grass and cereal crop biome type. There is an excellent agreement between VIIRS and GLASS LAI values during the whole year, with maximum differences of only 0. 26  between MODC5 and VIIRS LAI values were negative during the whole year. Just like VIIRS LAI values, the GLASS FAPAR values were higher than those of VIIRS only from Julian days 153 to 209. However, most of the differences between MODC5, MODC6 and VIIRS FAPAR values were positive for this year. The similarities between the series of VIIRS and GLASS LAI/FAPAR values (99.1311 and 98.2412, respectively) are slightly higher than those between the series of VIIRS and MODC5, MODC6 LAI/FAPAR values at this site.  The temporal profiles of LAI/FAPAR differences for the Argo site are shown in Figure 7b. The  The temporal profiles of LAI and FAPAR differences for the Counami site with broadleaf forest biome type are shown in Figure 7c. Some MODC5 and MODC6 LAI/FAPAR values are missing due to residual cloud contamination in MODIS surface reflectance data. The MODC5 and MODC6 LAI values, except for those contaminated by residual cloud, were higher than those of VIIRS (up to 2.1 LAI unit) after day 201 for this year. At this site, the VIIRS LAI values were larger than those of GLASS (up to 1.0 LAI unit) and the VIIRS LAI values were higher than those of GLASS (up to 0.08 units) except for Julian days 265 to 297. The similarities between the series of VIIRS and GLASS LAI/FAPAR products at the Counami site are less than 50% due to the fact that the number of good observations was very limited in tropical rain forests. Nevertheless, the similarities between the series of VIIRS and GLASS LAI/FAPAR products are still larger than those between the series of VIIRS and MODC5, MODC6 LAI/FAPAR products. Figure 7d shows the temporal trajectories of LAI and FAPAR differences for the Larose site with needleleaf forest biome type. The differences demonstrate that the MODC6 LAI/FAPAR values were in better agreement with the VIIRS LAI/FAPAR values than with the MODC5 LAI/FAPAR values. The MODC5 LAI values were between 1.0 and 3.1 LAI units higher than those of VIIRS especially during the vegetation growing season. The MODC5 and MODC6 FAPAR values were in good agreement with the VIIRS FAPAR values during the vegetation growing season. However, large discrepancies between these FAPAR products could be observed during the nongrowing season, when the MODC5 FAPAR values were between 0.2 and 0.4 units higher than those of VIIRS.
The temporal profiles of LAI and FAPAR differences for the Laprida sites with savannah biome type are shown in Figure 7e. The differences between the GLASS and VIIRS LAI values were lower than 0.11 LAI units and the differences between the GLASS and VIIRS FAPAR values were less than 0.03 units for the entire year, which demonstrated that VIIRS LAI/FAPAR values achieved very good agreement with those of GLASS at this site. The VIIRS LAI values were higher than those of MODC5 and MODC6 (up to 1.0 LAI units), whereas the VIIRS FAPAR values were lower than the MODC5 and MODC6 FAPAR values (up to 0.2 units) throughout the entire year.
The temporal trajectories of LAI and FAPAR differences for the Donga site with shrub biome type in 2014 are shown in Figure 7f. The VIIRS LAI values achieved good agreement with those of MODC5 and MODC6 during the nongrowing season. The MODC5 and MODC6 FAPAR values, except for those contaminated by residual cloud, were higher than the VIIRS FAPAR values (up to 0.2 units) during the whole year. The similarities between the series of VIIRS and GLASS LAI/FAPAR products were higher than 98%, which indicate that the VIIRS LAI/FAPAR values agree well with those of GLASS.
The similarities between the series of these LAI/FAPAR products in Table 3 demonstrate that the VIIRS and GLASS LAI/FAPAR products have the best agreement among these LAI/FAPAR products. The similarities between the series of the VIIRS and GLASS LAI/FAPAR values for each site are larger than those between the series of the VIIRS and MODC5, MODC6 LAI/FAPAR values. The similarities for LAI and surface reflectance in Table 3 also demonstrate that reconstruction of surface reflectance using the VIRR method and the proposed GRNN approach in this study are an important factor to ensure the high agreement between the VIIRS and GLASS LAI/FAPAR values. One the one hand, at the Argo and Konza sites, the series of the VIIRS and GLASS LAI/FAPAR values have large similarities when the similarities between the series of the reconstructed VIIRS and MODIS surface reflectance data are large. At the Counami site, the similarities between the series of the reconstructed VIIRS and MODIS surface reflectance for the red, NIR and SWIR bands are lower than those at other sites (no more than 41%), and the similarities between the series of the VIIRS and GLASS LAI/FAPAR values at this site are also lower than those at other sites (no more than 50%). Therefore, reconstruction of surface reflectance using the VIRR method is an important factor that contributes to the high agreement between the VIIRS and GLASS LAI/FAPAR values. On the other hand, the similarity between the GLASS and VIIRS LAI time series can still reach 49.92% although the similarities between the series of the reconstructed VIIRS and MODIS surface reflectance data especially in the red and SWIR bands are small at the Counami site. What is more, the similarities between the GLASS and VIIRS LAI/FAPAR time series are larger than the similarities between the series of the reconstructed VIIRS and MODIS surface reflectance data for each site. Therefore, another important factor to ensure the high agreement between the VIIRS and GLASS LAI/FAPAR values is the proposed GRNN approach in this study which uses the same GRNN structure as the GLASS LAI retrieval method and uses GLASS LAI values as the LAI training data of the GRNNs for VIIRS LAI retrieval.

Conclusions
This paper contributes to the research community for data continuity of LAI/FAPAR (leaf area index/fraction of absorbed photosynthetically active radiation) across Advanced Very High Resolution Radiometer (AVHRR), Moderate Resolution Imaging Spectroradiometer (MODIS), and Visible Infrared Imaging Radiometer Suite (VIIRS) sensors. To generate Global LAnd Surface Satellite (GLASS) LAI and FAPAR products after the era of the MODIS sensor, the GLASS LAI/FAPAR retrieval algorithms were adapted in this study to retrieve LAI and FAPAR values from VIIRS surface reflectance time-series data. After reprocessing of the VIIRS surface reflectance to remove remaining effects of cloud contamination and other factors, a database generated from the GLASS LAI product and the reprocessed VIIRS surface reflectance for all Benchmark Land Multisite Analysis and Intercomparison of Products (BELMANIP) sites was used to train general regression neural networks (GRNNs). The reprocessed VIIRS surface reflectance data from an entire year were entered into the trained GRNNs to estimate the one-year LAI values, which were then used to calculate FAPAR values.
A cross-comparison indicates that the VIIRS LAI/FAPAR values were generally consistent with the GLASS, MODIS and GEOV1 LAI/FAPAR values in their spatial patterns. The VIIRS and GLASS LAI/FAPAR products have the best agreement (R 2 = 0.8972 and RMSE = 0.3054; and R 2 = 0.9067 and RMSE = 0.0529, respectively) among these LAI/FAPAR products. The high agreement between the VIIRS and GLASS LAI/FAPAR values is partly attributed to the proposed GRNN approach which uses the same GRNN structure as the GLASS LAI retrieval method and uses GLASS LAI values as the LAI training data of the GRNNs for VIIRS LAI retrieval. Reconstruction of surface reflectance using the VIRR method is another factor to ensure the high agreement between the VIIRS and GLASS LAI/FAPAR values. Therefore, the method proposed in this study can be used to generate longer time-series GLASS LAI and FAPAR products from VIIRS surface reflectance data.
However, the proposed GRNN approach requires that the input of the GRNNs is one year of surface reflectance without missing values. Therefore, the VIRR method must be used to fill the missing surface reflectance values before deriving LAI/FPAR. In addition, validation of VIIRS LAI and FAPAR products is now limited by the scarcity of LAI/FAPAR ground measurements. In the near future, the authors hope to carry out more extensive validation and analysis of the LAI and FAPAR values derived from VIIRS reflectance data by searching for ground measurements with which to perform direct validation.