Evaluation of the Consistency of the Vegetation Clumping Index Retrieved from Updated MODIS BRDF Data

Abstract

The clumping index (CI) quantifies the foliage grouping within distinct canopies relative to randomly distributed canopies, which plays an important role in the vegetation radiative regime. Among the vegetation structure parameters, the global CI map can be retrieved by using multiangle remote sensing observations. The bidirectional reflectance distribution function (BRDF)/Albedo product (MCD43) of the Moderate-Resolution Imaging Spectroradiometer (MODIS) is the crucial input data of the global CI product, which provides validated spatiotemporal continuous directional reflectance data. To determine the impacts of updating the MCD43 product on the consistency and performance of global CI products, CIs retrieved from different MCD43 versions (Collection V005/V006, C5/6) were compared on a global scale and validated with field-measured CI data. The results showed that the global and seasonal comparisons of C5 and C6 CI data are generally consistent. Compared to C5 CI data, C6 CI data have improved quality with more main algorithm retrievals and fewer case of missing data. The comparisons over the field measurements indicate that both versions of CI data agree with field-measured CI data in terms of values and seasonal variations, while C6 CI data (R² = 0.89, RMSE = 0.05, bias = 0.02) are closer to field CIs than C5 CI data (R² = 0.80, RMSE = 0.07, bias = 0.03), indicating a higher accuracy for C6 CI data. The monthly CI is recommended for characterizing the overall seasonal patterns of surface CIs.

1. Introduction

Naturally, vegetation leaves are not randomly distributed; they are organized to different degrees according to vegetation structures, such as the distribution of vegetation landscapes and the configuration of crowns, branches, and leaves. The clumping index (CI) quantifies the relative degree of leaf grouping compared with a random distribution [1]. The clumping of leaves impacts the interaction between radiation and vegetation [2,3,4], which in turn affects vegetation growth and its carbon and water cycles [5,6,7]. For remote sensing applications, the CI reflects additional information on vegetation structure for mono-angle observations, which is the most extensive data acquisition mode to date. It not only influences the estimation of the leaf area index (LAI) by using mono-observed gap fraction information [8,9,10,11], but also affects the separation of sunlit and shadowed leaves [12], which play important roles in the interpretation and estimation of vegetation photosynthesis and evapotranspiration [6,13,14]. The CI is also very useful for ecological and meteorological models [15,16].

The multiangle capability of remote sensing instruments provides a means to obtain surface anisotropic reflectance that embeds structural information of vegetation surfaces [12,17,18]. To explore the relationship between the CI and reflectance anisotropy, three main steps must be taken. First, multiangle surface reflectance data are acquired. The anisotropic reflectance of surface scattering is usually described by the bidirectional reflectance distribution function (BRDF) [19], which is determined by surface characteristics and the illumination and viewing positions [20]. Multiangle satellite observations can be coupled with semiempirical models to retrieve the BRDF of a locality [21]. Second, the surface anisotropic reflectance is characterized by fewer feature parameters. Several angular vegetation indexes were proposed by using BRDF models to capture anisotropic surface features, e.g., the anisotropic flat index (AFX) calculated by normalizing the net scattering magnitude with isotropic scattering describes the dome-bowl BRDF patterns of the surface [22], and the normalized difference between hotspots and dark spots (NDHD) highlights the mutual shadow information in forward scattering [23]. Third, establishing the relationship between the angular index and the CI. As the hotspot and dark spot are sensitive to the mutual shadowing of vegetation, the NDHD has been successfully related to the CI under the circumstances of both model simulations and field measurements [12,23]. This relationship is essential to retrieve the CI from remote sensing observations and has been applied to obtain CI maps at global and regional scales [24,25,26,27,28,29,30]. Additionally, the AFX has also been related to the CI for the dome-bowl BRDF shape information contained in this index, which is used as a backup algorithm for CI-NDHD equations [26].

The Moderate-Resolution Imaging Spectroradiometer (MODIS) BRDF/Albedo product (MCD43) is produced at a 500 m spatial resolution by inverting multidate, multi-angular, cloud-free, atmospherically corrected, surface reflectance observations acquired by MODIS instruments onboard the Terra and Aqua satellites over a 16-day period [21,31]. The product is derived using a kernel-driven semiempirical BRDF model consisting of isotropic scattering, volumetric scattering, and geometric-optical (GO) scattering components [32,33,34]. When the multiangle observations are sufficient to invert model parameters, a high-quality full inversion can be conducted; otherwise, a low-quality magnitude inversion is performed by using available observations and a priori BRDFs based on recent high-quality MODIS retrievals at the same location [35]. The MCD43 product has been produced since the launch of Terra in 2000 and has been updated in several versions. The first MODIS MCD43 product was produced in 2000 at a 1-km spatial resolution every 16 days [21], and Collection V005 (C5) product was retrieved at a 500 m spatial resolution every 8 days. Collection V006 (C6) MCD43 product was released in 2016 and this product version was produced daily by weighting 16-day period observations as a function of data quality, angular geometry, and temporal distance from the day of interest [36]. The C6 product is expected to benefit from the improvement in upstream inputs, such as the various calibration improvements of MODIS Level-1B products [37,38]. The improvements in the MCD43 product from C5 to C6 mainly include the upgrades of the temporal processing strategy (retrieved every 8 days to daily), snow/snow-free status determination, quality assurance, and uncertainty identification and a priori BRDF database updating of the backup algorithm [35].

The C5/C6 MODIS BRDF/Albedo products has been widely used for modeling other parameters [39,40,41,42] and monitoring environmental changes [43,44,45,46]. As a major upstream input, the MODIS BRDF/Albedo products have also been used to map global and regional CI data [7,25,26,28,29,30,47,48,49], and these satellite-derived CI maps have been used to estimate a number of biophysical parameters, such as LAI [50,51,52], the fraction of absorbed photosynthetically active radiation by vegetation (FAPAR) [53], fractional vegetation cover (FVC) [54], gross primary production (GPP) [55,56] and solar induced fluorescence (SIF) [57]. The CI data have also been used in land surface models for ecological and meteorological applications [13,58,59]. The C6 MODIS BRDF/Albedo product presents a major improvement for environmental modeling and monitoring with higher temporal-resolution and data quality than that of the C5 product [35]. Wei, et al. [30] used the C6 MODIS BRDF/Albedo product to estimate global CI map from 2001 to 2017 and explored the improved capability of the CI time series in modeling and monitoring vegetation structural characteristics on a global scale. Obviously, the updating of MODIS BRDF/Albedo product will influence downstream products’ subsequent applications. Several comparisons have been conducted between C5 and C6 data for directional reflectance [60], snow BRDF Atlas [61], and fractional cover product [62]. However, the potential influences on MODIS CI retrievals from the update of C5 to C6 MODIS BRDF/Albedo products have not been comprehensively explored, despite the demand for widespread and reasonable application of MODIS CI data. In this paper, CIs retrieved from two versions of MCD43 (C5/C6) were compared and validated on a global scale by using MODIS data and field measurements. Section 2 briefly introduces the CI retrieval algorithms and dataset used in this study. Section 3 details the consistency and performance of different versions of CI retrievals. Section 4 discusses the uncertainty of the results. The conclusions are presented in Section 5.

2. Materials and Methods

2.1. MODIS BRDF/Albedo Product

