On Line Validation Exercise (OLIVE): A Web Based Service for the Validation of Medium Resolution Land Products. Application to FAPAR Products

3.2.1. GEOLAND Version1 (GEOV1)

GEOV1 [47] is derived from SPOT/VEGETATION sensors with a 10-day temporal sampling and 1/112° (about 1km at equator) ground sampling distance, in a Plate Carrée projection for the period 1998 (December) up to present [48].

Algorithm inputs are the same as for CYCLOPES [3], i.e., atmospherically corrected red, near-infrared and shortwave-infrared reflectances; normalized to a standard view-illumination geometry (nadir viewing, sun zenith angle corresponding to the median value in the compositing period) and with automatic outlier rejection [49]. The directional normalization is achieved by a parametric model inversion over data acquired during a 30 days compositing window period shifted every 10 days [50]. GEOLAND version 1 products are estimated using a neural network trained on values issued from the fusion of CYCLOPES [3] and MODIS Collection 5 products [6] to take advantage of their specific performances while limiting the situations where they show deficiencies. The learning data base was composed of the BELMANIP2.1 extracts of CYCLOPES L3A normalized reflectance products and the weighted sum of corresponding CYCLOPES and MODIS FAPAR products. FAPAR products correspond to the black sky values at 10:00 solar time.

A quality assessment flag value indicates when the retrieval algorithm fails (QA = 1), i.e., when there are less than two cloud- and snow-free observations in the compositing period (30 days) or when the retrieval is out of the valid range.

3.2.2. JRC-FAPAR (MGVI₀)

MGVI (MERIS Global Vegetation Index) is a FAPAR product derived from MERIS sensor with a 10-day temporal sampling for the period (June) 2002 up to (April) 2012 [10,51]. Algorithm inputs are the daily L1b MERIS top of atmosphere (TOA) reflectances in the blue (band 2), red (band 8), near-infrared (band 14). MGVI is computed from the formulae provided in [10]. The 10 days synthesis is performed according to [52] and consists in selecting the MGVI value that is the closest to the average valid value over the 10 day synthesis. Valid values are identified using the MERIS L1b flag. FAPAR products correspond to the instantaneous black sky value at the time of satellite observations (about 10:30 solar time).For practical considerations, MERIS L1b TOA reflectances for BELMANIP2.1 and DIRECT sites were extracted and projected in the UTM (WGS84) at 1 km resolution while the MGVI computation and synthesis was performed by EMMAH, INRA. Note that the results obtained were checked with JRC to be compliant with the official MERIS MGVI instantaneous product value. However, we were not able to check the compositing algorithm because of the unavailability of time series extracts over BELMANIP2.1 or DIRECT sites. For this reason, this FAPAR product will be named as MGVIo in the following of the study as they do not correspond to the official product. We also associate a quality assessment flag (QA) value indicated when the retrieval algorithm fails, i.e., when no valid value is available within the 10 day synthesis. Therefore, for MGVI, QA = 0 indicates that the data was valid and QA = 1 indicates that the data was not produced.

3.2.3. Terra MODIS FPAR Collection 5 Product (MODC5)

MODC5 products [6] are available at WWW4. Collection 5 was produced between February 2000 through present at an 8-day time step and 1km resolution over an integrized sinusoidal grid.

MODIS main algorithm is based on a three-dimensional radiative transfer canopy model inversion based on Look Up Tables (LUT) [53,54]. Red and near-infrared atmospherically corrected MODIS reflectance [55,56] and the corresponding illumination-view geometry are used as inputs. The output is the mean FAPAR (and LAI) values averaged over the set of acceptable solutions for which simulated MODIS surface reflectances agree within specified level of model and measurement uncertainties. When the main algorithm fails, a backup algorithm based on NDVI (Normalized Difference Vegetation Index) relationships, calibrated over the same radiative transfer model simulations is used. MODIS FAPAR retrieval is executed irrespective of cloud and snow state but the majority of retrievals under these conditions or with residual atmospheric contamination are performed with the back-up algorithm [57,58].

The main and backup algorithms are defined for 6 biome types as considered by the MODIS land cover map product [59,60]. MODIS algorithm accounts for vegetation clumping at the canopy and leaf (shoot) scales through 3D radiative transfer formulations. The clumping at the landscape scale is partly addressed via mechanisms based on radiative transfer theory of canopy spectral invariants [53,61,62].

MODISC5 FAPAR product corresponds to the instantaneous value at the time of the satellite overpass. FAPAR values are first produced daily. Then, the 8 days composite corresponds to the values of the product when the maximum FAPAR value within the eight day period is observed. Note that FAPAR values are not retrieved over bare and very sparsely vegetated area, permanent ice and snow area, permanent wetland, urban and water bodies.

