A Phenologically Simplified Two-Stage Clumping Index Product Derived from the 8-Day Global MODIS-CI Product Suite

Abstract

The clumping index (CI) is a key structural parameter that quantifies the nonrandomness of the spatial distribution of vegetation canopy leaves. Investigating seasonal variations in the CI is crucial, especially for estimating the leaf area index (LAI) and studying global carbon and water cycles. However, accurate estimations of the seasonal CI have substantial challenges, e.g., from the need for accurate hot spot measurements, i.e., the typical feature of the bidirectional reflectance distribution function (BRDF) shape in the current CI algorithm framework. Therefore, deriving a phenologically simplified stable CI product from a high-frequency CI product (e.g., 8 days) to reduce the uncertainty of CI seasonality and simplify CI applications remains important. In this study, we applied the discrete Fourier transform and an improved dynamic threshold method to estimate the start of season (SOS) and end of season (EOS) from the CI time series and indicated that the CI exhibits significant seasonal variation characteristics that are generally consistent with the MODIS land surface phenology (LSP) product (MCD12Q2), although seasonal differences between them probably exist. Second, we divided the vegetation cycle into two phenological stages based on the MODIS LSP product, ignoring the differences mentioned above, i.e., the leaf-on season (LOS, from greenup to dormancy) and the leaf-off season (LFS, after dormancy and before greenup of the next vegetation cycle), and developed the phenologically simplified two-stage CI product for the years 2001–2020 using the MODIS 8-day CI product suite. Finally, we assessed the accuracy of this CI product (RMSE = 0.06, bias = 0.01) via 95 datasets from 14 field-measured sites globally. This study revealed that the CI exhibited an approximately inverse trend in terms of phenological variation compared with the NDVI. Globally, based on the phenologically simplified two-stage CI product, the CI_LOS is smaller than the CI_LFS across all land cover types. Compared with the LFS stage, the quality for this CI product is better in the LOS stage, where the QA is basically identified as 0 and 1, accounting for more than ~90% of the total quality flag, which is significantly higher than that in the LFS stage (~60%). This study provides relatively reliable CI datasets that capture the general trend of seasonal CI variations and simplify potential applications in modeling ecological, meteorological, and other surface processes at both global and regional scales. Therefore, this study provides both new perspectives and datasets for future research in relation to CI and other biophysical parameters, e.g., the LAI.

1. Introduction

The vegetation canopy represents the most direct and dynamic interface between vegetation and the external environment, significantly influencing ecosystem material and energy exchanges as well as biodiversity [1]. In most natural vegetation canopies, leaves are not randomly distributed but exhibit varying degrees of aggregation at different levels, such as in tree clusters, crowns, branches, and leaf clusters [2]. The clumping index (CI) is a crucial structural parameter that quantifies the degree to which the leaf distribution deviates from randomness within the vegetation canopy [3,4]. The CI is very important for the accurate inversion of other biophysical parameters, e.g., the leaf area index (LAI). According to Chen and Black, the CI is defined as the ratio of the effective leaf area index (LAIe) to the true leaf area index (LAI), represented by the formula CI = LAIe/LAI [5]. Thus, the CI can be utilized to improve the accuracy of leaf area index inversion [6]. When leaves are randomly distributed, the CI equals 1; when they are regularly distributed, the CI is greater than 1; and when they are clumped, the CI is less than 1. CI and LAI are fundamental for distinguishing between sunlit and shaded leaves within the canopy, making CI critical for improving the accuracy of terrestrial carbon cycle models [7]. For example, the CI can improve the accuracy of global gross primary productivity (GPP) estimates. When the LAI is used directly without considering clumping effects, the global GPP is overestimated by 12%, whereas using the LAIe without accounting for clumping effects results in a 9% underestimate of the global GPP [8]. The clumping effect also influences the interception and distribution of light and precipitation within the canopy, thereby affecting photosynthesis [9,10,11] and transpiration [12]. Additionally, the amount of radiative energy available for evapotranspiration in vegetation is directly influenced by the degree of canopy clumping [13]. Therefore, CI plays a crucial role in the study and modeling of ecological processes, including light propagation, surface hydrology, photosynthesis, and vegetation productivity, at both regional and global scales [14,15,16].

Numerous studies have demonstrated that multi-angle data provide substantial structural information about vegetation canopies [5,7,17,18,19]. Consequently, many researchers have estimated the CI via multi-angle satellite data, including those from POLDER [20,21,22], MODIS [23,24,25,26,27], and the multiangle imaging spectroradiometer (MISR) [28]. Chen et al. [20] developed a linear regression model linking the normalized difference hotspot and darkspot index (NDHD) to CI based on a four-scale geometrical optical model under various conditions. Employing this model, a global CI product with a 6 km resolution was first generated, covering eight months of POLDER-1 data [20]. Building on these results, the terrain effects on CI inversion were corrected via POLDER-3 data [21]. However, because the operational MODIS algorithm underestimates the hotspot effect [26], Zhu et al. [24] refined the algorithm and generated a 500 m CI product for China. Similarly, the normalized difference vegetation index (NDVI) and solar zenith angle were applied to correct the MODIS hotspot reflectance, producing a global CI product for 2006 with a 500 m resolution, but only a single CI median value for the entire year was provided and released, posing challenges for time-sensitive models [23]. Pisek [28] subsequently applied the hotspot correction method proposed by He et al. [23] to MISR data for CI inversion, but a global CI product was not provided because of the MISR multiangular acquisition mode and capability. Wei et al. [27,29] evaluated the impact of different BRDF model configurations and solar zenith angles on CI inversion and then produced a global CI product at a 500 m resolution for the years 2001–2017 based on MODIS BRDF products. Jiao et al. corrected the hotspot effect of the MODIS RossThick-LiSparseReciprocal (RTLSR) model with greater accuracy, which has the potential to improve the accuracy of CI inversion [30]. They developed a backup algorithm based on the AFX and BRDF archetypes to limit outliers in CI inversion results, producing global MODIS CI products with 500 m resolution for monthly and 8-day intervals from 2001 to 2020 [31]. These products provide valuable references and data support for applications in land surface and ecological models, as well as for studying the spatiotemporal characteristics of CI.