The MODIS instruments are onboard the sun-synchronous Terra (launched in 1999) and Aqua (launched in 2002) satellites. Terra passes the equator in the morning, while Aqua passes it in the afternoon. The view swath of MODIS is 2330 km wide and covers the Earth every 1–2 days. The ground track repeat cycle is every 16 days. Thus, MODIS instruments can accumulate multiangle observations of the Earth. The obtained raw data are then radiometrically calibrated and atmospherically corrected to yield surface reflectance. High-quality multidate, multiangle surface reflectance data over a 16-day retrieval period are used to obtain the MCD43 product by fitting a kernel-driven semiempirical BRDF model to establish the surface anisotropic reflectance. The number of high-quality observations available, the root mean square error (RMSE) of the fit, and the weight of determination (WOD), which indicates whether the observations over the viewing and illumination geometry are well distributed to capture the surface anisotropy, are used to generate a high-quality full inversion retrieval. A magnitude inversion, using available observations with an a priori BRDF database, is applied if a high-quality full inversion BRDF cannot be retrieved due to poor quality input or insufficient sampling of the viewing hemisphere. While this magnitude of inversion often provides reliable results, it is always flagged as a poor quality retrieval.

The RossThick-LiSparseReciprocal (RTLSR) kernel-driven semiempirical BRDF model is used in the MCD43 operational processing system [21]. This type of model basically consists of three components, i.e., isotropic scattering, volumetric scattering, and GO scattering, and each component presents the product of a weight f_i (i.e., f_iso, f_vol and f_geo) and a kernel K_i (i.e., K_iso, K_vol and K_geo) [32,33,63]. The kernels are functions of the viewing and illumination geometry, which indicates physically distinct BRDF shapes. The weights are the spectrally dependent parameters inversed by multi-angular observations.

$$R\left({\theta }_v,{\theta }_s,\Delta \phi ,\lambda \right)=f_{iso}\left(\lambda \right)+f_{vol}\left(\lambda \right)K_{vol}\left({\theta }_v,{\theta }_s,\Delta \phi \right)+f_{geo}\left(\lambda \right)K_{geo}\left({\theta }_v,{\theta }_s,\Delta \phi \right)$$

(1)

where R(θ_v,θ_s,Δφ,λ) refers to the bidirectional reflectance at the viewing and illumination geometry of view zenith θ_v, illumination zenith θ_s, and relative azimuth Δφ in the waveband λ. The RossThick kernel derived from a single scattering approximation to radiative transfer theory for the assumption of dense leaf canopies describes the volume scattering contribution [33,64]. The LiSparseReciprocal kernel derived from a geometric-optical mutual shadowing model for the assumption of a sparse distribution of surface objects describes the GO scattering component [34,65].

The MODIS C5 and C6 500 m BRDF model parameter product (MCD43A1) and its corresponding quality product (MCD43A2) were used to estimate the MODIS CI in this study. The C5 product is produced every eight days using high-quality surface reflectance observations acquired by both the Terra and Aqua satellites over a 16-day period (up to four observations per day from each sensor), while the C6 product is retrieved daily by weighting all clear-sky surface reflectance observations acquired over a 16-day period [21,31,66]. The day of interest for the C6 product is the middle day of the period, i.e., the 9th day, while the day of interest for the C5 product is the first day of the period.