TERRA MODIS Collection 5 BELMANIP2.1 and DIRECT site extraction was performed at EMMAH, INRA. As the original MODIS projection is sinusoidal, the products were re-projected and extracted in Plate carrée (1/112° resolution) using the MODIS reprojection tool [63]. A value of FAPAR = 0 was set over pixels identified as bare and very sparsely vegetated area. Note that the initial TERRA MODIS quality flags were recombined to associate QA to the OLIVE extracts: QA = 0 indicates that the data were retrieved using the main algorithm without saturation, QA = 1 that they were retrieved using the main algorithm with saturation, QA = 2 using the back-up algorithm while a QA = 3 indicates that no product was retrieved.

3.3.1. Spatio-Temporal Domain Used for the Validation

The temporal domain used for the validation is defined by the starting and ending year for each product along with the QA threshold selected for considering a product as valid (Table 2). A pixel at a given date is considered valid if its QA value is below the QA threshold value given in Table 2.

The spatio-temporal sampling over the BELMANIP2.1 and DIRECT sites was adapted to each evaluation criteria to optimize statistics computation (Table 3): for those involving independent evaluation of the products (continuity, temporal consistency and stability, distributions), OLIVE uses the individual valid pixels extracted over a 49 × 49 pixels window centered over each BELMANIP2.1 site with their original temporal sampling (Table 2). This extent of 49 × 49 pixels was chosen to keep a good compromise between the size of the area that should be large enough to be the most representative, the homogeneity of the sites, and the corresponding data volume to be uploaded on the OLIVE platform. For the qualitative assessment of the temporal profiles of each site, the median value of the valid pixels over the central 3 × 3 pixels window is computed if at least 5 out of the 9 pixels are valid. If less than 5 pixels are valid, the data are declared invalid for the given date. The 3 × 3 pixel sampling is used to account for the geo-location uncertainties and the point spread function associated to the projected products making the spatial support more consistent across the investigated products [24].

For the pairwise comparison of products where both the spatial and temporal supports must be consistent, the median of the 3 × 3 valid pixels is first computed, as for the temporal profiles. Then, a monthly value centered on the 15th of the month is calculated using a weighted average of valid median values. The weighting function takes into account the temporal frequency and temporal resolution (R_t) of each product: the weight attributed to each day d corresponding to a valid product within a month is proportional to the number of days within the t ± R_t day period that belongs to the month.

Finally, for accuracy assessment based on the DIRECT sites, the median over the 3 × 3 valid pixels is first computed. Then, OLIVE provides two methods to make ground-based measurements and satellite product dates concomitant. First, the valid product that is the closest to the date of ground measurement acquisition (within a ±15 day period) is selected. This allows comparison with product values as similar as possible to the original ones. Second, the product is interpolated at the measurement date using the LOWESS algorithm [64] for a 2 month window span. This last process allows smoothing of the high frequency variability of the product while exploiting ground measurements that are not exactly located at the date of the product availability.

For both pairwise comparison of products and accuracy assessment, we compute the slope and intercept of the linear regression between the y and x axis. We also compute the associated root mean square error (RMSE) and the coefficient of determination (R²), and specify the number of data (N) that were used to calculate these quantities.

3.3.2. Continuity Assessment

Continuity is assessed as the fractions of valid (and non-valid) data as well as their distribution in space and time. It is computed over 49 × 49 valid pixels of the BELMANIP2.1 sites with the original temporal sampling (Table 2). The fraction of valid data for the considered period is provided at several levels including: global level, site level and continental level (results not reported here for the sake of brevity but available in the validation report). In addition to providing information on the continuity of the products, it also allows evaluating the degree of completion and reliability of the validation exercise.

The global distribution is detailed per QA level (Figure 2) and shows that GEOV1 and MODC5 provide product values in most situations, whereas MGVIo shows a large fraction of invalid data. Note that the back-up algorithm (QA = 2) for MODIS is used in about 5% of the cases, whereas the main algorithm is under saturation situations (QA = 1) in about 20% cases.

The spatio-temporal distribution of the fraction of invalid data is computed for each product at the monthly time step and per 10° latitude strips (Figure 3). As expected, it decreases around the equator because of large occurrence of clouds, as well as for the higher latitudes for GEOV1 and MGVIo due to short day length and poor illumination conditions in addition to the cloud occurrence. This is not observed for MODIS because of the exploitation of the replication of views for the higher latitudes that are fully exploited, in agreement with the results reported by [65].