As seasons change, vegetation undergoes a series of phenological variations, including sprouting, leaf development, flowering, fruiting, and defoliation, which is not only accompanied by changes in vegetation greenness but also by changes in canopy structure. Currently, vegetation phenological remote sensing monitoring technology has become highly advanced after decades of development. Due to its ability to significantly eliminate noise caused by observation conditions such as solar zenith angle, atmospheric and soil background, NDVI and EVI (Enhanced Vegetation Index), which characterize vegetation greenness, have become the primary parameters used in vegetation phenological remote sensing monitoring. As a crucial structural parameter reflecting the degree of nonrandomness in the spatial distribution of canopy foliage, the CI theoretically exhibits distinct phenological variations. However, owing to the limited data available and the relatively recent initiation of quantitative research on CI, there has been some debate regarding its phenological characteristics. Some researchers argue that the CI does not exhibit seasonal variations and thus is considered constant throughout the year [32,33,34], which impacts LAI estimation and the accuracy of ecological modeling [35]. However, more researchers believe and have started studying the seasonal variations in CI. Previous studies have shown that the timing of maximum clumping generally occurs earlier than the timing of maximum greenness of vegetation leaves [36]. This is likely related to the gradual increase in chlorophyll concentration, which represents greenness, during vegetation growth, while canopy leaves tend to reach their clumping state earlier [37,38]. Long-term field measurements of CI indicate that driven by vegetation phenology, CI in vegetation types such as evergreen needleleaf forests [39,40], deciduous broadleaf forests [2,41], savannas [42], grasslands [41], and crops [43,44] exhibit pronounced seasonal variations. Additionally, some researchers have studied the seasonal variations in CI via multiangle remote sensing data. Studies indicate that the CI exhibits significant seasonal variations at both regional and global scales. For example, Zhu et al. [45] demonstrated that CI in China clearly exhibits seasonal variations. Hu et al. [46] reported that while the interannual variation pattern of the CI in the Sanjiang Plain is not pronounced, the seasonal variations are significant. Wei et al. [27] reported significant seasonal variations in CI on a global scale. Yin et al. [47] compared the temporal and spatial characteristics of various global CI products and reported that both types of CI data exhibit seasonal variations. In summary, despite numerous studies that have analyzed CI over time to capture its seasonal and interannual variations, most research has been limited to examining the fluctuation characteristics of the average CI at regional or global scales. In general, CI values with distinct seasonal variation characteristics that are useful for ecological modeling have yet to be developed. Therefore, there is an urgent need to develop a CI product that more precisely reflects seasonal variations, providing direct and effective support for modeling various ecological, meteorological, and other surface processes.

This work employs the discrete Fourier transform and an improved dynamic threshold method to demonstrate that the CI exhibits seasonal variations, which provides evidence that leads to the development of a phenologically simplified two-stage MODIS CI product that reflects the general trend of CI seasonal variations. Section 2 comprehensively describes the various datasets and methods employed in this study, including a detailed explanation of the algorithm and band settings of the phenologically simplified two-stage MODIS-CI product. Section 3 elaborates on the seasonal variations in the CI, presenting the spatiotemporal characteristics and validation results of the phenologically simplified two-stage MODIS CI product. Section 4 discusses the uncertainties associated with the results. Finally, Section 5 presents the main conclusions of this study.

2. Datasets and Preprocess

2.1. MODIS CI Product

The global MODIS CI product used in the study was adapted from Jiao et al. [31], with a spatial resolution of 500 m, and a temporal resolution of 8 days and monthly intervals from January 2001 to December 2020. The product was derived from the MODIS BRDF and its quality assurance (QA) products (MCD43A1 and MCD43A2) and the MODIS land cover products (MCD12Q1) via primary and backup algorithms. Areas without vegetation cover, such as ice and snow (IGBP15), barren or sparsely vegetated areas (IGBP16), and water bodies (IGBP17), were masked [31]. The first band contains CI values with a valid range of 3300–10,000 and a scale factor of 0.0001; thus, CI values range from 0.33 to 1.0. The second band contains QA data, with values ranging from 0 to 3 that indicate quality from high to low. This dataset can be accessed at http://www.geodata.cn (accessed on 1 November 2023).

2.2. MODIS Land Surface Phenology Product

The MODIS land surface phenology (LSP) product (MCD12Q2, V061) provided by NASA is an annual dataset with a spatial resolution of 500 m that estimates phenological parameters via spline fit based on the 2-band enhanced vegetation index calculated from MODIS nadir BRDF adjusted surface reflectances (NBAR-EVI2). Previous comparisons and validations have shown that MCD12Q2 exhibits high accuracy and good generalizability [48,49,50,51] and has been applied in many studies worldwide. In this study, the “greenup” and “dormancy” bands of MCD12Q2 were used to represent the SOS and EOS, respectively.

2.3. MODIS NDVI Product

The NDVI dataset from MOD13A1 comprises two bands, the NDVI and EVI, which are derived from atmospherically corrected bidirectional reflectance. It has a temporal resolution of 16 days, dividing the year into 23 periods, with a spatial resolution of 500 m, reflecting the spectral and phenological features of the vegetation [52]. In this study, the maximum value composite method [53] was used to generate the monthly NDVI.

2.4. Land Cover Products

The MODIS land cover type product (MCD12Q1) provides annual global land cover maps with a 500 m spatial resolution for 2001–present. It includes eight land cover classification schemes derived via a supervised decision-tree classification method [54]. The International Geosphere–Biosphere Programme (IGBP) classification from the MCD12Q1 product was used in this study for consistency with the MODIS CI product.

Furthermore, the first global land cover datasets at 30 m resolution (GlobeLand30), obtained from the official online platform at https://www.webmap.cn/ (accessed on 15 January 2024), were selected as auxiliary data for the selection of typical pixels in this study. Additionally, these datasets were employed to further reduce the uncertainty in the land cover types of the CIs via a classification scheme that includes 10 first-level types, such as croplands, forests, and grasslands [55]. The 2000, 2010, and 2020 versions have been published to date.

2.5. Field-Measured CIs

