Mapping global leaf inclination angle (LIA) based on field measurement data

Leaf inclination angle (LIA), the angle between the leaf surface normal and zenith directions, is a vital trait in radiative transfer, rainfall interception, evapotranspiration, photosynthesis, and hydrological processes. Due to the difficulty of obtaining large-scale field measurement data, LIA is typically assumed to follow the spherical leaf distribution or simply considered to be constant for different plant types. However, the appropriateness of these simplifications and the global LIA distribution are still unknown. This study compiled global LIA measurements and generated the first global 500 m mean LIA (MLA) product by gap-filling the LIA measurement data using a random forest regressor. Different generation strategies were employed for noncrops and crops. The MLA product was evaluated by validating the nadir leaf projection function (G(0)) derived from the MLA product with high-resolution reference data. The global MLA is 41.47°±9.55°, and the value increases with latitude. The MLAs for different vegetation types follow the order of cereal crops (54.65°) > broadleaf crops (52.35°) > deciduous needleleaf forest (50.05°) > shrubland (49.23°) > evergreen needleleaf forest (47.13°) ≈ grassland (47.12°) > deciduous broadleaf forest (41.23°) > evergreen broadleaf forest (34.40°). Cross-validation shows that the predicted MLA presents a medium consistency (r=0.75, RMSE = 7.15°) with the validation samples for noncrops, whereas crops show relatively lower correspondence (r=0.48 and 0.60 for broadleaf crops and cereal crops, respectively) because of the limited LIA measurements and strong seasonality. The global mean G(0) is 0.68±0.11. The global G(0) distribution is out of phase with that of the MLA and agrees moderately with the reference data (r=0.62, RMSE = 0.15). This study shows that the common spherical and constant LIA assumptions may underestimate the interception of most vegetation types. The MLA and G(0) products derived in this study could enhance our knowledge of global LIA and should greatly facilitate remote sensing retrieval and land surface modeling studies.

The global MLA and G(0) products can be accessed at https://doi.org/10.5281/zenodo.12739662 (Li and Fang, 2025).

1 Introduction

Vegetation regulates terrestrial carbon and water cycles through a series of biophysical processes such as photosynthesis, respiration, and transpiration (Foley et al., 1996; Chen et al., 2019). These biophysical processes are mainly carried out by leaves, and the characterization of leaves within canopies is vital for remote sensing and Earth system modeling (Ross, 1975; Lawrence et al., 2019). Leaf inclination angle (LIA) denotes the inclination of the leaf or needle to the horizontal plane or the angle between the leaf surface normal and zenith directions (Wilson, 1960). LIA is a key canopy structural trait that determines radiative transfer, rainfall interception, evapotranspiration, photosynthesis, and hydrological processes (Sellers, 1985; Ross, 1981; Mantilla-Perez and Salas Fernandez, 2017; Xiao et al., 2000; Maes and Steppe, 2012). LIA has been commonly used in radiative transfer modeling (RTM), remote sensing inversion, and land surface modeling (LSM) studies (Tang et al., 2016; Wang and Fang, 2020; Lawrence et al., 2019; Ross, 1975).

At the canopy scale, the probability density of LIA or the fraction of leaf area per unit LIA is expressed as the leaf angle distribution (LAD) (de Wit, 1965). De Wit (1965) summarized six theoretical LADs, including planophile, erectophile, extremophile, plagiophile, uniform, and spherical distributions. Specifically, the spherical distribution assumes that the relative probability density of the LIA is proportional to the area of the corresponding sphere surface element and its mean leaf inclination angle (MLA) equals 57.3° (MLA = 57.3°) (de Wit, 1965). Furthermore, Ross (1981) defined the inclination index (χ_L) to describe the departure of LAD from the spherical distribution. For the planophile distribution, χ_L=1; for the erectophile distribution, ${\mathit{\chi }}_{\mathrm{L}}=-\mathrm{1}$; and for the spherical distribution, χ_L=0. In the radiative transfer regime, LIA is generally represented by the leaf projection function (G(θ)), which is defined as the average projection ratio of unit leaf area in the illumination or viewing direction θ (Ross, 1981; Nilson, 1971). The spherical distribution is characterized by an isotropic leaf projection function (G≡0.5) (de Wit, 1965).

In the field, LIA can be measured directly based on the leaf's geometrical structure or using indirect optical methods (Lang, 1973; Ryu et al., 2010; Norman and Campbell, 1989; Weiss and Baret, 2017). Using these methods, several LIA measurements have been carried out and some LIA datasets were constructed (Kattge et al., 2020; Chianucci et al., 2018; Hinojo-Hinojo and Goulden, 2020; Pisek and Adamson, 2020). These field methods are usually time-consuming and labor-intensive, and it is typically difficult to acquire large-scale LIA (Li et al., 2023). In addition, the existing LIA datasets have not been comprehensively analyzed. LIA has also been estimated from satellite imagery through empirical relationships or radiative transfer model inversions (Zou and Mõttus, 2015; Bayat et al., 2018; Goel and Thompson, 1984). Remote sensing methods are used primarily for crops in local regions, and the generality of these algorithms is limited (Li et al., 2023). Due to the difficulty of large-scale LIA measurements and estimations, our knowledge of the global LIA remains lacking.

Because our understanding of the global LIA is limited, different LIA simplification strategies have been adopted in various studies. For example, LIA is typically assumed to follow a spherical distribution (Tang et al., 2016; Zhao et al., 2020; Wang and Fang, 2020). However, this assumption may decrease the accuracy of radiative transfer modeling, significantly underestimate the radiation interception (Stadt and Lieffers, 2000), and cause large errors (>50 %) in leaf area index (LAI) measurements and inversions (Yan et al., 2021). The spherical LIA assumption may introduce greater error in the nadir direction than inother viewing geometries (Yan et al., 2021), considering the large G variation in this direction (Wilson, 1959). The lack of global LIA knowledge also limits the retrieval of other vegetation structural parameters (Li et al., 2023). In many LSMs, LIA is commonly treated as a fixed value for different plant functional types (PFTs) (Lawrence et al., 2019; Majasalmi and Bright, 2019). Field LIA measurements have demonstrated that the spherical distribution is not appropriate for forests, and the PFT-dependent LIA ignores LIA variation within the PFT (Pisek et al., 2013; Yan et al., 2021; Majasalmi and Bright, 2019).