The fraction of valid data is also evaluated per continent and site (Figure 4) and per biome type. The GLOBCOVER biome classification is used [42,43] with the classes being aggregated into larger biome types comprising of Evergreen Broadleaf Forest (EBF), Evergreen Needleleaf Forests (ENF), Deciduous Broadleaf Forest (DBF) and Non Forest (NF). GLOBCOVER was chosen in order to be consistent with the BELMANIP2 site selection and homogeneity verification. However, considering the broad categorization of the vegetation after aggregation, the impact of using a different land cover map should be marginal. The results on the distribution of the fraction of valid observations per biome type showed good consistency with the spatio-temporal distribution (Figure 3) and is detailed in the OLIVE validation report. It can be noted also that in recent years, smoothing and gap filling methods were developed to overcome the problem of missing data [16,66–68].

3.3.3. Temporal Profiles

The temporal profiles are produced for each BELMANIP2.1 and DIRECT site. They allow identifying qualitatively the main features of each product. Figure 5 presents a selection of 4 sites corresponding to the four main biomes considered in OLIVE (BELMANIP2.1, Sites number 5, 65, 165 and 320), the other temporal profiles being available in the full OLIVE report (WWW1). The temporal profiles show a general agreement of the seasonality for the non-forest and deciduous broad leaf forest sites (Figure 5). For the needleleaf forest and evergreen broadleaf forest, we observe some discrepancies, particularly for MGVIo, which shows more missing observations confirming the previous results on continuity assessment. MODIS and MGVIo show also high frequency variations associated to the shorter temporal resolution as compared to GEOV1. However, one of the main features observed is the difference in magnitude for MGVIo that shows almost systematically lower FAPAR values. This was already reported in previous studies [23,36–39] and could be explained mainly by differences in algorithm design where different assumptions on canopy structure and leaf optical properties are used.

3.3.4. Temporal Consistency: Smoothness and Stability

Considering the product temporal frequency, the temporal course of FAPAR and hence the products are expected to be smooth because the corresponding biophysical variable changes incrementally except in the case of disturbance. Smoothness was characterized by using the absolute difference δ between the product value at date t, P(t), and the mean value between the two bracketing dates δ = |1/2(P(t + Δt) + P(t)) − P(t)|, where Δt is the temporal sampling interval [24]. The difference δ is computed only if the two bracketing product values (P(t− Δt) and P(t + Δt)) exist. The smoother the temporal evolution, the smaller the δ difference should be. The smoothness δ is computed separately for each product for all the acquisition dates over the whole period considered (Table 3), and the 49 × 49 pixels centered on the BELMANIP2.1 sites using the original temporal sampling (Table 2). The results show that GEOV1 is generally smoother than MODIS and MGVIo, confirming the qualitative observations over the temporal profiles (Figure 6). This is mainly due to the longer temporal compositing window (30 days).

To better assess the temporal stability of the products, a global average is computed for each date over the 49 × 49 pixels of all the BELMANIP2.1 sites. The results (Figure 7) confirm the fair agreement between MODIS and GEOV1 FAPAR products with consistent year to year seasonality. Conversely, MGVIo shows both high and low frequency variability, for example, in June 2004, 2008 and 2009 MGVIo displays an anomaly in FAPAR seasonality that is not observed with MODC5 and GEOV1. A yearly climatology is finally computed to provide a global average seasonality (Figure 8). The difference in magnitude is the largest in July, corresponding to the maximum FAPAR value, with MODC5 > GEOV1 > MGVIo. Conversely, the spatial and inter-annual averages of the three products are in relatively good agreement in winter. Note however that the spatial representativeness may differ between the three products, depending on the fraction of valid data (Figures 2 and 3).

3.3.5. Distribution Per Biome and Continent

The histograms (probability density function) of each product were computed for the 4 main biome types (NF, DBF, ENF and EBF) over the valid pixels within the 49 × 49 pixel window centered over all the BELMANIP2.1 sites for the valid dates and for the considered period (Table 3) using the original temporal sampling (Table 2). Results show that MODC5 and GEOV1 have similar distributions while some differences are observed for Evergreen Needleleaf Forests (Figure 9). MODC5 shows an unexpected large number of almost bare soil situations. This is probably due to a misclassification of pixels into bare area, for which the MODIS algorithm is not triggered, and that were set to null FAPAR value in OLIVE. Conversely, GEOV1 shows a pronounced mode for high values of FAPAR. MGVIo modes are systematically shifted towards the low FAPAR values except for the Non Forest sites (Figure 9). Further, maximum FAPAR values for the forest biomes hardly reach FAPAR = 0.8 although such biomes, which have an almost complete cover fraction in many situations, should show higher values of FAPAR. In addition, the significant fraction of almost bare soil situations seen in GEOV1 and MODC5 for the Non Forest sites (Figure 9) is not seen by MGVIo.