The field-measured CI data derived from ground-based observation for validation (GBOV) were selected to validate the accuracy of the phenologically simplified two-stage MODIS CI product in this study. GBOV data are openly accessible from the GBOV portal (https://gbov.acri.fr, accessed on 1 February 2024) [56]. To ensure the quality of field-measured CIs and the reliability of continuous measurements, we controlled the data quality and errors and used MCD12Q2 to remove field-measured CIs that could not cover more than 3 months across different phenological stages. Ultimately, we obtained a total of 95 datasets from 14 field-measured sites for accuracy assessment with the MODIS time-share two-stage CI product. The distribution of the field-measured CI data is shown in Figure 1, with specific information provided in Appendix A.

3. Methods

The aim of this study is to conduct an in-depth analysis of the seasonal variations in the CI and to develop a phenologically simplified two-stage MODIS CI product that can provide global long time series pixel-by-pixel CIs with typical values of seasonal variations based on the MODIS 8-day CI product, thereby offering an effective reference for the modeling of various types of ecological, meteorological, and other surface processes. The general flow chart is shown in Figure 2. The study includes three parts: (1) due to the uncertainties in the CI inversion process, this study uses the more stable NDVI as a reference for comparison, estimates the phenological parameters of the CI and NDVI time series in the mid-to-high latitudes of the Northern Hemisphere (23.5°N–60°N) to identify CI seasonal variations relative to the NDVI; (2) determining the algorithm for developing the MODIS time-share two-stage CI product and evaluating its accuracy; and (3) exploring the seasonal differences in CI values and their spatiotemporal characteristics on the basis of the MODIS time-share two-stage CI product.

3.1. Study on the Seasonal Variation Characteristics of CI

3.1.1. Typical Pixel Selection Rules

Based on the MCD12Q1 and GlobeLand30 datasets, 90 typical pixels were selected in the mid-to-high latitudes of the Northern Hemisphere (23.5°N–60°N), where vegetation exhibits significant seasonal variation, remote sensing data quality is high, and the region has strong ecosystem representativeness, to validate whether CI exhibits seasonal variation characteristics. Six typical pixels were selected for each IGBP class, including annual and biannual croplands differentiated based on MCD12Q2 (Figure 1), to further explore the ability of the CI to capture the growth status of annual and biannual vegetation. Hypothesis testing on the sample size of the selected typical pixels revealed that the 90-pixel sample size in the study area is sufficient to meet the requirements of this study at the 95% confidence level. This study defines the typical pixel selection rules as follows: (1) the IGBP class in the MCD12Q1 product remains consistent from 2001 to 2020; and (2) for 500 m pixels, the percentage of GlobeLand30 land cover types corresponding to the same IGBP class in the MCD12Q1 product should be 70% or greater. These criteria ensure the representativeness of the selected typical pixels and help avoid the problem of heterogeneity, with mixed pixels identified at a high spatial resolution, which can obscure seasonal variations. Additionally, to minimize noise and reduce the uncertainty of the CI inversion results due to the nearly 20-year time span, this study calculates the mean values of the CI and NDVI values in a 3 × 3 pixel window surrounding each typical pixel. The distribution of these typical pixels is shown in Figure 1.

3.1.2. Discrete Fourier Transform

Owing to the large amounts of noise in vegetation index time series caused by cloud cover, atmospheric aerosols, the solar irradiation angle, and the bidirectional reflectance of features [57,58], selecting a fitting method that suits the characteristics of the time series when estimating phenological parameters is crucial. The Fourier transform treats time series as a fit of a series of sine and cosine functions with different frequencies. Given that vegetation growth exhibits periodic variation over time, the Fourier transform has been widely used in the fields of phenological parameter estimation and vegetation index time series reconstruction [59,60,61,62]. The CI and NDVI values of typical pixels are 8-day discrete datasets suitable for the use of the discrete Fourier transform method. Owing to the limitations of the inversion process (e.g., lack of hotspot observations, especially in a cross principal plane (CPP)), the CI time series contains more potential noise. The discrete Fourier transform, as a finite time series algorithm, is theoretically significant and can more effectively mitigate noise. Therefore, the discrete Fourier transform (DFT) is chosen for the fast and effective fitting of the CI and NDVI time series. The DFT is given by Equation (1):

$$F_{(u)}={\displaystyle \frac{1}{N}}{{\displaystyle \sum _{x=0}^{N-1}{f}_x\times e}}^{-2\Pi ux/T},$$

(1)

which can be decomposed into a real part [Equation (2)] and an imaginary part [Equation (3)]:

$$F_{C(u)}={\displaystyle \frac{1}{N}}{\displaystyle \sum _{x=0}^{N-1}\left[f\left(x\right)\times \cos \left(2\pi \frac{ux}{T}\right)\right]},$$

(2)

$$F_{S(u)}={\displaystyle \frac{1}{N}}{\displaystyle \sum _{x=0}^{N-1}\left[f\left(x\right)\times \sin \left(2\pi \frac{ux}{T}\right)\right]},$$

(3)

with the magnitude calculated as:

$$F_{(\mathrm{Magnitude})u}=\sqrt{F_{C\left(u\right)}^2+F_{S\left(u\right)}^2},$$

(4)

and the phase is calculated as:

$$F_{\left(\mathrm{Phase}\right)u}=\tan \left({\displaystyle \frac{F_{C\left(u\right)}}{F_{S\left(u\right)}}}\right).$$

(5)

where f(x) is the xth sample value, u is the number of Fourier components or harmonics, x is an index representing the current sample number, and T is the length of the period covered.

3.1.3. Improved Dynamic Threshold Method

The dynamic threshold method adopts a dynamic ratio that is closely tied to the amplitude of seasonal variations in the time series, effectively reducing the uncertainty in phenological parameter estimations [63]. When the CI time series descends or ascends to a specific percentage of the annual CI amplitude, it is marked as the SOS or EOS, respectively, whereas for the NDVI, the process is the reverse, which indicates that a larger NDVI value tends to have a smaller CI value for a vegetation canopy. As a vegetation structural parameter, the CI is affected mainly by hotspot reflectance [23,24] and mixed pixel heterogeneity [20], leading to complex noise anomalies in the time series, even after fitting, compared with the NDVI. To address this, an improved dynamic threshold method was developed based on the characteristics of the CI time series. For noise with a CI fluctuation amplitude less than 0.1 of the ipsilateral amplitude, the SOS or EOS is defined as the mean value of multiple moments corresponding to a certain percentage of the current year’s amplitude during ascent or descent, further reducing noise. In addition, Jonsson et al. [64], who proposed the dynamic threshold method, suggested that the threshold value of the SOS is approximately 20%, but the selection of the threshold value varies according to the vegetation index, land cover type and other factors. Therefore, this study broadens the threshold range, from 10% to 90%, based on the characteristics of CI time series fluctuations, using the MCD12Q2 product for validation to determine the optimal threshold for estimating phenological parameters from the CI and NDVI time series.

3.2. Development of the MODIS Time-Share Two-Stage CI Product

3.2.1. Statistical Analysis of Phenological Stage Classification

As a crucial structural parameter reflecting the degree of nonrandomness in the spatial distribution of canopy foliage, the CI, in theory, should exhibit distinct phenological variations. Previous studies have confirmed that the MODIS LSP product (MCD12Q2, V061) has high accuracy and good generalizability in monitoring land surface phenology [48,49,50,51]. A diagram of the phenometrics retrieved for a single hypothetical vegetation cycle of the MODIS LSP product is shown in Figure 3. Therefore, in this study, a vegetation cycle was divided into four phenological stages based on the four phenological dates (greenup, midgreenup, midgreendown, and dormancy) provided by the MODIS LSP product: the leaf-off season (LFS, after dormancy and before greenup of the next vegetation cycle), the leaf-on season (LOS, from greenup to dormancy), the leaf-scattering season (LSS, from greenup to midgreenup and from midgreendown to dormancy), and the leaf-gathering season (LGS, from midgreenup to midgreendown), where the LOS includes the LSS and LGS. Li et al. estimated seasonal CI values via this classification method, which improved the high GPP estimates of the TL-LUE model to some extent [65]. In this study, we estimated the mean CI values of different IGBP classes for each phenological stage globally in 2020 and used one-way analysis of variance (ANOVA) and Tukey’s post hoc test to comparatively assess the variability of the mean CI values within each phenological stage, which in turn comprehensively considered the feasibility of dividing a vegetation cycle into two phenological stages (LFS and LOS) or three phenological stages (LFS, LSS and LGS). The results are expressed as the means ± standard deviations. In addition, “*” indicates significant differences in the mean CIs for the same IGBP classes and different phenological stages (* p < 0.05, ** p < 0.01, *** p < 0.001).

3.2.2. MODIS Time-Share Two-Stage CI Product

The processing flow and algorithm framework of this MODIS time-share two-stage CI product are shown in Figure 2 and

Table 1. Development algorithm and quality assurance (QA) system of the MODIS time-share two-stage CI product.

Most Frequent QA in 8-Day CI	Output QA	Output CI
		The First Algorithm	The Second Algorithm
0	0	$C{I}_t={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^{n}C{I}_i},\left(CI\_Q{A}_i=0\right)$	$C{I}_t={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^{n}C{I}_i},\left(CI\_Q{A}_i=0\right)$
1	1	$C{I}_t={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^{n}C{I}_i},\left(CI\_Q{A}_i=1\right)$	$C{I}_t={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^{n}C{I}_i},\left(CI\_Q{A}_i\le 1\right)$
2	2	$C{I}_t={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^{n}C{I}_i},\left(CI\_Q{A}_i=2\right)$	$C{I}_t={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^{n}C{I}_i},\left(CI\_Q{A}_i\le 2\right)$
3	3	$C{I}_t={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^{n}C{I}_i},\left(CI\_Q{A}_i=3\right)$	$C{I}_t={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^{n}C{I}_i},\left(CI\_Q{A}_i\le 3\right)$
32,765–32,767 (partial)	4	$C{I}_t={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^{n}C{I}_i},\left(CI\_Q{A}_i\le 3\right)$	$C{I}_t={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^{n}C{I}_i},\left(CI\_Q{A}_i\le 3\right)$
32,765–32,767 (all)	32,767	32,767	32,767

. The MODIS LSP and MODIS 8-day CI products are used as input data to retrieve the SOS and EOS of each year, based on the “greenup” and “dormancy” bands of the MODIS LSP product. All CI values and corresponding QAs of the MODIS 8-day CI product in LOS and LFS in turn and the QA with the highest occurrence frequency in the MODIS 8-day CI will be set as the QA of the MODIS time-share two-stage CI. If the QAs ranging from 32,765 to 36,767 (representing snow, barren or sparsely vegetated areas, and water bodies) occur most frequently and contain other QAs in the MODIS 8-day CI, then the QA of the MODIS time-share two-stage CI is set to 4. If all the QAs are 32,765–36,767 8-day MODIS CIs, then both the MODIS time-share two-stage CI and its corresponding QA are set to the fill value of 32,767. Additionally, two algorithms were designed for the MODIS time-share two-stage CI. The first algorithm averages only the 8-day CIs corresponding to the QA with the highest frequency. The second algorithm averages the 8-day CIs corresponding to the QAs with frequencies less than and equal to the highest QA frequency. Finally, the two algorithms were applied to produce the global MODIS time-share two-stage CI products for 2008 and 2018, which were statistically evaluated to identify the optimal algorithm, and the time-share two-stage CIs obtained from the optimal algorithm were validated against field-measured CIs, resulting in the final MODIS time-share two-stage CI product. Furthermore, time-share two-stage CIs for each of the two vegetation cycles were derived based on the MODIS LSP product’s distinction between annual and biannual vegetation. As a result, the MODIS time-share two-stage CI product produced in this study consists of eight bands: the first to fourth bands represent the estimated time-share two-stage CIs for each phenological stage of the vegetation cycle, and the fifth to eighth bands correspond to the QA values for each stage. The band settings are provided in

Table 2. Band settings for the MODIS time-share two-stage CI product.

Band	Vegetation Cycle	Stage	Value	Range
1	1	LOS	CI	3300–10,000, 32,767
2		LFS	CI	3300–10,000, 32,767
3	2	LOS	CI	3300–10,000, 32,767
4		LFS	CI	3300–10,000, 32,767
5	1	LOS	QA	0–4, 32,767
6		LFS	QA	0–4, 32,767
7	2	LOS	QA	0–4, 32,767
8		LFS	QA	0–4, 32,767

.

4. Results

4.1. Estimation of Vegetation Phenology Parameters from CI and NDVI Time Series

Figure 4 and Figure 5 show the accuracy evaluation results, comparing the SOS and EOS estimated from the CI and NDVI time series for typical pixels with those from the MODIS LSP product. The accuracy evaluation metrics are the root mean square error (RMSE) and bias. For the CI time series, Figure 4a,b indicate that, except for PWe and CL¹, the RMSE of the estimated SOS using the optimal threshold is less than 40, with best thresholds concentrated at approximately 40% and 80%, respectively. For the EOS, except for Wsa, PWe, CL¹, and the CVM, the RMSE associated with the optimal thresholds is less than 35, and the optimal thresholds are concentrated between 80% and 90%. Figure 4c,d demonstrate that the estimated SOS using the optimal threshold shows a balance between overestimation and underestimation, where CL^2-1 and CL^2-2 represent the first and second vegetation cycles of biannual cropland, respectively, as indicated in the legend of Figure 4. For the EOS, underestimation is more common, indicating that the estimated EOS tends to appear earlier. For the NDVI time series, Figure 5a,b show that the optimal thresholds for estimating the SOS are mostly between 40% and 70%, whereas those for the EOS are mainly between 10% and 40%. In addition, Figure 5c,d indicate that the estimated SOS using the optimal threshold shows a balance between overestimation and underestimation. The bias of the EOS is relatively large, mostly reflecting overestimation. In general, the accuracy of estimating phenological parameters via the NDVI time series tends to have a small degree of uncertainty relative to that of the CI time series when the improved dynamic threshold method is used. This is probably due to CI being influenced by multiple factors during the inversion process, especially the high standards needed for accurately measuring hotspots of typical features of the BRDF shapes in the principal plane, which leads to reduced stability. However, the errors remain within acceptable limits, indicating that the CI exhibits a distinctly and approximately inverse trend of phenological variations compared with the NDVI. Additionally, the optimal thresholds for estimating the SOS and EOS across different IGBP classes based on CI and NDVI time series varied significantly. This finding indicates that the optimal thresholds for estimating phenological parameters are not only related to the phenological stages but also differ depending on the vegetation indices and land cover types.

4.2. Development of the Global MODIS Time-Share Two-Stage CI Product

4.2.1. Comparison of CIs for Different Phenological Stages

The analysis of these studies indicates that the CI exhibits distinct phenological variations, similar to those of the NDVI, although some differences between them deserve further research. Based on this finding, the global mean CIs for different phenological stages in 2020 were calculated and analyzed via ANOVA followed by Tukey’s post hoc test. The results (

Table 3. CI estimation of the different phenological stages for each IGBP class globally.

Vegetation	Estimation of CI
	LFS	LOS	LGS	LSS
ENF ***	0.549 ± 0.055 a	0.539 ± 0.038 c	0.534 ± 0.042 d	0.546 ± 0.043 b
EBF ***	0.626 ± 0.109 c	0.654 ± 0.107 b	0.654 ± 0.119 b	0.655 ± 0.104 a
DNF ***	0.604 ± 0.062 a	0.572 ± 0.027 c	0.566 ± 0.032 d	0.583 ± 0.031 b
DBF ***	0.722 ± 0.073 b	0.716 ± 0.064 c	0.711 ± 0.075 d	0.728 ± 0.069 a
MF ***	0.674 ± 0.078 d	0.683 ± 0.067 b	0.676 ± 0.078 c	0.697 ± 0.075 a
CSh ***	0.833 ± 0.083 a	0.801 ± 0.080 c	0.790 ± 0.085 d	0.819 ± 0.085 b
OSh ***	0.854 ± 0.044 a	0.803 ± 0.061 c	0.791 ± 0.073 d	0.817 ± 0.065 b
Wsa ***	0.723 ± 0.071 b	0.722 ± 0.073 c	0.715 ± 0.084 d	0.734 ± 0.077 a
Sav ***	0.765 ± 0.068 a	0.754 ± 0.071 c	0.749 ± 0.083 d	0.762 ± 0.072 b
GL ***	0.810 ± 0.074 a	0.789 ± 0.080 c	0.783 ± 0.086 d	0.797 ± 0.082 b
PWe ***	0.778 ± 0.065 b	0.776 ± 0.077 c	0.775 ± 0.085 d	0.779 ± 0.080 a
CL ***	0.809 ± 0.063 a	0.792 ± 0.059 c	0.787 ± 0.069 d	0.798 ± 0.064 b
UB ***	0.792 ± 0.070 ad	0.787 ± 0.069 b	0.784 ± 0.076 c	0.792 ± 0.072 a
CVM ***	0.780 ± 0.079 b	0.775 ± 0.072 c	0.771 ± 0.081 d	0.781 ± 0.074 a

Different letters in the same row indicate significant differences between the mean CIs of different phenological stages. (Tukey’s test, p < 0.05). LFS: leaf-off season; LOS: leaf-on season; LGS: leaf-gathering season; LSS: leaf-scattering season. “*” indicates significant differences in the mean CIs for the same IGBP classes and different phenological stages (* p < 0.05, ** p < 0.01, *** p < 0.001).

) indicated that there were highly significant differences (p < 0.001) in the CIs across different phenological stages within the same IGBP class globally. Tukey’s post hoc test further revealed that, except for CI_LOS (the CI in the LOS) and CI_LGS (the CI in the LGS) in EBF and CI_LFS (the CI in the LFS) and CI_LSS (the CI in the LSS) in UB, the differences in CIs between phenological stages in all other IGBP classes were highly significant (p < 0.001). This confirms once again that CI exhibits significant seasonal variations. In addition, except for EBF and MF, where CI_LFS is smaller than CI_LOS, the other IGBP classes presented CI_LFS values greater than CI_LOS values, with CI_LGS being the smallest. For EBF, which is mainly distributed in the tropics, seasonal variations are less pronounced, and frequent cloud cover often interferes with CI inversion, possibly explaining why the CI_LFS is smaller than the CI_LOS. According to the IGBP classification system, mixed forests (MFs) are defined as combinations of deciduous needleleaf forests (DNFs), evergreen needleleaf forests (ENFs), deciduous broadleaf forests (DBFs), and evergreen broadleaf forests (EBFs) [54]. However, in the inversion process, the CI of the MF did not account for the differences between needleleaf and broadleaf forests in the mixed pixels. Instead, the broadleaf forest CI-NDHD equation was directly applied for estimation, neglecting the influence of needleleaf forests in MF pixels [66]. This resulted in higher CIs, particularly for LOS, leading to the observed pattern of CI_LFS being smaller than the CI_LOS for MF. Li applied two or three CIs of LFS, LOS, LGS, and LSS to light use efficiency (LUE) models. The results demonstrated that the improvement in the estimated GPP depending on two CIs (8.97%) of the leaf-on season and the leaf-off season is nearly as high as the estimated GPP that relies on three CIs (9.76%) of the leaf-scattering, leaf-gathering, and leaf-off seasons, which may be due to the rapid growth of leaves at the beginning of the growing season, which has already developed the appearance of the canopy. Consequently, the CIs of the LGS and LSS exhibited only minor differences [65]. This is consistent with the results of our study. Therefore, based on the impact of different numbers of time-share two-stage CIs on improving GPP estimation accuracy, we decided to divide a vegetation cycle into two stages, the leaf-on season (LOS) and the leaf-off season (LFS), to make the seasonal variations in the CI more concise and effective for various meteorological and ecological models. A global MODIS time-share two-stage CI product was developed, considering the quality assurance band of the MODIS 8-day CI product.

4.2.2. Comparison of Different Developed Algorithms

Table 4. Percentages of CI_LOS < CI_LFS and CI_LOS > CI_LFS in 2008 and 2018.

Year	Range	The First Vegetation Cycle				The Second Vegetation Cycle
		The First Algorithm		The Second Algorithm		The First Algorithm		The Second Algorithm
		CILOS < CILFS	CILOS > CILFS	CILOS < CILFS	CILOS > CILFS	CILOS < CILFS	CILOS > CILFS	CILOS < CILFS	CILOS > CILFS
2008	NH	66.08%	33.92%	68.57%	37.43%	50.34%	49.66%	51.69%	48.31%
	Trop	65.37%	34.63%	66.73%	33.27%	58.31%	41.69%	59.45%	40.55%
	SH	73.76%	26.24%	74.15%	25.85%	56.93%	43.07%	57.16%	42.84%
	Global	66.22%	33.78%	68.23%	31.77%	55.04%	44.96%	56.14%	43.86%
2018	NH	58.16%	41.84%	60.11%	39.89%	46.83%	53.17%	47.32%	52.68%
	Trop	62.63%	37.37%	64.57%	35.43%	53.95%	46.05%	55.07%	44.93%
	SH	68.26%	31.74%	68.66%	31.34%	66.87%	33.13%	67.35%	32.65%
	Global	60.10%	39.90%	61.98%	38.02%	52.91%	47.09%	53.71%	46.29%

NH (Northern Hemisphere, 23.5–90°N), SH (Southern Hemisphere, 23.5–90°S), and Trop (Tropics, 23.5°S–23.5°N); CI_LOS (time-share CI in LOS), CI_LFS (time-share CI in LFS).

shows that the percentage of CI_LOS < CI_LFS in the second algorithm is greater than that in the first algorithm for both the first and second vegetation cycles. Since the leaves exhibit a clumped state during the LOS, CI_LOS should be smaller than CI_LFS. This suggests that the second algorithm generally performs better, capturing the seasonal variations in CI more accurately across the southern and northern hemispheres, as well as the tropics. Therefore, this study selected the second algorithm to generate a global MODIS time-share two-stage CI product. According to the second algorithm, the percentage of CI_LOS < CI_LFS during the first vegetation cycle shows little variation between the southern and northern hemispheres and the tropics, with all values exceeding 60%. In 2008 and 2018, the global percentages were 68.23% and 61.98%, respectively. During the second vegetation cycle, except for the northern hemisphere in 2018, the percentage of CI_LOS < CI_LFS exceeded 50%, with global percentages of 56.14% and 53.71% in 2008 and 2018, respectively. This suggests that CI typically displays seasonal variation, but its ability to capture seasonal variation is somewhat reduced during the second vegetation cycle, particularly for biannual vegetation. Additionally, during the first vegetation cycle, the southern hemisphere had the highest percentage of CI_LOS < CI_LFS, with 74.15% in 2008 and 68.66% in 2018, significantly surpassing the percentages in the tropics and the Northern Hemisphere. This suggests that the seasonal variation in CI is most pronounced in the southern hemisphere, likely because of its significantly lower number of vegetation-covered pixels compared to the northern hemisphere and the tropics. The southern hemisphere is less influenced by topographic diversity, human activities, and latitudinal variation. In contrast, needleleaf forests, which exhibit relatively stable temporal trends due to their lower renewal frequency and are predominantly located in the Northern Hemisphere (~96%) [35,67], have contributed to the reduction in the percentage of CI_LOS < CI_LFS in the northern hemisphere.

4.2.3. Accuracy Evaluation for the MODIS Time-Share Two-Stage CI Product

Figure 6 shows the accuracy evaluation results when the time-share two-stage CIs are compared with the field-measured CIs. A total of 95 datasets were involved in the evaluation, and the time-share two-stage CIs and the field-measured CIs were roughly in agreement (RMSE = 0.06, bias = 0.01). Among them, nine data points were overestimated by more than 0.1, and six data points were underestimated by more than 0.1, all of which are marked with gray circles. These outliers accounted for only 15.6% of the total data. The values in parentheses in Figure 6 represent the accuracy evaluation results after removing the gray-marked data, showing a significant improvement in accuracy (RMSE = 0.04, bias = 0.00). Generally, the MODIS time-share two-stage CI product shows high accuracy but has a slight tendency toward overestimation.

Figure 7 shows the interannual variations in the time-share two-stage Cis and field-measured CIs at several sites. Across the four sites, the time-share two-stage CIs exhibited obvious CI_LOS < CI_LFS rules, further confirming that the CI can effectively capture the seasonal variation in vegetation. Of these, the SCBI site, classified as DBF, displayed the most pronounced difference between the CI_LOS and CI_LFS, with the CI_LOS consistently being smaller than the CI_LFS throughout the 2001–2020 time series. The field-measured CIs are also closely aligned with the time-share two-stage CIs. Additionally, while there are some discrepancies between the field-measured CIs and the time-share two-stage CIs, these differences are mostly within one standard deviation, suggesting that the errors fall within an acceptable range. These errors may stem from the pooled field-measured CI data and uncertainties in the CI inversion process. Even though we selected field-measured CIs covering more than 3 months within the LOS or LFS, errors persist due to the large scale and the possibility that field-measured CIs may not be continuous data on a monthly scale.

4.3. Spatial Distribution and Temporal Variation in the Global Time-Share Two-Stage CI

Figure 8 shows the global distributions of the multiyear average time-share two-stage CIs for both the LOS and LFS during the first vegetation cycle. In general, the spatial distributions of the global time-share two-stage CIs for both the LOS and LFS shows a similar pattern to that of IGBP classes (Figure 1), with lower CI values for forests and higher CI values for grasslands. Needleleaf forests in North America and Eurasia are distributed near 60°N latitude, and evergreen broadleaf forests in tropical rainforest regions present lower CIs, generally below 0.65. Grasslands in Inner Mongolia, China, the Midwestern United States, and equatorial Africa display higher CIs, approximately 0.78. Open shrublands in southern Africa and northern Eurasia also presented higher CIs, typically above 0.80. Therefore, the CI is not only related to phenological stage but also strongly correlated with land cover type and spatial distribution. Additionally, a comparison between Figure 8a and Figure 8b reveals that the CI_LOS is generally smaller than the CI_LFS worldwide, indicating pronounced seasonal variation in the CI. This further supports the conclusion that CI is lower and that the clumping effect is stronger during LOS, whereas CI is greater and the clumping effect is weaker during LFS.

Figure 9 compares the multiyear average CI_LOS and CI_LFS across different land cover types from 2001 to 2020. Generally, the CI_LOS is smaller than the CI_LFS across all land cover types, which is consistent with the characteristics shown in Figure 8, suggesting that CIs across all land cover types exhibit significant seasonal variations at the global scale. Among these forests, needleleaf forests presented the lowest CI values, with minimal differences between the CI_LOS and CI_LFS, indicating a high degree of clumping during both the LOS and LFS. The CI values of broadleaf forests are slightly greater than those of needleleaf forests, with a greater difference between the CI_LOS and CI_LFS, indicating more distinct seasonal variations. This is because the clustered structure of needleleaf forests leads to a stronger clumping effect, and needle leaves have a lower renewal frequency [35,67], resulting in more conservative seasonal variations in CI. Osh has the highest CI values, with a multiyear average CI_LOS of 0.80 and CI_LFS of 0.86. CSh, GL, PWe, CL, UB, and CVM also presented relatively high CI values, whereas broadleaf forests, MF, Wsa, and Sav presented relatively low CI values.

As shown in Figure 10a,b, the QA distribution of the time-share two-stage CI product during LOS is better than that during LFS, with QA predominantly at 0 and 1, accounting for over 90%, which is significantly greater than that during LFS (~60%). During LFS, the QA of the time-share two-stage CI product is lower, and the time-share two-stage CIs with low QA are primarily distributed in high-latitude regions (50–90°N), with a significant presence of fill values. In addition, a considerable number of fill values are present in the mid-to-low latitude regions (−30° to 30°N) (Figure 10a,b), primarily due to the contamination of clouds and aerosols [47]. Figure 10c,d illustrate that the QA of the time-share two-stage CI product from 2001 to 2020 was relatively consistent. During LOS, except for a slightly lower proportion of QA = 0 in 2001–2002, the percentage remained at approximately 73% in the following years, whereas the proportion of QA values between 3 and 4 remained consistently low, at approximately 8%. During LFS, the quality of the time-share two-stage CI product showed a notable decline, with the proportion of QA = 0–1 remaining at approximately 60%, whereas the proportion of QA between 3 and 4 increased to approximately 35%. This is attributed primarily to the shorter daylight hours in high latitudes of the Northern Hemisphere (50–90°N) during the LFS, which results in fewer observation angles, affects the quality of multiangle input data and leads to relatively lower overall data accuracy.

5. Discussion

5.1. Seasonal Variability in the Global CI

Quantitative research on CI started relatively late, resulting in considerable debate concerning conclusions concerning its seasonal variations. To address these issues, this study focuses on resolving two current controversies regarding CI: (1) whether CI exhibits seasonal variations, and (2) based on the conclusions, MODIS time-share two-stage CI products were developed to further explore the specific seasonal variation characteristics of the CI.

The results indicate that although constrained by the inversion process, the CI has a large uncertainty relative to the NDVI in capturing the precise growth status of vegetation. However, CI still exhibits more pronounced seasonal variations, aligning with previous findings [27,35,47,68,69,70], making it a valuable parameter for phenological monitoring. Building on these conclusions and those of previous studies [65], this study used the MODIS LSP product to divide the vegetation cycle into LOS and LFS, taking into account both the seasonal variation characteristics of CI and reducing the influence of low-quality data and background reflectance changes through quality screening and averaging. The global 500 m MODIS time-share two-stage CI product from 2001 to 2020 was then derived from the MODIS 8-day CI product. This approach provides reliable CI values with seasonal variations for global- and regional-scale modeling of ecological, meteorological, and other surface processes, ultimately improving the accuracy of various models. Globally, the CIs of different land cover types exhibit a strong clumping effect during the LOS and a weak clumping effect during the LFS. Needleleaf forests have the lowest CI values, with a slight difference between the CI_LOS and CI_LFS. Broadleaf forests have slightly higher CI values than do needleleaf forests, with more distinct seasonal variations. Shrublands presented the largest CI values, with the most pronounced seasonal variations. Thus, the seasonal variations in the CI are sensitive to land cover type, suggesting that the CI may assist in the accurate classification of land cover.

Additionally, some scholars [71,72] tend to assume that the BRDF varies slowly over time and even consider it constant throughout the year in their applications. In contrast, the MODIS 8-day CI product, used in this study to derive the MODIS time-share two-stage CI product, is produced based on the NDHD–CI relationship. Since NDHD is calculated via reflectance in the directions of hotspots and darkspots—two typical features of BRDF shapes in the principal plane—the results of this study imply that the BRDF exhibits seasonal variations.

Moreover, there are still some limitations when validating whether CI exhibits seasonal variation characteristics. Ground-based phenological observation data provides accuracy and objectivity and is an important part of the validation method, ensuring the reliability of the conclusions. However, due to various constraints, ground-based phenological observation data cannot be applied on a large scale, and continuous long-term ground-based phenological observation data is difficult to obtain. Furthermore, it is challenging to find ground-based phenological observation data that matches the typical pixels selected in this study under conditions such as land cover restrictions. Therefore, this study uses the MODIS LSP product for indirect validation. Studies have shown that this product can effectively monitor vegetation phenology and is widely used as a direct application or auxiliary validation data in various phenology monitoring studies [1,2]. Nevertheless, we hope that in the future, the results obtained from CI can be compared with ground-based observed phenological data to better understand and validate the seasonal variation characteristics of CI.

5.2. Uncertainty of the MODIS Time-Share Two-Stage CI Product

In this study, the MODIS time-share two-stage CI product is derived from the MODIS 8-day CI product, which is strongly influenced by the quality of the input data (MCD43A1/A2 and MCD12Q1) and the inversion algorithm. The MODIS 8-day CI product used in this study is the C6 version, which uses the enhanced and more accurate C6 MCD43 product compared with the C5 version [73], so the C6 version of the MODIS 8-day CI product has a greater proportion of inversion of the main algorithm and a lower proportion of filling values, which determines the high accuracy of the MODIS time-share two-stage CI product to some extent. However, since the basic structural characteristics of the vegetation canopy are determined by the IGBP classification schemes in MCD12Q1 as input data, the error of the land cover types also affects the quality of the MODIS 8-day CI product [23]. Furthermore, the quality of MODIS CI products is influenced by factors such as sparse vegetation cover, high background reflectance [20,23,35], and surface shading due to topographic relief [20,21,24]. Therefore, special attention should be given to these factors when MODIS time-share two-stage CI products are used.

In the current version, the MODIS time-share two-stage CI product is related mainly to the MODIS LSP (MCD12Q2) product in the vegetation phenological cycle, which divides a vegetation cycle into LOS and LFS using the “greenup” and “dormancy” bands. The MCD12Q2 product estimates phenological parameters based on the 2-band enhanced vegetation index (EVI2). In theory, the CI characterizes the degree of clumping in the vegetation canopy; however, the EVI/NDVI mainly characterize the changes in the greenness of vegetation leaves [74,75]. Despite the probable differences in the vegetation phenological cycles identified in this study, the assumption that the foliage clumping state indicated by the CI is synchronously related to the foliage greenness indicated by the EVI/NDVI may unavoidably generate uncertainties to some degree for the MODIS time-sharing two-stage CI product [36]. Previous studies have analyzed the timing of maximum clumping and maximum greenness of vegetation leaves, indicating that the timing of maximum clumping generally precedes that of maximum greenness by approximately 0.73 months (22 days) in the Northern Hemisphere [36]. This implies that such a difference has a marginal influence on the uncertainties of this time-share two-stage CI product [37,38].

Although the MCD12Q2 product derives vegetation phenological parameters based on the EVI, which represents the greenness of vegetation leaves, in this study, we selected the MODIS CI product with an 8-day temporal resolution as the input data. Therefore, when phenological stages are divided based on MCD12Q2, the 22-day time difference between maximum aggregation and maximum greenness has a negligible effect on the accuracy of the MODIS time-share two-stage CI product. However, we will also aim to explore more detailed seasonal variation characteristics of CI based on this conclusion in future research.

In addition, the filling values in the MCD12Q2 product may generate uncertainties, especially for the filling values in the MODIS time-share two-stage CI product. Owing to technical constraints, retrieving the values of MCD12Q2 in arid and semiarid ecosystems has faced challenges [76,77,78], especially in Australia, where vegetation phenology is highly correlated with rainfall [78,79]. This has resulted in substantial missing remote sensing information on vegetation phenology in this region, leading to numerous filling values in the MODIS time-share two-stage CI product. Indeed, how to consider and reduce the uncertainties generated by the upper-level input product is a challenge for all remote sensing quantitative products.

6. Conclusions

As a significant vegetation structural parameter, the seasonal variation in CI plays an important role in LAI mapping and various ecosystem models. However, a phenologically simplified time-sharing two-stage CI product helps to further extend the CI applications mentioned above. In this study, we estimated the key phenological parameters from CI and NDVI time series as a comparison based on the discrete Fourier transform and the improved dynamic threshold method, respectively, and confirmed that the CI exhibits seasonal variations. Additionally, we obtained the SOS and EOS on a global scale via a pixel-by-pixel approach from the MODIS LSP product (MCD12Q2, V061) and developed the MODIS time-share two-stage CI product based on the quality bands (QA) of the MODIS 8-day CI product. The main conclusions are as follows:

(1)

The study shows that the CI exhibits an approximately inverse trend of phenological variation compared with the NDVI. The optimal thresholds for estimating phenological parameters are not only related to phenological stages but also vary with vegetation indices and land cover types. Additionally, the optimal thresholds for the SOS generally range from 40% to 80%, whereas those for the EOS range from 80% to 90%, with errors remaining within acceptable limits.

(2)

The accuracy evaluation results of the MODIS time-share two-stage CIs using the field-measured CIs indicate that the time-share two-stage CI is highly accurate (RMSE = 0.06, bias = 0.01), although a slight overestimation is generally observed.

(3)

Globally, based on the MODIS time-share two-stage CI product, the CI generally shows distinct seasonal variations, with the CI_LOS being smaller than the CI_LFS across all land cover types. Needleleaf forests presented the smallest CI values, and the difference between the CI_LOS and CI_LFS was insignificant. Compared with needleleaf forests, broadleaf forests, with slightly higher CI values, present more pronounced seasonal variations. In addition, Osh presented the largest CI values, whereas CSh, GL, PWe, CL, UB, and CVM presented relatively high CI values. Broadleaf forests, MF, Wsa, and Sav had relatively low CI values.

(4)

Compared with the LFS stage, the quality of the MODIS time-share two-stage CI product is better in the LOS stage, where the QA values are basically 0 and 1, accounting for more than 90% of the total, which is significantly greater than that in the LFS stage (~60%).

As research has advanced, the seasonal variation characteristics of CI have been realized and gradually clarified, which should not be overlooked in studies, especially from subsequent CI applications. The MODIS time-share two-stage CI product in this study provides seasonal, large-scale, pixel-by-pixel CI values for reference and comparison in potential studies of CI derived from VIIRS and other data. The product presents the seasonally simplified variations in CI in a concise fashion, providing new insights for the derivation of CI phenology and opportunities for subsequent CI applications. In addition, exploring the extent to which CI characterizes the relationship between leaves and other canopy structures in a localized area, as well as investigating more detailed seasonal variation characteristics of CI, will be valuable in the future.