This study aims to generate the first global MLA map from existing LIA field measurements using a data-driven gap-filling method. This method involves spatial expansion and upscaling of LIA measurements, as well as a random forest regressor using input spectral, climate, and PFT data. Based on the global MLA map, we tested whether the spherical LIA assumption is appropriate at the global scale. The new MLA map was validated by comparing the nadir G (G(0)) derived from the MLA with high-resolution reference data. Section 2 outlines the materials and methods employed to generate and evaluate the global MLA. Section 3 presents the global LIA measurements, global MLA and G(0), and evaluation results. Section 4 discusses the performance of the global MLA and G(0), the usage of the new MLA map, and the limitations of the study. Section 6 presents the main conclusions.

2 Materials and methods

2.1 LIA measurement data

2.1.1 TRY LIA dataset

TRY is a network of vegetation scientists headed by Future Earth, the Max Planck Institute for Biogeochemistry, and the German Centre for Integrative Biodiversity Research, providing a global database of curated plant traits (the TRY database) (https://www.try-db.org/TryWeb/Home.php, last access: 17 April 2023). Since its establishment in 2007, the TRY database has continuously evolved and has become one of the most widely used vegetation trait databases. The latest V6 version (released on 13 October 2022) employed in this study contains 15 409 681 trait records covering 305 594 plant taxa (Kattge et al., 2020). In this database, LIA was recorded as a numerical or categorical variable. After data extraction and checking, 31 043 valid records were used, which include numerical LIA, locations, and species. Many measurements lack location information, whereas, for some locations, there are many measurements for individual species. The spatial distribution map appears to be relatively sparse despite a large volume of data (Fig. 1). The LIA measurements of most species are located in the Northern Hemisphere, while the LIA measurements in South America are mainly from palms.

Figure 1 Locations of global leaf inclination angle measurements collected from TRY and the literature. DBF: deciduous broadleaf forest, DNF: deciduous needleleaf forest, EBF: evergreen broadleaf forest, ENF: evergreen needleleaf forest, CRO-B: broadleaf crops, CRO-C: cereal crops, GRA: grassland, SHR: shrubland.

2.1.2 LIA data from the literature

To fully utilize distributed and considerable LIA measurement data in the published literature, several keyword searches (leaf angle, leaf inclination angle, and leaf tilt angle) were performed in the Web of Science, Google Scholar, Google, and Chinese documentary databases. Subsequently, the LIA, location, and species information were manually extracted from the literature (Fig. 1). Several LIA measurements were already included in the TRY database (Chianucci et al., 2018; Pisek and Adamson, 2020). After aggregating LIA measurements for the same species at the same location, 780 LIA records were accessed from previous studies (Hinojo-Hinojo and Goulden, 2020; Pisek et al., 2022; Chen et al., 2021).

2.1.3 Manual LIA extraction

After TRY and a literature search, only a few measurements in the northern tundra region were obtained, and the measurements in tropical regions are dominated by palm trees (Fig. 1). Therefore, LIA data for the northern tundra and tropical regions were further extracted from horizontal side-view photographs searched from Google (Fig. S1 in the Supplement). ImageJ software (https://imagej.nih.gov/ij/, last access: 17 April 2023) was used to process the leveled photographs and derive LIA following the method of Pisek et al. (2011). The TRY species location data (848 919, Fig. S2b) (3 January 2022) were used to obtain the dominant species information in tropical rainforests and the northern tundra. The species location points in these two vegetation types were spatially filtered and the frequency of occurrence for each species was counted. The species with a high frequency of occurrence were selected to measure the LIA. For each species, more than 75 leaves perpendicular to the viewing direction were selected and processed based on visual judgment to ensure the stability and reliability of the MLA (Pisek et al., 2013). In total, the MLA of 104 species was manually derived.

In this study, most LIA measurements were obtained with a protractor and level digital photogrammetry, especially for needleleaf species. Therefore, the distinction between branches and leaves was considered. The diverse LIA records from different sources were sorted to match the TRY species and to get the PFT based on the TRY Categorical Traits Dataset 2018 (https://www.try-db.org/TryWeb/Data.php#3, last access: 17 April 2023). LIA measurements from different sources were unified into canopy-level MLA with average operation by leaf number (see Appendix A). If there were multiple LIA records for the same species, the mean value was computed for the same location and species. In total, 5554 LIA records of 1194 species were collected, covering the growing season from 2001 to 2022. LIA location replicates per species range from 1 to 330, and there are fewer than 50 location replicates for most species (98 %). Considering the different numbers of records for each species, the LIA data were further aggregated by species.

2.2 Remote sensing data

2.2.1 Ancillary data used for MLA mapping

The ancillary data used for global MLA mapping and analysis are listed in Table 1. Most Earth observation data were accessed and processed in Google Earth Engine (GEE) (https://earthengine.google.com/, last access: 12 January 2025). The PFT classification system in the MODIS global 500 m land cover type product (MCD12Q1 C6) was used and mode-aggregated from 2001 to 2022 to match the LIA measurements (Fig. S3) (Sulla-Menashe et al., 2019). The 2001–2022 Landsat surface reflectance (Level 2, Collection 2, Tier 1) (Crawford et al., 2023), including Landsat 5 (2001–2012), Landsat 7 (2012–2013), and Landsat 8 (2013–2022), was utilized to generate a global 30 m PFT map (Sect. 2.3.1), which was subsequently employed for LIA upscaling. Considering the sensitivity of directional reflectance variation to LIA (Jacquemoud et al., 2009; Li et al., 2023), the 2001–2022 MODIS bidirectional reflectance distribution function (BRDF) model parameter dataset (MCD43A1 C6.1) (Schaaf and Wang, 2015) and nadir BRDF-adjusted reflectance dataset (MCD43A4 V6 NBAR) (Schaaf and Wang, 2015) produced daily using 16 d of Terra and Aqua MODIS data at 500 m resolution were utilized as predictive variables. We used MCD43A1 C6.1 as well as MCD12Q1 and MCD43A4 C6 for MLA mapping as these data were available on GEE when this study was conducted. For MCD43A1, only minor calibration changes and polarization correction were adopted in the upgrading from Collection 6 to 6.1, while the MCD12Q1 and MCD43A4 algorithms remain the same (https://landweb.modaps.eosdis.nasa.gov/data/userguide/MODIS_Land_C61_Changes.pdf, last access: 12 January 2025). In addition, the multiyear aggregation of these products (Table 2) further mitigates the version impact. Due to the scarcity of crop LIAs and the lack of location information for existing crop LIA measurements, fine-resolution ($\mathrm{10}/\mathrm{30}$ m) crop-type maps (Table 1) in 2018 were employed to support crop LIA mapping. Other data include the ERA5-Land reanalysis data, the ALOS digital elevation model (AW3D30 V3.2), and the 2001–2022 MODIS LAI product (MCD15A2H) (Myneni et al., 2015). The LAI product was averaged and aggregated from 2001–2022.

Table 1 Remote sensing data for global MLA mapping. BRDF: bidirectional reflectance distribution function.

Table 2 Predictive features in global MLA mapping.

2.2.2 High-resolution reference data

The high-resolution reference datasets provided by Ground Based Observations for Validation (GBOV; https://land.copernicus.eu/global/gbov/dataaccessLP/, last access: 17 April 2023) and DIRECT 2.1 (https://calvalportal.ceos.org/lpv-direct-v2.1, last access: 17 April 2023) were used to evaluate the generated global MLA (Fig. 2). These datasets provide high-resolution ($\mathrm{20}/\mathrm{30}$ m) LAI, effective LAI (LAI_e), and fractional vegetation cover (FVC, the proportion of the vertical projection area covered by green vegetation; Gitelson et al., 2002) data over a 3 km × 3 km area centered on each site generated using empirical relationships between various vegetation indices and ground measurements (Li et al., 2022; Brown et al., 2020). GBOV has provided continuous high-resolution reference data since 2013 (Fig. 2).

Figure 2 Locations of GBOV and DIRECT 2.1 sites used in this study. CRO: cultivated crops, MF: mixed forest, PAS: pasture/hay, WET: woody wetlands. See Fig. 1 for other acronyms. The black frame indicates sites with more than five continuous records.

The global MLA map was indirectly evaluated by comparing the nadir leaf projection function derived from MLA with reference G(0) because of the lack of high-resolution reference MLA. The high-resolution reference G(0) was derived from high-resolution LAI, FVC, and the clumping index (CI) (LAIe $/$ LAI) with the Beer–Lambert law (Fig. S4) (Nilson, 1971):

$$\begin{array}{}\text{(1)}& P\left(\mathit{\theta }\right)={\exp }^{-\frac{G\left(\mathit{\theta }\right)\cdot \mathrm{LAI}\cdot \mathrm{CI}\left(\mathit{\theta }\right)}{\cos \left(\mathit{\theta }\right)}},\end{array}$$

where P(θ), CI(θ), and G(θ) denote the gap fraction, CI, and G in direction θ, respectively. Specifically, the gap fraction in the nadir direction can be expressed by FVC.

$$\begin{array}{}\text{(2)}& P\left(\mathrm{0}\right)=\mathrm{1}-\mathrm{FVC}\end{array}$$

Therefore, the reference G(0) was derived using the following formula:

$$\begin{array}{}\text{(3)}& G\left(\mathrm{0}\right)\mathit{\_}\mathrm{CI}\left(\mathrm{0}\right)=-{\displaystyle \frac{{\ln }^{(\mathrm{1}-\mathrm{FVC})}}{\mathrm{CI}\left(\mathrm{0}\right)\cdot \mathrm{LAI}}}.\end{array}$$

By using the whole CI as the nadir CI (CI(0)) in the above equation (Fang et al., 2021; Li et al., 2022), G(0) was calculated as follows:

$$\begin{array}{}\text{(4)}& G\left(\mathrm{0}\right)\mathit{\_}\mathrm{CI}\approx -{\displaystyle \frac{{\ln }^{(\mathrm{1}-\mathrm{FVC})}}{\mathrm{CI}\cdot \mathrm{LAI}}}.\end{array}$$

2.3 Mapping global LIA

2.3.1 Data preparation

Many studies have treated LIA as a species-specific static trait and ignored within-species variations when LIA measurements are limited (Pisek et al., 2022; Toda et al., 2022; Raabe et al., 2015). Following this rationale, the spatial coverage of LIA measurements was expanded, and records without location information were also utilized (Sect. 2.1.1). Under this assumption, the LIA measurements were expanded through TRY species location data with species name matching. The species location data comprise trait measurements for common species representing an area of hundreds of square meters around the location. The dominant species was artificially identified by investigators and thus the spatial representativeness was considered. When a species had multiple LIA observations at different locations, the nearest LIA was assigned to the TRY species location. Visual inspections were conducted to remove potential TRY location biases, especially for non-vegetated points such as water bodies and deserts. After spatial expansion, the number of LIAs reached 12 328 and the spatial distribution became more uniform (Fig. S2c).

Figure 3 Leaf inclination angle (LIA) upscaling (a) and global mean LIA (MLA) mapping (b) strategies.

In this study, the scale gap between field measurements and satellite remote sensing data was fully considered. The canopy-level MLA measurement was regarded as equal to 30 m MLA considering its spatial representativeness. To upscale the MLA measurements from the canopy level to the satellite resolution (500 m), a 30 m PFT map was first derived from Landsat reflectance using a random forest classification method. The random forest was trained at a 500 m scale using the mode-aggregated MODIS PFT classification map as training samples to generate a 30 m PFT map by hierarchically selecting homogeneous pixels (with a coefficient of variation in reflectance < 0.2). The classification features were the same as those in the MODIS classification algorithm (Sulla-Menashe et al., 2019). For a 500 m pixel with multiple PFTs (Fig. 3a), when one PFT had no MLA measurement, the MLA of the PFT was assigned the value of its nearest neighbor within 100 km with the same PFT. This distance setting (100 km) was based on a previous study which derived global maps of various leaf traits from a limited number of field measurements, remote sensing, and climate data (Moreno-Martínez et al., 2018). In field measurement, the entire canopy MLA is commonly calculated as the average of all measured leaf LIAs weighted by leaf area (see Appendix A) (Zou et al., 2014; de Wit, 1965; Yan et al., 2021). Leaves with larger areas have higher weights. Upscaling MLA from 30 to 500 m follows the same rationale as that from the leaf to canopy scale. For a 30 m pixel with a higher LAI, the weight of the pixel is higher. Therefore, the 500 m MLA was computed as the weighted average of the enhanced vegetation index (EVI2) assuming a linear relationship between LAI and EVI2 (see Appendix A) (Dong et al., 2019; Alexandridis et al., 2019). Although previous studies have reported that the vegetation index may be nonlinearly correlated with LAI because of the saturation effect at medium and high LAI conditions, EVI2 is highly resistant to the saturation effect (Gao et al., 2023). The errors caused by this slight nonlinearity were further analyzed in Sect. 4.4.

$$\begin{array}{}\text{(5)}& {\mathrm{MLA}}_{\mathrm{500}\mathrm{m}}={\displaystyle \frac{\sum {\mathrm{MLA}}_{\mathrm{30}\mathrm{m}}\times \mathrm{EVI}{\mathrm{2}}_{\mathrm{30}\mathrm{m}}}{\sum \mathrm{EVI}{\mathrm{2}}_{\mathrm{30}\mathrm{m}}}}\end{array}$$

The 500 m upscaled MLA samples were further refined to select the most representative samples following three criteria: (1) the coefficient of variation of the 30 m EVI2 in the 500 m pixel is less than 0.2, (2) the vegetation proportion in the 500 m pixel is greater than 0.8, and (3) the proportion of PFTs represented by the MLA measurements in the 500 m pixel is greater than 0.4. The final number of samples after refinement is 3013 with a more even spatial distribution (Fig. 4).

Figure 4 Distribution of global mean leaf inclination angle samples after screening. See Fig. 1 for acronyms.

2.3.2 Global MLA mapping

Different mapping strategies were employed for noncrops and crops (Fig. 3b) considering the small number of valid crop samples (Fig. 4) and the lack of location information for most crop samples. For noncrops, the upscaled 500 m MLA samples were used to train a random forest regressor to predict the global MLA from different features (Table 2). All input features were unified to the 500 m resolution. Therefore, the derived MLA map corresponds to the average MLA at the 500 m scale. Notably, this study used all MODIS BRDF and spectral reflectance data including low-quality data that may be contaminated by clouds. The multiyear aggregation (Table 2) can partly mitigate the influence induced by low-quality observations (Sulla-Menashe et al., 2019). The normalized difference vegetation index (NDVI) was used as the predictive feature because it is strongly coupled with LIA, especially under low- and medium-vegetation-density conditions (Dong et al., 2019; Zou and Mõttus, 2015). To reduce computational complexity and potential overfitting, a feature selection process was conducted based on the variable importance (the sum of the decrease in the Gini impurity index over all trees in the forest) computed by the model, and only the 40 most important variables were used in the final prediction. During the training process, the out-of-bag error was minimized to obtain the optimal hyperparameters. The prediction performance of the random forest regressor was evaluated using a 10-fold cross-validation approach with upscaled MLA samples.

For crops, the measured MLA values were averaged for different crop types as a typical MLA (Table S2). After assigning typical MLAs for different crops with high-resolution crop maps (Table 1), the high-resolution crop MLAs were upscaled to 500 m as training samples (Eq. 5). Only the samples with a crop area ratio > 80 % within a 500 m pixel were selected for training. The crops were further divided into broadleaf crops and cereal crops and processed with the same procedure used for noncrops (Fig. 3b). All procedures were conducted on GEE under the WGS-84 geographic coordinate system.

Two quality layers were added to represent the quality of input data and the prediction model. The input data quality was denoted by the proportion of high-quality BRDF inversions for each pixel. The prediction model quality was represented qualitatively for each pixel considering whether the MLA was predicted by extrapolating beyond the range of the training samples. The random forest model is typically regarded as a black box, and its uncertainty is difficult to quantify in the present study.

2.4 Evaluation of global MLA

The global MLA map was indirectly evaluated using the nadir leaf projection function. The global G(0) was derived from the MLA and evaluated with a high-resolution reference (Sect. 2.2.2) following the upscaling scheme recommended by the Land Product Validation (LPV) Subgroup of the Committee on Earth Observation Satellites (CEOS) (http://lpvs.gsfc.nasa.gov/, last access: 17 April 2023).

Assuming a single-parameter ellipsoidal leaf angle distribution (Campbell, 1990; Wang et al., 2007), the parameter χ, the ratio of the horizontal and vertical axes of an ellipsoid, was first derived from MLA in radians. Compared to other models, the single-parameter ellipsoidal leaf angle distribution is a relatively more accurate and simpler model and has been used in many remote sensing studies (Campbell, 1990; Wang et al., 2007; Kuusk, 2001; Verhoef et al., 2007).

$$\begin{array}{}\text{(6)}& \mathit{\chi }=-\mathrm{3}+{\left({\displaystyle \frac{\mathrm{MLA}}{\mathrm{9.65}}}\right)}^{-\mathrm{0.6061}}\end{array}$$

The G(θ) value in the nadir direction (θ=0°) was calculated using an analytical formula (Leblanc and Fournier, 2017).

$$\begin{array}{}\text{(7)}& G\left(\mathit{\theta }\right)={\displaystyle \frac{\sqrt{\left({\mathit{\chi }}^{\mathrm{2}}+{\tan }^{\mathrm{2}}\mathit{\theta }\right)}\cos \mathit{\theta }}{\mathit{\chi }+\mathrm{1.774}(\mathit{\chi }+\mathrm{1.182}{)}^{-\mathrm{0.73}}}}\end{array}$$

The MLA product was first upscaled to 3 km through a weighted averaging method using the MODIS LAI to derive G(0) (Eq. 7). The reference LAI, FVC, and CI were also upscaled to 3 km through simple averaging to compute the reference G(0) (Eq. 4). The MLA-derived G(0) and the reference G(0) were compared at the 3 km × 3 km area around each site. The correlation coefficient (r), bias, and root mean square error (RMSE) were calculated as the evaluation metrics, as follows:

$$\begin{array}{}\text{(8)}& r=\sqrt{\mathrm{1}-{\displaystyle \frac{\sum _{i=\mathrm{1}}^n{\left({\widehat{y}}_i-y_i\right)}^{\mathrm{2}}}{\sum _{i=\mathrm{1}}^n{\left(y_i-\stackrel{\mathrm{‾}}{y}\right)}^{\mathrm{2}}}}},\end{array}$$

$$\begin{array}{}\text{(9)}& {\displaystyle }& {\displaystyle }\mathrm{Bias}={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n\left(\widehat{y}-y_i\right),\text{(10)}& {\displaystyle }& {\displaystyle }\mathrm{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n{\left(\widehat{y}-y_i\right)}^{\mathrm{2}}},\end{array}$$

where ${\widehat{y}}_i$, y_i, $\stackrel{\mathrm{‾}}{y}$, and n denote the MLA-derived G(0), reference G(0), mean of the reference G(0), and number of G(0) values, respectively.

3 Results

3.1 Global measured LIA values

The species-aggregated LIA was employed in the analysis of global LIA measurements. Figure 5 shows the distributions of global measured LIA values for different PFTs. The global measured MLA is 40.74° and generally follows the order of CRO-C > GRA > ENF > CRO-B > EBF > SHR > DNF > DBF (Table 3). Cereal crops exhibit the highest MLA (59.11°), whereas DBF has the most horizontal leaves (MLA = 34.94°). GRA and EBF show large LIA variations (SD = 20.44 and 17.17°, respectively), whereas CRO-B exhibits a small range. The DNF LIA measurements are only for one species and show very little variation (Fig. 5).

Figure 5 Distribution of global mean LIA (MLA) for different plant functional types (see Fig. 1 for acronyms). The last shape shows the global average. Statistics are conducted for each species as represented by points in the figure.

Table 3 Statistics of leaf inclination angle measured for different plant functional types (PFTs). SD is the standard deviation. The inclination index (χ_L) is converted from mean leaf inclination angle (MLA) (${\mathit{\chi }}_{\mathrm{L}}=\mathrm{2}\cos \left(\mathrm{MLA}\right)-\mathrm{1}$) (Lawrence et al., 2019). The bold format highlights high and low values.

3.2 The relationships between MLA and other variables

Figure 6 shows the importance of the top 40 variables in the MLA prediction obtained from the random forest regression model. The importance of these 40 variables accounts for 78 % of the total importance among all 76 variables. Spectral features account for 30 % of the importance, which is higher than that of other features. Among the spectral features, NDVI, near-infrared (NIR) band, and red band reflectance are most critical for MLA prediction. The importance of BRDF features is comparable to that of climatic variables (21 % vs. 20 %), followed by terrain features (7 %). Among the BRDF features, the NIR BRDF information shows a higher contribution than the red band, with importance in the following order: geometrically scattered kernel > isotropic scattering kernel > volumetric scattering kernel. The importance ranking of the climatic variables follows the order of precipitation ≈ solar radiation > temperature. In addition, elevation, slope, and aspect significantly impact the MLA prediction.

Figure 6 The importance of variables in the mean leaf inclination angle prediction. NIR, red, green, and blue denote the nadir reflectance in the near-infrared, red, green, and blue bands, respectively; geo, iso, and vol represent kernel coefficients of geometric optical surface scattering, isotropic scattering, and volumetric scattering, respectively. The suffixes p (followed by a number), mean, and “std” represent the percent quantile, mean, and standard deviation, respectively.

Figure 7 illustrates the relationships between the upscaled MLA samples and the 16 most important variables. Overall, MLA decreases with the increase in NDVI, NIR reflectance, and NIR BRDF kernel parameters, whereas it increases with the standard deviation of NDVI. MLA is negatively correlated with solar radiation, precipitation, and temperature. Additionally, MLA increases with an increasing standard deviation of solar radiation (corresponding to middle- to high-latitude regions), while it decreases with an increase in the standard deviation of precipitation (corresponding to tropical and subtropical regions with high precipitation). MLA increases slightly with altitude and then decreases.

Figure 7 Relationships between mean leaf inclination angle (MLA) and different predictive variables. See Fig. 6 for different variables.

3.3 Global MLA and G(0) maps

Figure 8 shows the spatial distribution of the global 500 m MLA product. Central Asia (grasslands), southern India (cereal crops), and the central United States (grasslands and cereal crops) show higher MLAs of approximately 60°, whereas the rainforests in South America, central Africa, and Southeast Asia have more horizontal leaves with MLAs of around 30° (Figs. 8 and S2). MLA increases with latitude, from 32.93±7.03° around the Equator (∼1.5° N) to 53.48±3.20° in the northern tundra (∼76.5° N). Variation in MLA decreases as latitude increases (Fig. 8). Among different PFTs, cereal crops show the highest MLA (54.65±6.28°), while evergreen broadleaf forest has the lowest MLA (34.40±6.42°), and PFTs follow the order CRO-C > CRO-B > DNF > SHR > ENF ≈ GRA > DBF > EBF (Table 4). Grassland, broadleaf forest, and evergreen needleleaf forests show larger MLA variations than other PFTs, whereas deciduous needleleaf forests show minimal variation. The global vegetation MLA is 41.47°, with a standard deviation of 9.55°, which is comparable to the MLA of DBF (41.23±6.58°) (Fig. 9a and Table 4).

Figure 8 The global mean leaf inclination angle (MLA) map. The right panel shows the MLA latitudinal mean (solid line) and the standard deviation values (shaded area) weighted by leaf area index.

Table 4 Statistics of global predicted mean leaf inclination angle (MLA), nadir leaf projection function (G(0)), and inclination index (χ_L) for different plant functional types (PFTs). SD is the standard deviation. The χ_L is converted from MLA (${\mathit{\chi }}_{\mathrm{L}}=\mathrm{2}\cos \left(\mathrm{MLA}\right)-\mathrm{1}$) (Lawrence et al., 2019). The bold format highlights high and low values.

Figure 9 Statistics (a) and probability density distributions (b) of the global mean leaf inclination angle (MLA) for different plant functional types. The error bars in (a) represent the standard deviation. The dashed black line and shaded area in (b) indicate the global MLA mean and standard deviation. The dashed gray line represents the MLA (=57.30°) of the spherical leaf angle distribution. The mean, standard deviation, and probability density values are weighted by leaf area index. See Fig. 1 for the acronyms.

The global MLA exhibits an asymmetric probability density distribution toward the lower MLA (Fig. 9b). It presents roughly three peaks, with the highest peak (∼51°) containing DNF, ENF, CRO, GRA, and SHR. The moderate peak (∼35°) is mainly composed of EBF and DBF, while the third peak (∼58°) is dominated by crops. The MLAs of crops and some grasslands are close to the MLA of the spherical distribution (57.30°). The global MLA (41.47°) is 15.83° (38 %) smaller than the MLA of the spherical distribution because the vegetation MLA is mostly less than 57.30° (Fig. 9b).

Figure 10 displays the spatial distribution of global G(0) generated from MLA. Overall, the global G(0) shows a pattern opposite to the global MLA. The G(0) values in central Asia (grasslands, Fig. S3), southern India (cereal crops), and the central United States (grasslands and cereal crops) are relatively lower than those in tropical rainforests and DBFs in the eastern United States. G(0) generally decreases slowly with latitude, from 0.78±0.08 at the Equator (∼1.5° N) to 0.52±0.04 in the northern tundra (∼76.5° N).

Figure 10 The global nadir leaf projection function (G(0)) map. The right panel shows the G(0) mean (solid line) and standard deviation values (shaded area) weighted by leaf area index.

Among different PFTs, EBF has the highest G(0), at approximately 0.76±0.06 (Fig. 11a, Table 4), whereas cereal crops show the lowest value, at approximately 0.52±0.08. The DBF G(0) is comparable to the global average. The G(0) of broad-leaved forests is greater than that of other PFTs (Fig. 11a, Table 4). The global G(0) probability density distribution peaks at 0.52–0.65, with an asymmetric distribution (Fig. 11b). The proportion on the right side of the peak is larger than that on the left. The peak of the global G(0) distribution mainly contains DNF, ENF, CRO, GRA, and SHR. The left side of the peak is mainly composed of crops, while the right side is dominated by broad-leaved forests and some shrubs. The spherical distribution G(0) (0.50) is mainly represented by crops and a small amount of grassland, where G(0) also shows a large variation (∼0.35). The spherical distribution G(0) is 0.18 (26 %) less than the global average G(0) (0.68), as most vegetation G(0) is greater than 0.50 (Fig. 11b).

Figure 11 Statistics (a) and probability density distributions (b) of the global nadir leaf projection function (G(0)) for different plant functional types. The error bars in (a) represent the standard deviation. The dashed black line and shaded area in (b) indicate the globalG(0) mean and standard deviation. The dashed gray line represents the G(0) (=0.50) of spherical leaf angle distribution. The mean, standard deviation, and probability density values are weighted by leaf area index. See Fig. 1 for the acronyms.

Figure 12 demonstrates the global distributions of the MLA quality indicators. The global mean proportion of high-quality BRDF inputs is 68.03 %. The tropical regions have a low proportion of high-quality inputs (20 %) because of cloud contamination (Fig. 12a). Considering the large number of observations for each pixel (7904 from 2001 to 2022), this percentage (20 %) of high-quality observations is sufficient to map MLA. In addition, 80.39 % of the global MLA map was derived within the feature ranges of training samples, and the remaining 19.61 % is mainly located in high-latitude regions and Africa. For the latter areas, the MLA map was predicted with extrapolation, and caution should be taken when using the data in these areas (Fig. 12b).

Figure 12 Global distributions of quality indicators. Panels (a) and (b) show the proportion of high-quality BRDF inversions and whether the prediction is within the ranges of training samples, respectively.

3.4 Evaluation of global MLA

Figure 13 shows the comparison between the predicted MLA and upscaled MLA samples using the 10-fold cross-validation method. For noncrops, the predicted MLA is moderately consistent with the upscaled sample MLA (r=0.75, RMSE = 7.15°), with 83 % of samples having residuals < 10° and 94 % of samples having residuals < 15°. For DNF and SHR, the predicted MLA compresses the variation range of sample MLA (Fig. 13a). For crops, the predicted MLA of CRO-C shows higher consistency (r=0.60) than that of CRO-B (r=0.48) (Fig. 13b and c).

Figure 13 Comparisons between predicted MLA and sample MLA for noncrop (a), broadleaf crops (b), and cereal crops (c) (see Fig. 1 for the acronyms). The error bar in (a) represents the standard deviation.

Figure 14 compares G(0) derived from the MLA and high-resolution reference data. The MLA-derived G(0) shows moderate consistency with the reference G(0) (r=0.62); 65 % of the estimated G(0) residuals are <0.15 and 84 % of the residuals are <0.20. The estimated G(0) generally overestimates (bias = 0.11), especially when G(0) is low (<0.60), mainly for crops, pasture, woody wetlands, and shrubs, whereas grasslands show better consistency. The estimated G(0) is temporally more stable than the reference G(0), which is generally greater than 0.50 and displays seasonal variation (horizontally distributed bars in Fig. 14).

Figure 14 Comparisons of G(0) derived from mean leaf inclination angle and high-resolution reference data for different plant functional types (see Fig. 2 for the acronyms). The error bar represents the standard deviation of reference G(0) in different seasons.

4 Discussion

4.1 Global MLA and G(0)

This study compiled global LIA field measurements and generated the first global 500 m MLA and G(0) maps (Figs. 8 and 10). These maps show the average MLA and G(0) conditions during the growing seasons from 2001 to 2022. Overall, the global MLA is lowest around the Equator and increases with latitude (Figs. 8 and 10). This is in accordance with the MLA latitude variation derived from model simulations (Huemmrich, 2013). Crops have higher MLA than broadleaf forests, the leaves of which are relatively horizontal. The global MLA and G(0) maps enhance our understanding of the global distribution of MLA and G(0) and should be useful in radiative transfer modeling, remote sensing of vegetation parameters, land surface modeling, and ecological studies.

The global MLA shows good consistency with validation samples (Fig. 13) and the statistics of LIA field measurements (Tables 3 and 4), demonstrating its reliability. The globally derived MLA is 41.47°, which is consistent with the LIA measurements (40.74°, Tables 3 and 4). However, the derived MLAs of DBF, DNF, CRO-B, and SHR are approximately 10° greater than the measured MLAs. It is noted that the number and spatial distribution of LIA measurements for these biomes are limited. For example, the global CRO-B areas are dominated by soybeans with higher LIA (Table S2), and the LIA measurements for soybeans are limited, which possibly caused the CRO-B MLA in the global map to be greater than that in the measurement statistics (Tables 3 and 4). The poor crop MLA prediction (Fig. 13b) is mainly caused by a small number of samples and the strong seasonal variation. It is difficult to consider within-crop LIA variation when typical MLA values are assigned to different crops.

Due to the lack of high-resolution reference MLA, the global MLA was evaluated through a comparison of the MLA-derived G(0) with the high-resolution reference G(0) (Fig. 14). This practice was adopted because both MLA and G(0) are closely related. G(0) is typically calculated from the LIA distribution function based on Nilson's algorithm (Nilson, 1971). We calculated G(0) from MLA assuming an ellipsoidal LIA distribution (de Wit, 1965) and found that the calculated G(0) is highly consistent with the reference G(0) calculated from Nilson's algorithm for different theoretical LIA distributions (Fig. S5). The MLA-calculated G(0) also shows a monotonic decreasing relationship with MLA (Fig. S6).

Figure 14 shows medium consistency but MLA-derived G(0) overestimates at low values (<0.60), especially for CRO, PAS, SHR, and WET. The overestimation may be partly caused by the underestimation of MLA at high values that is related to the errors introduced in the sample expansion and upscaling. These errors are mainly caused by a lack of LIA measurements, vegetation structural complexity, and seasonal variation. In addition, the uncertainties in the reference G(0) may have contributed to the overestimation. The reference G(0) was derived from the Beer–Lambert law (Eq. 1), which assumes that the canopy is a turbid medium. The turbid medium assumption is unrealistic for complex vegetation (Widlowski et al., 2014). The angular variation of CI and the mixture of branches and leaves in generating high-resolution G(0) can also lead to the overestimation. Previous studies have shown that CI increases with the view zenith angle (Fang, 2021), which means that the whole CI > CI(0) and can lead to the underestimation of the reference G(0) (Eqs. 3 and  4). Woody materials may introduce biases into the reference G(0) as they were not separated in the high-resolution FVC and LAI products. The mixture of woody materials and leaves may have caused the underestimation of the reference G(0) because trunks usually have higher inclination angles (Liu et al., 2019). The MODIS LAI product used for LIA upscaling in the G(0) validation (Sect. 2.4) is known to have issues such as internal inconsistency, backup algorithm accuracy, and spatiotemporal gaps (Kandasamy et al., 2013; Pu et al., 2023; Zhang et al., 2024). In the future, new improved MODIS LAI can be used in the G(0) validation (Pu et al., 2024; Yan et al., 2024). Compared with the previous G(0) derived from global vegetation biophysical products (Eq. 4) (R²=0.11, RMSE = 0.53) (Li et al., 2022), the MLA-derived G(0) performs better (r=0.62, RMSE = 0.15). In addition, the G(0) data obtained from our study can be used to derive the G(θ) for any arbitrary angle. One method of getting G(θ) is based on a single-parameter ellipsoidal leaf angle distribution (Campbell, 1990) (Eq. 7). Another method is to make use of both G(0) and G (57.3°) (≡0.5) to derive G(θ) using a simple linear ($G\left(\mathit{\theta }\right)=a\cdot \mathit{\theta }+b$) or sinusoidal ($G\left(\mathit{\theta }\right)=a\cdot \sin \left(\mathit{\theta }\right)+b)$) interpolation method. Since G(θ) varies most significantly in the nadir direction for different MLA (Wilson, 1959), the uncertainty of G(θ) derived from the global MLA in other directions is expected be smaller than that of G(0).

4.2 The relationship between MLA and other variables

Analysis of the relationships between MLA and other features in the MLA mapping process reveals that MLA is negatively correlated with NDVI, NIR reflectance, and NIR BRDF kernel coefficients (Fig. 7). These findings are consistent with other simulation and experimental studies (Zou and Mõttus, 2015; Liu et al., 2012; Dong et al., 2019; Jacquemoud et al., 1994). Higher MLA generally means lower canopy interception and higher transmission for high solar altitude, and more soil background can be detected in the nadir direction (Liu et al., 2012). This results in lower (higher) canopy NIR (red) reflectance because of the generally lower (higher) NIR (red) soil reflectance than that of the leaf components (Siegmund and Menz, 2005) and negative correlations between MLA and NIR reflectance and NDVI (Liu et al., 2012). The negative relationships between MLA and radiation, precipitation, and temperature (Fig. 7) are related to the vegetation adaptation mechanism. Under suitable climate conditions (radiation, precipitation, and temperature), horizontal leaves are formed to absorb more radiation and increase the photosynthesis rate (van Zanten et al., 2010; King, 1997). The positive correlation between MLA and the standard deviation of radiation and temperature (Fig. 7) indicates that the MLA is more vertical in areas with significant seasonal changes in radiation and temperature (mid to high-latitude areas) because vertical leaves maximize intercepted radiation under low solar altitudes at middle- to high-latitude areas (Huemmrich, 2013).

Plant functional types were initially used as a predictive variable (Tables 1 and 2), but relatively low importance was found for LIA prediction (Fig. 6, ranked 47 out of 76). This may be because the biome information is implicitly included in the spectral features as the former is frequently derived from the latter (Sulla-Menashe et al., 2019). Previous studies have demonstrated that the LIA variation within PFTs may be larger than that between PFTs. This indicates that the PFT is not a good predictor (Prentice et al., 2024). To avoid overfitting, only the most important 40 features were used for MLA prediction (Fig. 6). To explore the regional differences of the variable importance, an analysis was conducted for the tropical (23.5° S–23.5° N), northern temperate (23.5–60° N), northern polar (60–90° N), and southern temperate (23.5–60° S) zones. The 40 most important variables are similar among different regions although minor differences exist (Fig. S7). Among the 40 variables for tropical, northern temperate, northern polar, and southern temperate zones, 32, 35, 30, and 31 of them, respectively, are the same as the 40 global variables (Fig. S7). Climate and spectral variables are significant among all regions, whereas BRDF features are the most important in the southern temperate zone. The 40 most important variables in the global MLA prediction account for ∼80 % of total importance among different regions, which is similar to that in the global prediction.

4.3 Use of the new MLA map

The spherical LAD assumption has been widely adopted in the literature (Tang et al., 2016; Zhao et al., 2020; Wang and Fang, 2020). This study demonstrates that the spherical assumption is valid only for cereal crops, but not for broadleaf forests (Tables 3 and 4). This finding is consistent with previous local LIA measurements (de Wit, 1965; Pisek et al., 2013; Yan et al., 2021). For crops, the spherical assumption may become invalid because of seasonal change and species diversity (Table S2, Figs. 5 and 9). In addition, most (72 %) of the reference G(0) values are greater than 0.50 (Fig. S8); in this case, the spherical distribution would underestimate the radiation and rainfall interception because of the overestimated LIA and underestimated G(0) for most conditions (Figs. 9 and 11) (Stadt and Lieffers, 2000). In current LSMs, a constant LIA is commonly assigned for each PFT (Majasalmi and Bright, 2019). For example, the Community Land Model V5 (CLM5) (Table S4) (Lawrence et al., 2019) uses lower inclination indices and higher LIA values than our results (Tables 3 and 4) and may thus underestimate canopy interception. The global LIA map generated in this study provides a more reasonable LIA parameterization strategy for the application communities.

4.4 Limitations and prospects

This study is mainly limited by the small number of LIA measurements, especially continuous measurements. First, within-species LIA variations were neglected in the spatial expansion due to limited spatial coverage of existing LIA-measured data (Sect. 2.3.1). This may introduce some errors, especially for crops. Second, three different sources of LIA measurements were gathered with different measurement schemes, and uncertainty may arise because of these differences. The random forest algorithm is robust to these differences because some of the samples and features were randomly selected and the algorithm ensembled the predictions from multiple decision trees (Svetnik et al., 2003). We manually inspected all field LIA data and are confident in their quality. Third, for forests, the contribution of the understory was not considered. Typically, the understory is characterized by more horizontal leaves, and ignoring the understory may lead to an MLA overestimation (Utsugi et al., 2006). Nevertheless, a previous study showed that the relative contribution of the understory to the overall MLA is less than 10 % (Li et al., 2022). Finally, only the growing season MLA was calculated, whereas the seasonal and long-term variations of MLA were not considered due to the lack of continuous LIA measurements.

We assumed a linear LAI–EVI2 relationship ($\mathrm{LAI}=a\cdot \mathrm{EVI}\mathrm{2}$) to upscale MLA from the canopy to 500 m scale (Sect. 2.3.1 and Appendix A). Global analysis of MODIS LAI and EVI2 products shows a slight nonlinear relationship between them (Fig. S9). The nonlinear relationship was also used to upscale MLA (Eq. A2) in a side experiment, where the derived MLA was found to be consistent with the original one (Fig. S10) because of the homogeneity of the 500 m pixel after rigorous sample screening (Sect. 2.3.1). This demonstrates the suitability of the linear assumption.

In the future, more efficient LIA observation systems should be developed to provide continuous LIA data (Kattenborn et al., 2022). LIA measurements can be integrated into existing ground observation networks, such as the National Ecological Observatory Network (NEON) (Kao et al., 2012), Integrated Carbon Observation System (ICOS) (Gielen et al., 2018), and Terrestrial Ecosystem Research Network (TERN) (Karan et al., 2016), to enhance temporal LIA measurements in a larger spatial extent, especially for DNF and crops. Using standard LIA measurement protocols will certainly improve the LIA data consistency (Li et al., 2023). In addition, canopy structure parameters are interrelated, and introducing other structure parameter products, such as LAI, FVC, CI, and canopy height, as predictive variables may improve the MLA prediction. Multiangle reflectance (Jacquemoud et al., 2009; Goel and Thompson, 1984; Jacquemoud et al., 1994) and light detection and ranging (Zheng and Moskal, 2012; Bailey and Mahaffee, 2017; Itakura and Hosoi, 2019) are encouraging remote sensing tools that can help to derive temporally continuous and high-resolution MLA data.

5 Data availability

The global MLA and G(0) products (CAS-GLA V1.1) are available at https://doi.org/10.5281/zenodo.12739662 (Li and Fang, 2025). The related GEE code can be accessed for GEE user at https://code.earthengine.google.com/?accept_repo=users/SiJia/MTA (created by Sijia Li, last access: 12 January 2025).

6 Conclusion

This study compiled existing global LIA measurements and generated the first global 500 m MLA and G(0) products by gap-filling the LIA measurement data using a random forest regressor. The mean of global LIA measurements is 40.74° and cereal crops show the highest MLA (59.11°). The global estimated MLA shows an explicit spatial distribution and the value increases with latitude. The global MLA is 41.47°±9.55° and follows the order of CRO-C > CRO-B > DNF > SHR  > ENF ≈ GRA > DBF > EBF. The predicted MLA presents a medium consistency (r=0.75, RMSE = 7.15°) with the validation samples for noncrops. For crops, the results are relatively poorer (r=0.48 and 0.60 for broadleaf crops and cereal crops) because of limited LIA measurements and strong seasonality. The estimated G(0) is moderately consistent with the reference G(0) (r=0.62).

The MLA and G(0) products obtained in this study could enhance our understanding of global LIA and assist remote sensing retrieval and land surface modeling studies. These products provide a more realistic parameterization strategy than the commonly used spherical LAD and PFT-specific MLA assignment. Note that the global MLA and G(0) products mainly represent the typical state during the growing season. These products can be further improved and temporal MLA data can be obtained through continuous measurement and remote sensing retrieval.

Appendix A: Upscaling LIA from leaf to canopy and ecosystem scales

From leaf to canopy scale, the entire canopy MLA is commonly calculated as the average of all measured leaf LIAs weighted by leaf area (Eq. A1) (Zou et al., 2014; de Wit, 1965; Yan et al., 2021). In practice, because of the difficulty of leaf area measurement, especially for a large number of leaves, the variability of leaf areas within a canopy is often ignored and the areas of all leaves are assumed to be similar. In this case, the canopy LIA can be simplified as the average LIA weighted by leaf number (Eq. A1) (Ryu et al., 2010; Pisek et al., 2011; Chianucci et al., 2018):

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}{\mathrm{MLA}}_{\mathrm{canopy}}& ={\displaystyle \frac{\sum _i{\mathrm{LIA}}_i\cdot {\mathrm{LA}}_i}{\sum _i{\mathrm{LA}}_i}}={\displaystyle \frac{{\mathrm{LA}}_{\mathrm{mean}}\cdot \sum _i{\mathrm{LIA}}_i}{{\mathrm{LA}}_{\mathrm{mean}}\cdot N}}\\ & ={\displaystyle \frac{\sum _i{\mathrm{LIA}}_i}{N}},\end{array}\end{array}$$

where MLA_canopy is the MLA at canopy scale, i is the ith leaf, LIA is leaf inclination angle, LA is single leaf area, LA_mean is the mean leaf area by ignoring the variation of leaf area within a canopy, and N is the number of leaves within a canopy.

From the canopy to 30 m scale, the canopy-level MLA is regarded as equal to 30 m MLA because for MLA measurements, the dominant species was artificially identified by investigators and the spatial representativeness at the extent of 30 m is ensured.

From 30 to 500 m, the 500 m MLA was formulated as the weighted average of 30 m MLA by the leaf area of the 30 m pixel (Eq. A2), the same as that from the leaf to canopy scale. The leaf area of a 30 m pixel can be deduced from the product of leaf area index (LAI) and the ground area of a 30 m pixel according to the definition of LAI (half of the green leaf area per unit of ground area) (Eq. A2) (Fang et al., 2019):

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}{\mathrm{MLA}}_{\mathrm{500}}& ={\displaystyle \frac{\sum _j{\mathrm{MLA}}_{\mathrm{30}\mathit{\_}j}\cdot {\mathrm{LA}}_{\mathrm{30}\mathit{\_}j}}{\sum _j{\mathrm{LA}}_{\mathrm{30}\mathit{\_}j}}}\\ & ={\displaystyle \frac{\sum _j{\mathrm{MLA}}_{\mathrm{30}\mathit{\_}j}\cdot {\mathrm{LAI}}_{\mathrm{30}\mathit{\_}j}\cdot S}{\sum _j{\mathrm{LAI}}_{\mathrm{30}\mathit{\_}j}\cdot S}}\\ & ={\displaystyle \frac{\sum _j{\mathrm{MLA}}_{\mathrm{30}\mathit{\_}j}\cdot {\mathrm{LAI}}_{\mathrm{30}\mathit{\_}j}}{\sum _j{\mathrm{LAI}}_{\mathrm{30}\mathit{\_}j}}},\end{array}\end{array}$$

where MLA₅₀₀ and MLA₃₀ represent MLA at 500 and 30 m scales, j is the jth 30 m pixel, LA_{30_j} is the total leaf area of a 30 m pixel, LAI_{30_j} is the leaf area index (m² m⁻²) of a 30 m pixel, and S is the ground area of a 30 m pixel.

Assuming $\mathrm{LAI}=a\cdot \mathrm{EVI}\mathrm{2}+b$ and b≈0 (as illustrated in Fig. S9), the MLA at 500 m scale can be calculated as

$$\begin{array}{}\text{(A3)}& {\mathrm{MLA}}_{\mathrm{500}}={\displaystyle \frac{\sum _j{\mathrm{MLA}}_{\mathrm{30}\mathit{\_}j}\cdot \mathrm{EVI}{\mathrm{2}}_{\mathrm{30}\mathit{\_}j}}{\sum _j\mathrm{EVI}{\mathrm{2}}_{\mathrm{30}\mathit{\_}j}}}.\end{array}$$