The observed distribution over continents (only three are displayed here) and for the global distribution (Figure 10) show similar features as those observed over the main biome types: relatively good agreement between MODC5 and GEOV1, while MGVIo shows higher fraction of low FAPAR values including very low values.

3.3.6. Pairwise Comparison between the Three Products

The same spatial and temporal support is used to provide scatterplots between products: the median over the valid 3 × 3 pixels centered over each site is computed and then interpolated for fixed dates every 15th of the month during the considered period (Table 2). Scatterplots were generated globally for each pair of products (all BELMANIP2.1 sites, Figure 11) and as a function of the main biome types (results not presented).

Results show that MODC5 and GEOV1 are fairly consistent over the whole range of FAPAR values. However, an offset around 0.1 is observed for MODC5 for the very low GEOV1 FAPAR values as already noticed by [18]. The dissymmetry in the scattered points under the 1:1 line is mainly explained by the high frequency variability observed in MODC5, partly due to cloud contamination that negatively impacts the estimated FAPAR value (Figure 11). MGVIo shows consistent values with GEOV1 for the lowest FAPAR values while MODC5 is offset by 0.1 in the same conditions, in agreement with the offset value already observed with GEOV1. For the medium to high FAPAR values, MGVIo shows systematically lower values, compared to GEOV1 and MODC5, with the differences increasing with FAPAR values and reaching up to 0.3 for the highest FAPAR values.

3.3.7. Accuracy Assessment Using Actual Ground Measurements

Most of the DIRECT sites were sampled using Digital Hemispherical Photographs (DHPs) and processed using can-eye (WWW6) for 10:00 local solar time. DHP provides FIPAR (Fraction of Intercepted Photosynthetically Active Radiation) rather than FAPAR, but the difference between both is marginal, except for highly homogeneous canopies [11,45,46]. Therefore the consistency between satellite products and ground validation measurements is ensured at the exception of non-green elements. Indeed for the forest sites the impact is expected to be small since they have mainly been sampled during the peak of the growing season, the leaves at the top of the crown masking the trunks and branches. For the other sites, cases with large amount of non-green vegetation were discarded.

The comparison with ground-based measurements was performed using the two methods proposed in Table 3 that is either using the product value period the closest to the ground measurement date (±15 days) or the value interpolated at the date of the ground measurement. It should be noted that evergreen needle leaf or broadleaf forest sites were treated as special cases to accommodate situations when the ground campaign is outside the validation exercise period (Table 2). This allows exploiting additional ground measurements: for each year of the validation exercise time period, the median of the product values is computed at the date in the year at which the measurement acquisition was achieved. Indeed, these specific sites are expected to show some limited inter-annual stability, and can be checked against the temporal profiles available for all the DIRECT sites. They are identified as open symbols in the results (Figure 12).

Results show that the use of the closest observation reduces the number of available matches between ground measurements and available products (Table 4). This is mostly obvious for MGVIo where only 18 DIRECT sites could be used because of the larger fraction of missing data as compared to 60 when the products values were interpolated with time. For MODC5 and GEOV1, the performances measured with the “closest” product are very similar to those obtained with the “interpolated” product values (Table 4), except that the RMSE of MODC5 is much higher with the “closest” method. This is explained by the high frequency variability observed on MODC5 temporal profiles that is partly smoothed out using the interpolation method. Emphasis here is therefore on the performances as measured by the “interpolated” method since it allows more data points as well as smoothing out the high temporal variability that is not necessarily targeted when addressing accuracy assessment. The high temporal variability is one of the main components of the precision of the products. Note that the 3 × 3 spatial support already smoothed out the random variability associated to the products at the pixel level.

Detailed inspection of the results obtained with the method based on temporal interpolation shows that GEOV1 achieves the best performances with a RMSE value of 0.12 (Table 4) followed closely by MODC5 (RMSE = 0.14). MODC5 shows an intercept of 0.18 for the low FAPAR values (Figure 12), confirming the previous observations derived from the inter-comparison exercise (Figure 11). The comparison with ground measurements confirms MGVIo underestimation as already observed by [33] and [69], particularly for the high FAPAR values.

The RMSE values obtained in this study appear to be slightly above the targeted uncertainties values required within the Global Climate Observing System GCOS [1], i.e., the value of 0.05% or 10% for FAPAR. However the uncertainties associated to the ground measurements also need to be accounted for. Regarding the number of data points that are considered (Table 4), and given their distribution in the [0, 1] possible range over several vegetation types, the validation stage 2 is achieved, making significant progress toward stage 3 validation.