The three model parameters (isotropic, volumetric, and geometric optical) of the red band (620–670 nm) in MCD43A1 and the band-dependent BRDF quality of the red band in MCD43A2 are used to estimate the MODIS CI. The BRDF_Albedo_Band_Quality layer of the C5 MCD43A2 product and the BRDF_Albedo_Band_Mandatory_Quality_Band1 layer of the C6 product are used to identify the different quality flags for the multiangle reflectance data and to construct the quality assurance (QA) index of CI retrievals. The Snow_BRDF_Albedo layer of the MCD43A2 product is used to distinguish snow and snow-free pixels. The MODIS product data were downloaded from the NASA Earthdata Search Site (https://search.earthdata.nasa.gov/, accessed on 5 November 2021). As the C6.1 product was not completely released until this study was completed, only the C5 and C6 products were involved in the global analysis in this study. The C6.1 product is recommended for users as soon as the update is completed.

2.2. MODIS Land Cover Product

The MODIS land cover type product (MCD12Q1) provides a suite of global land cover data with eight different classification schemes at 500 m resolution and annual time steps. This product is created using supervised classification based on reflectance data [67,68,69]. The International Geosphere-Biosphere Programme (IGBP) classification data of the MCD12Q1 product were used in this study to determine the type of canopy crown shape, which is required by the CI algorithm. The non-vegetation land cover types, such as barrens, snow, and ice, were masked in this study.

2.3. Retrieval Algorithm of the MODIS CI

The CI retrieval algorithm using the hotspot-adjusted BRDF model [70] was exploited to calculate different versions of MODIS CIs. This algorithm framework includes the main algorithm and backup algorithm and constrains the retrieved CIs to the closed interval [0.33, 1.0] according to the theory of BRDF variability for vegetation canopy in the principal plane and, based on the relationship of the NDHD and CI. There are three main steps to retrieving the CI from MODIS multiangle reflectance data. First, the hotspot-adjusted BRDF model with two optimized hotspot parameters (C1 and C2) was used to fit the BRFs. Then, the reflectances in the direction of the hotspot (45° in backward scattering) and dark spot (45° in forward scattering) were reconstructed to calculate the NDHD. Third, the CI was retrieved in conjunction with IGBP data based on the CI-NDHD linear relationships established by the 4-scale model [24,71,72]. Each retrieved CI value was assigned a CI quality flag that is closely associated with the MODIS BRDF quality scheme. In particular, if the derived CI is outside the closed interval [0.33, 1.0] regardless of the high-quality flag from the MODIS BRDF product, a backup algorithm based on the AFX is adopted to process the so-called outliers.

As the hotspot effect of surface anisotropic reflectance tends to be underestimated by the MODIS operational RTLSR model, a method has been developed to improve the estimation of the hotspot reflectance [73,74,75]. This method exploits an exponential approximation of the hotspot kernel to correct the hotspot effect of the volumetric scattering component for the RTLSR model. Two prior hotspot parameter values (i.e., C1 = 0.7, C2 = 3.2) were used in the improved model (i.e., the hotspot-adjusted version of the RTLSR model (RTCLSR)) to reconstruct the reflectance of hotspots (ρ_HS) and dark spots (ρ_DS). Then, the NDHD can be obtained by the following:

$$\mathrm{NDHD}=\frac{{\rho }_{HS}-{\rho }_{DS}}{{\rho }_{HS}+{\rho }_{DS}}$$

(2)

According to the linear CI-NDHD equations in the red band [24,71], the CI can be obtained in association with a priori classification information. For coniferous forests,

$$\mathrm{CI}=-0.47\times \mathrm{NDHD}+0.80$$

(3)

For the other vegetation classes,

$$\mathrm{CI}=-1.23\times \mathrm{NDHD}+1.34$$

(4)

For the backup algorithm, the AFX was used to find the most similar BRDF archetype of the same land cover type and provide a most likely estimation of the CI for pixels outside the close interval [0.33, 1.0] in terms of the main algorithm. The AFX is calculated using the BRDF model parameters as follows:

$$\mathrm{AFX}=\frac{WSA\left(\lambda \right)}{f_{iso}\left(\lambda \right)}=1+\frac{f_{vol}\left(\lambda \right)}{f_{iso}\left(\lambda \right)}\times H_{vol}+\frac{f_{geo}\left(\lambda \right)}{f_{iso}\left(\lambda \right)}\times H_{geo}$$

(5)

where WSA(λ) is the white sky albedo and H_vol and H_geo are bi-hemispherical integral values of volumetric and GO scattering kernels, respectively. The corresponding CI value of the backup algorithm can be found through a CI-AFX look-up table (LUT) established by using the CI-NDHD relationship based on the average BRDF shape (i.e., BRDF archetype) determined by the AFX.

To analyze the C5 and C6 CI products on a monthly scale, temporal compositing is utilized to select the best retrievals and generate monthly CIs from 8-day (C5) and daily (C6) retrievals. Compositing mainly consists of two steps. The first step is selecting retrievals according to quality flags: the main algorithm retrievals with the best quality MCD43A1 data have the highest priority, followed by the main algorithm retrievals with good-quality MCD43A1 data, and if no main algorithm retrievals are available, backup algorithm retrievals are selected. Then, the monthly CI is calculated as the average of the selected retrievals, and corresponding quality flags are used to identify monthly QA data.

2.4. Field-Measured CIs

To validate the two versions of CI retrievals, we compiled a series of field-measured CI data (Figure 1) from various land cover types.

First, 48 representative field-measured CIs (Table 1) were collected to validate the monthly MODIS CIs retrieved from the two product versions. These field-measured CIs were refined from the database of Jiao et al. [26], who compiled a series of references [25,49,76]. These field CI data were measured by the tracing radiation and architecture of canopies (TRAC) instrument along transects of tens to hundreds of meters [25,49,77], which are usually considered to be the relatively appropriate data for comparisons with satellite-derived CIs [26]. A spatial representativeness evaluation using the variogram model based on high-resolution ETM+ images [78] was conducted for the compiled database, and 48 CI field measurements were ultimately selected (Table 1) and used to evaluate the MODIS CIs in this study. The field-measured CIs from these 48 sites are spatially representative of the MODIS CIs at a 500-m spatial resolution.

Then, 7 sites with seasonal profiles of field-measured CIs (Table 2) were collected to validate the seasonal variation in the two versions of CI data. Two main sources of these field-measured sites are He et al. [47] and Fang [79]. The Tonzi site is located in California, USA, with rainfall that is typically concentrated between October and May. This site is classified as oak savanna woodland, with scattered blue oak trees (Quercus douglasii) in the overstory and grasses and herbs in the understory [80]. The Radiation transfer Model Intercomparison (RAMI) pine site is located in Järvselja, Estonia, with poor growth conditions. The dominant tree species is Scots pine (Pinus sylvestris), with herbs and moss in the understory [81]. The TP39 and TP74 sites are located in Norfolk, Ontario, Canada, and are dominated by white pine (Pinus strobus L.), which were planted in 1939 and 1974. TP39 trees were thinned in 1983, which caused lower tree density and more understory vegetation species in this stand [82]. The Yatir site is located at the edge of the Negev and Judean deserts, Israel, with rainfall events usually occurring between November and March. The dominant tree species of this site is Aleppo pine (Pinus halepensis Mill.), and the forest consists of individual canopies for limited available water [83]. In the understory, sparse vegetation usually develops during the rainy season and disappears shortly thereafter [84]. The Honghe and Hailun sites are located in Heilongjiang Province, China [85]. The Honghe site is covered with paddy rice (Oryza sativa L.) at Honghe farm, and the rice is usually transplanted in late May, with flowering in early July, grain filling in early August, and maturation in early September [86]. The Hailun site is planted with maize (Zea mays L.), soybean (Glycine max L.), and sorghum (Sorghum bicolor). Maize is usually planted in early May, with tasseling in late July, grouting in late August, and maturation in late September. Sorghum and soybean have similar growth patterns to maize, while the maturity stage of soybean normally occurs in early September [87]. The TRAC was used to measure field CIs along transects [77] at the TP39, TP74, and Yatir sites, while digital hemispherical photographs (DHP) taken by fisheye lens were used to retrieve field CIs at the Honghe [86], Hailun [87,88], and RAMI pine [89,90] sites. The upward-pointing digital images [91,92] were used to estimate field CIs at the Tonzi site.

Table 1. Characteristics of 48 validation sites and corresponding MODIS CIs. IGBP is the classification code of the International Geosphere-Biosphere Programme (IGBP) classification scheme. Ω_e denotes the clumping of foliage elements, and γ_e is the needle-to-shoot area ratio. CI_Field and CI_MOD refer to field-measured CIs and MODIS-derived monthly CIs, respectively.

ID	Sites	Lat.	Lon.	IGBP	Ωe	γe	CIField	CIMOD		Date	Source
								C5	C6
1	SRF	49.250	−82.050	1	0.88	1.71	0.51	0.52	0.49	June 2001	Leblanc et al. [90]
2	Krasnoyarsk2	57.233	91.583	1	0.85	1.53	0.56	0.60	0.58	Summers 2000, 2001	Leblanc et al. [90]
3	Okarito	−43.200	170.300	5	0.87	1.4	0.62	0.71	0.66	January 2003	Walcroft et al. [93]
4	SETRES	34.902	−79.486	1	0.9	1.21	0.74	0.80	0.79	August 2003	Iiames et al. [94]
5	Hertford	36.383	−77.001	1	0.94	1.21	0.78	0.77	0.68	August 2003	Iiames et al. [95]
6	SOJP	53.916	−104.692	1	0.85	1.42	0.6	0.60	0.57	Summers 2003–2005	Chen et al. [96]
7	HJP75	53.875	−104.045	1	0.93	1.44	0.65	0.67	0.67	Summers 2003–2005	Chen et al. [96]
8	SOBS	53.987	−105.117	1	0.9	1.36	0.66	0.60	0.65	Summers 2003–2005	Chen et al. [96]
9	Mer Bleue	45.400	−75.500	3	0.87	1.36	0.627	0.72	0.63	August 2005	Sonnentag et al. [97]
10	Howland ME	45.210	−68.740	1	0.98	1.6	0.61	0.68	0.55	June 2007	Richardson (unpublished)
11	RAMI spruce	58.295	27.256	1	0.84	1.42	0.59	0.58	0.56	July 2008	Pisek (unpublished)
12	Tonzi	38.431	−120.966	4	0.82	1	0.82	0.80	0.79	September 2008	Ryu et al. [98]
13	QYZ	26.751	115.060	1	0.81	1.45	0.56	0.59	0.47	April 2009	Zhu et al. [49]
14	QYZ	26.749	115.059	1	0.79	1.45	0.54	0.59	0.47	April 2009
15	QYZ	26.746	115.066	1	0.74	1.45	0.51	0.59	0.52	April 2009
16	QYZ	26.742	115.062	1	0.77	1.45	0.53	0.59	0.54	April 2009
17	QYZ	26.742	115.058	1	0.76	1.45	0.53	0.59	0.54	April 2009
18	QYZ	26.740	115.059	4	0.87	1	0.87	0.74	0.99	April 2009
19	MES	45.323	127.543	3	0.97	1.5	0.65	0.54	0.63	July 2009	Zhu et al. [49]
20	MES	45.322	127.548	3	0.93	1.5	0.62	0.54	0.63	July 2009
21	MES	45.308	127.553	4	0.86	1	0.86	0.69 *	0.61 *	July 2009
22	MES	45.297	127.541	4	0.76	1	0.76	0.62	0.76	July 2009
23	MES	45.297	127.544	4	0.89	1	0.89	0.74	0.76	July 2009
24	MES	45.297	127.496	5	0.97	1.3	0.75	0.56 *	0.61	July 2009
25	MES	45.296	127.540	4	0.81	1	0.81	0.68	0.76	July 2009
26	MES	45.295	127.499	1	0.95	1.5	0.63	0.51	0.59	July 2009
27	MES	45.294	127.515	1	0.97	1.5	0.65	0.64	0.59	July 2009
28	MES	45.267	127.577	3	0.97	1.5	0.65	0.57	0.66	July 2009
29	MES	45.266	127.578	3	0.93	1.5	0.62	0.62	0.59	July 2009
30	MES	45.308	127.559	4	0.89	1	0.89	0.81	0.80	July 2009
31	TTS	29.855	121.740	4	0.91	1	0.91	0.82	0.69 *	September 2009	Zhu et al. [49]
32	TTS	29.854	121.738	1	0.93	1.4	0.66	0.60	0.69	September 2009
33	TTS	29.854	121.696	2	0.88	1	0.88	0.81	0.62 *	September 2009
34	TTS	29.854	121.701	1	0.9	1.5	0.6	0.60	0.53	September 2009
35	TTS	29.853	121.707	4	0.88	1	0.88	0.81	0.81	September 2009
36	TTS	29.843	121.748	4	0.88	1	0.88	0.82	0.90	September 2009
37	TTS	29.842	121.746	4	0.78	1	0.78	0.82	0.90	September 2009
38	TTS	29.804	121.798	7	0.83	1	0.83	0.75	0.95	September 2009
39	TTS	29.802	121.788	2	0.75	1	0.75	0.81	0.68	September 2009
40	TTS	29.796	121.804	7	0.91	1	0.91	0.82	0.54 *	September 2009
41	TTS	29.796	121.732	4	0.8	1	0.8	0.81	0.36 *	September 2009
42	TTS	29.784	121.806	1	0.94	1.4	0.67	0.60	0.65	September 2009
43	TTS	29.784	121.802	1	0.81	1.5	0.54	0.60	0.57	September 2009
44	TTS	29.783	121.810	3	0.89	1.5	0.59	0.60	0.55	September 2009
45	TTS	29.810	121.789	4	0.84	1	0.84	0.81	0.57 *	September 2009
46	TTS	29.807	121.787	2	0.85	1	0.85	0.80	0.80	September 2009
47	TTS	29.785	121.808	7	0.84	1	0.84	0.80	0.79	September 2009
48	TTS	29.778	121.762	7	0.9	1	0.9	0.81	0.71 *	September 2009

* MODIS CI retrievals are somewhat underestimated relative to field-measured CIs.

Table 1. Characteristics of 48 validation sites and corresponding MODIS CIs. IGBP is the classification code of the International Geosphere-Biosphere Programme (IGBP) classification scheme. Ω_e denotes the clumping of foliage elements, and γ_e is the needle-to-shoot area ratio. CI_Field and CI_MOD refer to field-measured CIs and MODIS-derived monthly CIs, respectively.

ID	Sites	Lat.	Lon.	IGBP	Ωe	γe	CIField	CIMOD		Date	Source
								C5	C6
1	SRF	49.250	−82.050	1	0.88	1.71	0.51	0.52	0.49	June 2001	Leblanc et al. [90]
2	Krasnoyarsk2	57.233	91.583	1	0.85	1.53	0.56	0.60	0.58	Summers 2000, 2001	Leblanc et al. [90]
3	Okarito	−43.200	170.300	5	0.87	1.4	0.62	0.71	0.66	January 2003	Walcroft et al. [93]
4	SETRES	34.902	−79.486	1	0.9	1.21	0.74	0.80	0.79	August 2003	Iiames et al. [94]
5	Hertford	36.383	−77.001	1	0.94	1.21	0.78	0.77	0.68	August 2003	Iiames et al. [95]
6	SOJP	53.916	−104.692	1	0.85	1.42	0.6	0.60	0.57	Summers 2003–2005	Chen et al. [96]
7	HJP75	53.875	−104.045	1	0.93	1.44	0.65	0.67	0.67	Summers 2003–2005	Chen et al. [96]
8	SOBS	53.987	−105.117	1	0.9	1.36	0.66	0.60	0.65	Summers 2003–2005	Chen et al. [96]
9	Mer Bleue	45.400	−75.500	3	0.87	1.36	0.627	0.72	0.63	August 2005	Sonnentag et al. [97]
10	Howland ME	45.210	−68.740	1	0.98	1.6	0.61	0.68	0.55	June 2007	Richardson (unpublished)
11	RAMI spruce	58.295	27.256	1	0.84	1.42	0.59	0.58	0.56	July 2008	Pisek (unpublished)
12	Tonzi	38.431	−120.966	4	0.82	1	0.82	0.80	0.79	September 2008	Ryu et al. [98]
13	QYZ	26.751	115.060	1	0.81	1.45	0.56	0.59	0.47	April 2009	Zhu et al. [49]
14	QYZ	26.749	115.059	1	0.79	1.45	0.54	0.59	0.47	April 2009
15	QYZ	26.746	115.066	1	0.74	1.45	0.51	0.59	0.52	April 2009
16	QYZ	26.742	115.062	1	0.77	1.45	0.53	0.59	0.54	April 2009
17	QYZ	26.742	115.058	1	0.76	1.45	0.53	0.59	0.54	April 2009
18	QYZ	26.740	115.059	4	0.87	1	0.87	0.74	0.99	April 2009
19	MES	45.323	127.543	3	0.97	1.5	0.65	0.54	0.63	July 2009	Zhu et al. [49]
20	MES	45.322	127.548	3	0.93	1.5	0.62	0.54	0.63	July 2009
21	MES	45.308	127.553	4	0.86	1	0.86	0.69 *	0.61 *	July 2009
22	MES	45.297	127.541	4	0.76	1	0.76	0.62	0.76	July 2009
23	MES	45.297	127.544	4	0.89	1	0.89	0.74	0.76	July 2009
24	MES	45.297	127.496	5	0.97	1.3	0.75	0.56 *	0.61	July 2009
25	MES	45.296	127.540	4	0.81	1	0.81	0.68	0.76	July 2009
26	MES	45.295	127.499	1	0.95	1.5	0.63	0.51	0.59	July 2009
27	MES	45.294	127.515	1	0.97	1.5	0.65	0.64	0.59	July 2009
28	MES	45.267	127.577	3	0.97	1.5	0.65	0.57	0.66	July 2009
29	MES	45.266	127.578	3	0.93	1.5	0.62	0.62	0.59	July 2009
30	MES	45.308	127.559	4	0.89	1	0.89	0.81	0.80	July 2009
31	TTS	29.855	121.740	4	0.91	1	0.91	0.82	0.69 *	September 2009	Zhu et al. [49]
32	TTS	29.854	121.738	1	0.93	1.4	0.66	0.60	0.69	September 2009
33	TTS	29.854	121.696	2	0.88	1	0.88	0.81	0.62 *	September 2009
34	TTS	29.854	121.701	1	0.9	1.5	0.6	0.60	0.53	September 2009
35	TTS	29.853	121.707	4	0.88	1	0.88	0.81	0.81	September 2009
36	TTS	29.843	121.748	4	0.88	1	0.88	0.82	0.90	September 2009
37	TTS	29.842	121.746	4	0.78	1	0.78	0.82	0.90	September 2009
38	TTS	29.804	121.798	7	0.83	1	0.83	0.75	0.95	September 2009
39	TTS	29.802	121.788	2	0.75	1	0.75	0.81	0.68	September 2009
40	TTS	29.796	121.804	7	0.91	1	0.91	0.82	0.54 *	September 2009
41	TTS	29.796	121.732	4	0.8	1	0.8	0.81	0.36 *	September 2009
42	TTS	29.784	121.806	1	0.94	1.4	0.67	0.60	0.65	September 2009
43	TTS	29.784	121.802	1	0.81	1.5	0.54	0.60	0.57	September 2009
44	TTS	29.783	121.810	3	0.89	1.5	0.59	0.60	0.55	September 2009
45	TTS	29.810	121.789	4	0.84	1	0.84	0.81	0.57 *	September 2009
46	TTS	29.807	121.787	2	0.85	1	0.85	0.80	0.80	September 2009
47	TTS	29.785	121.808	7	0.84	1	0.84	0.80	0.79	September 2009
48	TTS	29.778	121.762	7	0.9	1	0.9	0.81	0.71 *	September 2009

* MODIS CI retrievals are somewhat underestimated relative to field-measured CIs.

Table 2. Characteristics of validation sites with seasonal profiles of field-measured CIs.

Site	Latitude	Longitude	IGBP	Species	Method	Measurement Dates	References
Tonzi	38.43	−120.96	8	Oak savanna woodland	DP 1	2009–2010	Baldocchi, et al. [80]
RAMI pine	58.311	27.297	1	Scots Pine	DHP 2	2011	Kuusk, et al. [81]
TP39	42.710	−80.357	1	White Pine	TRAC 3	2011, 2012	Peichl, et al. [82]
TP74	42.707	−80.348	1	White Pine	TRAC	2011, 2012	Peichl, et al. [82]
Yatir	31.35	35.03	1	Aleppo pine	TRAC	2005, 2012, 2013	Sprintsin, et al. [83]
Honghe	47.652	133.522	12	Rice	DHP	2012, 2013	Fang, et al. [86]
Hailun	47.415	126.818	12	Maize, soybean, sorghum	DHP	2016	Fang, et al. [87]

¹ DP: digital photography. ² DHP: digital hemispherical photography. ³ TRAC: tracing radiation and architecture of canopies instrument.

Table 2. Characteristics of validation sites with seasonal profiles of field-measured CIs.

Site	Latitude	Longitude	IGBP	Species	Method	Measurement Dates	References
Tonzi	38.43	−120.96	8	Oak savanna woodland	DP 1	2009–2010	Baldocchi, et al. [80]
RAMI pine	58.311	27.297	1	Scots Pine	DHP 2	2011	Kuusk, et al. [81]
TP39	42.710	−80.357	1	White Pine	TRAC 3	2011, 2012	Peichl, et al. [82]
TP74	42.707	−80.348	1	White Pine	TRAC	2011, 2012	Peichl, et al. [82]
Yatir	31.35	35.03	1	Aleppo pine	TRAC	2005, 2012, 2013	Sprintsin, et al. [83]
Honghe	47.652	133.522	12	Rice	DHP	2012, 2013	Fang, et al. [86]
Hailun	47.415	126.818	12	Maize, soybean, sorghum	DHP	2016	Fang, et al. [87]

¹ DP: digital photography. ² DHP: digital hemispherical photography. ³ TRAC: tracing radiation and architecture of canopies instrument.

2.5. Experimental Design

The experimental design is illustrated through the use of a general flow chart (Figure 2), which organizes the main components in this study. This design includes four main parts: (1) retrieval of global CI data from different product versions, (2) spatial distribution comparison based on global monthly CI data, (3) seasonal variation comparison between global C5 8-day CI data and C6 daily CI data, and (4) validation and comparison based on globally distributed field-measured CI data.

For data processing, the C5 CI data were retrieved in the local server by downloading the MODIS data to a local database, while C6 CI data were directly retrieved on the Google Earth Engine cloud-based platform without downloading large amounts of MODIS data input to the local machine. After the retrieved CI data were obtained, we compared the spatial distribution and seasonal variation of these two versions of CI data on a global scale. As a major improvement in the C6 BRDF product is the daily inversion rather than an 8-day inversion of the C5 product, the different temporal resolutions of these two versions of CI data are compared for temporal variation. Moreover, these two versions of CI data are synthesized into monthly data for a global distribution comparison in terms of CI value and retrieval quality. Finally, the MODIS CI data retrieved from different product versions are validated and compared based on the collected field-measured CI data. The monthly CI data are validated by comparison with 48 spatially representative field-measured CI data, and the seasonal variabilities in the 8-day and daily CI data are evaluated using 7 field data sites with seasonal profiles of field-measured CIs. These different versions of MODIS CI retrievals are geolocated at the surface in terms of latitude and longitude.

3. Results

3.1. Comparison of C5 and C6 CI Products

3.1.1. Spatial Evaluation

The worldwide spatial distribution of CIs in January and July 2015 retrieved from C5 and C6 MCD43 data is shown in Figure 3. Figure 3a–d present global maps of CIs, from green to yellow, indicating CI values from low (clumped distribution) to high (random distribution), and the white, blue, and gray colors denote excluded non-vegetation area, water, and missing value pixels, respectively. A comparison between the global maps of C5 and C6 CIs (Figure 3a–d) shows that these two data exhibit similar global spatial patterns that are somewhat related to the distribution of vegetation types [67,69]. Tropical forest regions and middle-high latitude forests of the Northern Hemisphere have lower CI values, while the grass, shrub, crop, and savanna regions display higher CI values [24,26,99]. In the Northern Hemisphere, the January CI is generally higher than that in July, especially in middle-high latitudes, where the CI responds significantly to mature and senescent vegetation stages. In the Southern Hemisphere, the CI in January is lower than that in July, particularly in southern Africa, where the predominant vegetation is shrubs. These results are confirmed against the latitudinal distributions of these two versions of CI data.

As shown in Figure 3e,f, the latitudinal distributions of the C5 and C6 CIs in January and July exhibit similar profiles. In the higher latitudes, both the C5 and C6 CIs present clear seasonal differences; in line with the vegetation phenology [45,46], the CI in the mature stage is lower than the CI in the senescent stage [30,47]. In the tropics, there is no obvious seasonal difference in these latitudinal belts for the year-around ample illumination, except that the C5 CI is slightly higher in January among 10°N–10°S, and the C6 CI is slightly higher in July within the 10°S latitude and the Tropic of Capricorn. Comparing the two data versions reveals that the C5 and C6 CI profiles in July match relatively well (R = 0.86 for July and R = 0.83 for January), while the C5 CI overestimates the C6 CI in January within the latitude range of 10°N–30°S. The higher C5 CI in January may result from the relatively large number of pixels retrieved when using the backup algorithm compared to the C6 data.

3.1.2. Spatial Coverage of the Main Algorithm

The quality flag data stored in the CI products indicate the retrieval status of each pixel. Comparing the quality flag data on the global scale can help us further understand the overall quality of the C5 and C6 CI products. Figure 4a–d display the global maps of CI QA data corresponding to the CI maps in Figure 3a–d. The quality flags of 0–1 indicate the main algorithm retrievals by using the best- and good-quality BRDF model parameters, respectively; flag 2 represents the backup algorithm retrievals obtained from CI-AFX LUT, while the main algorithm retrievals with good- and high-quality MCD43A1 data are out of the physical range of [0.33, 1.0]; flag 3 denotes pixels labeled in snow status according to the MCD43A2 data, and the backup algorithm retrievals are used to reduce the influence of cover-snow on the anisotropic reflectance [100]. Pixels of excluded non-vegetation, water, and missing values are the same as those on the CI maps, and the quality flag of the missing value refers to a lack of BRDF model parameters, which indicates a limited number of clear-sky observations for retrievals [21,35]. Figure 3e,f show the latitudinal distribution of the retrieval rates of different algorithm paths, including the main algorithm, backup algorithm and non-retrievals (missing value). In each bar chart, the left panel is for the C5 CI, and the right panel is for the C6 CI.

C5 and C6 show similar patterns for the latitudinal distribution of global CI QA, while more high-quality pixels are acquired in C6 on a global scale, indicating an improvement in the quality of C6 retrievals. The overall proportion of main algorithm retrievals is higher in July than January due to the large number of main retrievals in the Northern Hemisphere summer. The proportion of main algorithm retrievals is higher in the Southern Hemisphere in January and higher in the Northern Hemisphere in July. Pixels in mid and low latitudes (−30°–30°N) consist of relatively more missing values due to cloud and aerosol contamination. Pixels at high latitudes (60°–90°N) contain even more missing values due to their long-term snow or cloud cover and large solar zenith angle.

Compared to C5, C6 shows a higher rate of main algorithm retrieval and a lower rate of missing value in the summer hemisphere (Southern Hemisphere in January and Northern Hemisphere in July) overall, which benefits from the improved algorithm of the C6 MOD43 data. In the Southern Hemisphere in January, the rates of the main algorithm retrievals and missing value are 89% and 9% for C6 CI data, and 86% and 12% for C5 CI data, respectively. In the Northern Hemisphere in July, the rates of the main algorithm retrievals and missing value are 89% and 8% for C6 CI data, and 89% and 10% for C5 CI data, respectively. The C6 MODIS BRDF/Albedo products use all available clear-sky observations, rather than the four observations per day used by C5 products for the limitation of storage processing volume; thus, the former acquires more high-quality retrievals. On the other hand, the improvement in the cloud detection algorithm in MODIS C6 products, combined with the updated threshold of the quality of directional observations (measured by the RMSE and the WOD), partly leads to a lower rate of the main algorithm retrieval of C6 at low latitudes (−10°–10°N) in July. Moreover, the C6 CI shows a higher proportion of missing values than the C5 CI in high latitudes (50°–70°) in the winter hemisphere (Northern Hemisphere in January and Southern Hemisphere in July). In addition, 42% of the C6 CI data fail to be retrieved for the winter hemisphere in both January and July, while the corresponding proportions of C5 CI data are approximately 25% and 33% in the winter hemisphere in January and July, respectively. This difference may be attributed to the improved threshold of snow, cloud, and cloud shadow detection in the C6 BRDF/Albedo product.

The rate of backup algorithm retrieval of the C6 CI is lower than that of the C5 CI, except for pixels south of the Sahara in July. Most CIs retrieved by the main algorithm are larger than 1 in south of the Sahara from the C6 CI data in July, and the backup algorithm is adopted to generate the corresponding reasonable CI values for these pixels according to an averaged BRDF shape using the AFX. The retrieved CIs from the backup algorithm of the C6 data are generally consistent with the corresponding C5 CIs. The C5 BRDF data of these backup algorithm CIs consist of more full inversions and missing values than the C6 data, which is considered to be the main reason for the different performances of these two versions of CIs. The C6 BRDF/Albedo products use all available clear-sky observations for retrievals and result in fewer missing values. The updated thresholds for high-quality observations may lead to fewer full inversions at low latitudes with high aerosol or cloud contamination [35]. Previous studies reported that the accuracy of CI retrievals from the main algorithm is higher than that of the backup algorithm for snow-free surfaces [26,100]. Considering the higher rate of the main algorithm retrieval of the C6 CI mentioned above, the overall inversion quality of the C6 CI is higher than that of the C5 CI.

3.1.3. Temporal Evaluation

Figure 5 shows the seasonal variations in MODIS C5 and C6 CIs averaged over the Northern Hemisphere (NH, 30–90°N), Southern Hemisphere (SH, 30–90°S), and tropics (Trop, 30°S–30°N) for different IGBP land cover types. In general, the C5 and C6 CIs show similar seasonal variations in these three latitudinal zones, which is consistent with the regional vegetation phenology [45,46]. The CIs show opposite seasonal variations over the NH and SH, for both the C5 and C6 CIs, while the variations in the CIs in the Trop are relatively stable. In the NH, the CIs show relatively large values (less clumped) in January, reach a minimum (most clumped) in July, and exhibit relatively large values at the end of the year. In contrast, in the SH, the CIs show relatively small values in January, reach a maximum in June, and fall back to small values in December, except for the croplands that show a slight increase in October. These seasonal variations in the CIs are highly related to regional vegetation phenology, in which low CI values generally arise in the peak growing season with large LAIs and a closed canopy, while high CI values usually appear in the dormant period with lower LAIs and a sparse canopy [47,100].

The C5 and C6 CIs also show similar values and seasonal variations for different vegetation types, which are associated with the morphology and canopy properties, as well as the growing patterns of different types. For both the C5 and C6 CIs over these three latitude regions, the CIs of evergreen/deciduous needleleaf forests (ENF/DNF) are smallest (mostly clumped); the CIs of open shrublands (OSh) are largest (least clumped); the CIs of closed shrublands (CSh), grasslands (GL), croplands (CL), and cropland/natural vegetation mosaics (CVM) are relatively high; and the CIs of evergreen/deciduous broadleaf forests (EBF/DBF), woody savanna (WSa), and savanna (Sav) are relatively low. The shoots of coniferous trees enhance the clumping effect of needle leaves compared with other vegetation types with flat, broad leaves; thus, coniferous forests have lower CI values than other land cover types. Land cover types containing trees with high canopies (>2 m) show smaller CIs than shrubs, crops, and herbaceous vegetation. For seasonal variations, generally, broadleaf forests are more variable than needleleaf forests, and the variations in mixed forests (MF) are in the middle for both versions of the CIs. Most of the coniferous forests (~96%) that were distributed in the NH show relatively smooth seasonal variations in the low needle turnover rate [47,101]. EBFs are mainly distributed in the Trop (~95%), while DBFs are primarily in the NH (~71%), and broadleaf forests in mid and high latitudes are more variable than those in the Trop due to the more stable climate conditions in the tropics. The seasonal variations in shrubs and savannas are commonly smaller than those in broadleaf forests but larger than those in herbaceous vegetation (i.e., CL and GL).

Compared with 8-day and daily CIs, the monthly CIs could better characterize the overall seasonal patterns of surface CIs, especially for C6 daily data, which present more fluctuations due to a higher temporal resolution and uncertainty, most likely using the daily BRDF data. Through quality screening and averaging, monthly CIs smooth the fluctuations caused by poor quality data and background reflectance changes (e.g., ephemeral rainfall) [47] to provide relatively reasonable CIs of surface vegetation for ecological and meteorological modeling at global and regional scales. Compared with monthly data, the 8-day and daily CIs in the NH fluctuate more at the beginning and end of the year, and the SH CIs fluctuate more in the middle of the year. These small waves are largely related to the decline in the main algorithm retrievals in winter for increased snow/cloud cover. Temporal variations in the CIs in the Trop are relatively stable. The 8-day and daily CIs of shrublands and the CI in the Trop exhibit relatively larger fluctuations than the monthly CIs, especially during May to August, when the main algorithm retrievals decline for increased precipitation. The daily CIs of ENF and MF, mainly found in southeastern North America and southeastern Asia, present short-term increases during August and September. These short-term variations may be attributed to tropical storms at that time, which may affect the main algorithm retrievals by destroying tree canopies and influencing the background reflectance [102,103]. Although the C6 daily CI shows larger fluctuations than the C5 8-day CI, it provides more information to capture rapid surface changes, and the corresponding monthly CI exhibits better continuity than the C5 data. The temporal patterns of the C6 monthly CIs of the EBF in the NH, the MF/ENF in the Trop and the MF/DBF in the SH are more stable and reasonable than the C5 monthly CIs.

3.2. Comparison with Field Measurements

3.2.1. Direct Validation

Figure 6 compares the MODIS monthly C5 and C6 CIs with field measurements. All the field CI data (Table 1) are evaluated to be spatially representative for 500 m MODIS pixels by using a method developed by Román et al. (2009) [78] and Wang et al. (2012, 2014, 2017) [104,105,106], which rationalizes a direct comparison between field-measured CIs and the corresponding MODIS CIs [26].

In general, both versions of the MODIS CIs have a high correlation (R² ≥ 0.80) with the field CI data, while the C6 CI (R² = 0.89, RMSE = 0.051, bias = 0.02, MAE = 0.03) better agrees with the field measurements than the C5 CI (R² = 0.80, RMSE = 0.065, bias = 0.03, MAE = 0.05). Both CI versions slightly underestimate the surface CIs (bias = 0.02 and 0.03 for C6 and C5 CI), while the C6 CI shows an improvement over the C5 CI in the higher range of values (0.75–0.95). Much of the improvement in the C6 CI may be due to the improved data quality of the MODIS C6 BRDF/Albedo product. Compared with the MODIS C5 BRDF/Albedo product, the C6 implementation enhances the estimation of surface reflectance anisotropy and improves the inversion quality [35]. On the one hand, all valid clear-sky observations are used for each retrieval in the C6 BRDF/Albedo inversion algorithm, rather than up to four observations used for each retrieval with the C5 product [66]; thus, the C6 product provides more available observations for BRDF inversion. On the other hand, in the case of insufficient observations for a full inversion, the strategy of pixel-specific a priori knowledge updated daily is used to improve the magnitude inversion of the C6 BRDF/Albedo product, instead of the static seasonal land cover BRDFs strategy used in the C5 inversion algorithm [107].

3.2.2. Evaluation of Seasonal Variability

Figure 7 shows the seasonal variabilities in the C5 8-day CI (red), C6 daily CI (blue), and field CI (black) data for the 7 field sites with seasonal profiles of field-measured CIs (Table 2). All valid MODIS CI values from the field sites during the field measurement year are shown, including the main algorithm (circles and squares denote retrievals of QA = 0 and 1, respectively) and backup algorithm (triangles) retrievals. Corresponding scatter plot for comparison of C5 8-day CI and C6 daily CI data with the field CI data can be found in Supplementary Materials (Figure S1).

In general, seasonal variations in the C5 and C6 MODIS CI data agree with the field-measured CIs, which reflect the states of vegetation growth and its intrinsic structures to some extent. The C6 CI data are closer overall to the field-measured CIs than the C5 CI data, while the C6 CI data present more temporal details and fluctuations due to their higher temporal resolutions. Both versions of the MODIS CI data and field-measured CI data show lower CI values (0.4–0.8) and small seasonal variations at evergreen conifer forest sites (i.e., TP39, TP74, Yatir, and RAMI pine) and higher CI values (0.5–1.0) and strong seasonal variations at cropland sites (i.e., Honghe and Hailun) and woody savanna sites (i.e., Tonzi), which are consistent with previous studies [30,47,89,100]. The relatively small seasonal change in the CIs for coniferous forests may be attributed to the low needle turnover rate [47,101]. For the TP39, TP74, and RAMI pine sites, the field-measured CIs are stable at approximately 0.55, while the C5 and C6 CIs fluctuate at approximately 0.55, and the C6 CIs are closer to the field-measured CIs than the C5 CIs. For the two cropland sites, field data were mainly collected in the growing season (days of year (DOYs) 170–240), and both the C5 and C6 CIs showed seasonal variations similar to the field-measured CI data, which exhibited lower CI values during the peak growing period and higher CI values at the beginning and end of the growing season [30]. At the Honghe site, paddy rice is usually transplanted in May (DOYs 120–150) with higher CI values, and subsequently, the CI values decrease with vegetation growth and return to higher values after the rice harvest in late August (DOYs 230–240). At the Hailun site, crops are generally sown in May with relatively high CI values, and then, the CI values decrease in June (DOYs 150–180) with vegetation growth and increase in July and August (DOYs 200–230) during the tasseling and milking stages; finally, they slightly increase in the mature stage in September (DOYs 240–270). In the fallow period, except for the freeze period (late October to April), the two versions of the MODIS CIs show relatively strong fluctuations before (DOYs 70–100) and after (DOYs 260–300) the growing season, which may be associated with the uncertainties in the corresponding MODIS BRDF data. The quality of MODIS BRDF data may be influenced by the increased heterogeneity of surface reflectivity before and after the growing season [108] and the poor angular sampling away from the principal plane in winter [109].

While the MODIS CIs generally follow the same trend as the field-measured CIs, there are some deviations in values, especially for the Yatir and Tonzi sites which have relatively sparse canopy distributions. The MODIS CIs mainly respond to the top structure of vegetation, while the field-measured CIs are usually influenced by understory vegetation [25,99]. At the Yatir sites, the field-measured CIs are higher (~0.7) in the middle of the year and lower (~0.5) at the beginning and end of the year, while the C5 and C6 CIs undulate at approximately 0.6. This site is located between arid and semiarid climatic zones, and the forest usually consists of individual canopies for the limited available water [83]. The relatively low field-measured CI values at the beginning and end of the year may be associated with the vegetation developed during the rainy season. At the Tonzi site, both versions of the MODIS CI data and field-measured CIs are lower from April to June and higher at the beginning and end of the year. This site generally consists of two layers of vegetation, a tree overstory and grass understory, that operate in and out of phase with each other over the course of a year [80,110]. The CI values decrease during DOYs 1–120, keeping pace with the growth of tree leaves, and the values increase slowly from the dry season until November (beginning of the rainy season). While the MODIS-derived CIs capture the overall trend of the seasonal variation caused by vegetation phenology, the satellite CI retrievals somewhat overestimate field-measured CIs [47,89]. The inconsistency between field-measured CIs and MODIS-derived CIs may be caused by the large difference in spatial scale between field measurements and satellite observations, as well as the influence of the change in background reflectance.

4. Discussion

4.1. Uncertainty in the MODIS CI Retrievals

The quality of MODIS-derived CI data relies heavily on the quality of input data (i.e., the MCD43A1/A2 and MCD12Q1 products) and the configuration of the retrieval algorithm. By using the improved MODIS C6 surface reflectance and cloud mask products and the enhanced daily retrieval algorithm, the C6 MCD43 products present improved data quality compared with the C5 products [35]. Correspondingly, the C6 CI data exhibit a higher proportion of main algorithm inversion and a lower proportion of missing values than the C5 CI data globally. The C6 CI data are also more consistent with the field-measured CI data than the C5 CI data are. Moreover, the input MCD12Q1 IGBP land cover data determine the crown shape in the retrieval algorithm and thus influence the accuracy of the CI data [25]. The land cover map provides the basic distribution pattern of the global CI data, as the needleleaf forest comprises shoot level clumping rather than the leaf element level clumping in broadleaf vegetation types. As the accuracy of the BRDF product has been implied by the CI QA data, caution should be exercised when determining the uncertainty of land cover data in the utilization of the CI product.

In addition to the uncertainties of the input data, the configuration of the retrieval algorithm is also crucial to the accuracy of the MODIS CI data. As the RTLSR BRDF model used in the operational MODIS BRDF inversion algorithm reportedly underestimates the reflectance at the hotspot [25,89], the RTCLSR model [74] is used in this study. The RTCLSR BRDF model rebuilds the hotspot magnitude and width with two parameters based on an exponential approximation of a physical hotspot kernel, which corrects the hotspot signatures with high accuracy and without distinct scale dependence [26]. The linear CI-NDHD equations have been widely used to retrieve global CI maps [25,28,29,30,111] and are considered to be the best guess for the surface CI [24,26]. As the CI-NDHD relationships are affected by sparse vegetation coverage and high background reflectance, uncertainties in the CI data in these circumstances should be treated with some caution [24,25,47]. Furthermore, terrain-induced shadows can also increase the uncertainty of CI retrievals, which deserves further investigation [24,49,111].

4.2. Uncertainty in Field-Measured CI Data

The field-measured CI data used in this study are collected from different sources, and uncertainties arising from different observational instruments and methods should be considered an error source in validating the MODIS CI data. The field measurement of CI data is usually based on observations of transmitted light of the canopy by using instruments, such as TRAC [77,112] and DHP [113,114]. Various methods have been developed to estimate the CI data from field measurements, such as the method based on gap size distribution [77] and its corrected method [112], the logarithmic average method [115], and the combined gap size and logarithm method [90]. While previous studies reported a high correlation among different methods, observations with different view zenith angles, the assumption of the leaf projection function and the segment size may influence the assessment of field CI values [76,116]. Moreover, the measurement conditions, such as sunlight, terrain, and observation height, may also influence the correct estimation of the field CI [76,117].

The different spatial scales of field-measured CI data and MODIS-derived CI data may introduce scaling errors when validating the MODIS CI data [25,78,89]. Generally, field-measured data are optimally integrated with high-resolution data to validate moderate-resolution products [118,119]. As there are no adequate high-resolution CI data to conduct this validation strategy, the homogeneity and representativeness of field sites are usually used to assess the field-measured CI data and reduce spatial scaling errors [25,26,29,49]. While these assessments of sites are limited by the use of high-resolution single-angular reflectance data, further studies with high-resolution multiangle data and retrieved CI data are urgently needed to perform comprehensive validations of moderate-resolution CI data.

4.3. Seasonal Variability in the CI

The results of this study indicate that the seasonal variability in the CI data has the potential to reflect the change in vegetation status, which confirms the results of previous studies [30,47,70,99,100]. While the variation in the CI is closely related to the LAI, which exhibits low CI values in the growing season and high values in the dormant period, the tree cover, canopy arrangement, understory vegetation, and background reflectance all contribute to the seasonal variation in the CI [7,30,47]. Compared with 8-day and daily CIs, the monthly CIs reduce the influence of poor-quality data and ephemeral background reflectance changes that can provide relatively reasonable surface CIs for ecological and meteorological modeling at global and regional scales.

Moreover, the spatial and temporal mismatch of field-measured CIs and retrieved MODIS CIs most likely lead to inconsistencies in the evaluation of seasonal variabilities. As seen in 7 sites with seasonal profiles of field-measured CIs, compared to field-measured CI data, satellite-derived CIs with larger spatial extents are more susceptible to changes in the surrounding and background reflectance in theory [47]. In addition, satellite data mainly respond to the top structure of vegetation, while the field-measured CIs are probably influenced by understory vegetation. Field-measured CIs influenced by understory vegetation may overestimate the forest CIs, as the CIs of shrubs are usually higher than those of trees [25,99]. These issues propose a huge challenge for the validation of satellite CI retrievals using field CI measurements.

The different temporal scales of the C5 8-day CIs, C6 daily Cis, and field-measured CIs should be considered sources of uncertainty in the comparison of seasonal variability. The temporal resolutions of the C5 8-day and C6 daily CI data are the same as those of the corresponding BRDF products. Both the C5 and C6 BRDF products are produced by using observations accumulated in a 16-day period, while the C5 algorithm retrieves data every 8 days with the first day being the day of interest, and the C6 algorithm retrieves data every day with the central day being the day of interest [35]. These temporal differences between the C5 CIs, C6 CIs, and field CIs may introduce temporal bias in the comparison of CI seasonal variations.

5. Conclusions

The CI is a crucial vegetation structure parameter that plays an important role in LAI estimations and the building of ecological, hydrological, and land surface models. Its seasonality is also critical in LAI mapping and ecosystem modeling. The MODIS BRDF/Albedo product is the most widely used input data for CI retrieval. As the newly released C6 product presents significant improvements in data quality and temporal resolution, it is imperative for the application of MODIS CI products to determine the impacts of updating the MCD43 product as an important input for the consistency and performance of CI retrievals. In this study, CI retrievals from the C5 and C6 MCD43 products were compared and validated in terms of spatial distributions and seasonal variations. The results are described as follows:

(1)

The C5 and C6 CI data show similar spatial distributions globally. These two versions of data exhibit similar spatial patterns and latitudinal distributions globally in January and July, which relate to the distribution of vegetation types and leaf-on/leaf-off seasons. Forest areas have lower CI values, while grass, shrub, crop, and savanna regions display higher CI values. In the Northern Hemisphere, the January CI is generally higher than that in July, and the opposite is true in the Southern Hemisphere.

(2)

The C5 and C6 CI QA data show similar patterns globally, while the C6 CI data have an improved quality with more main algorithm retrievals and fewer missing values. For both versions of data, the overall proportion of main algorithm retrievals is higher in July than in January, while the C6 data show a higher rate of main algorithm retrievals and a lower rate of missing value in the summer hemisphere and a lower rate of backup algorithm retrieval in most of the world.

(3)

In general, the C5 and C6 CI data show similar seasonal variations in the three latitude zones (NH, SH, and Trop) and different land cover types, which is consistent with the vegetation phenology. Through quality screening and averaging, the monthly CI data have higher data quality and are recommended to characterize the overall seasonal patterns of the surface CIs well, with less uncertainty than using C5 8-day, monthly, and C6 daily CI data.

(4)

Through a comparison with field-measured CI data, both versions of the MODIS CI data agree with the field-measured CIs and their seasonal variations, while the C6 CI data (R² = 0.89, RMSE = 0.05, bias = 0.02) show better consistency with the field measurements than the C5 CI data (R² = 0.80, RMSE = 0.07, bias = 0.03).

Remote sensing CI retrieval products are susceptible to the performance of upstream data. Comprehensive evaluation and analysis of updated derived data as input are crucial for choosing the correct and reasonable CI products for downstream applications for user communities. As improvements in remote sensing instruments and processing algorithms are ongoing, further studies are needed for updated products, such as the MODIS C6.1 product and even the Visible Infrared Imaging Radiometer Suite (VIIRS) BRDF/Albedo product [120]. In addition, the within-pixel variability of CI measurements in satellite CI product validation and the influence of spatial resolution on CI retrievals from different satellite sensors deserve to be further explored in near future.