## Materials and Methods
### Mapping Heterogeneity of Canopy Properties  
To evaluate the upscaling of field LAI measurements to satellite products, Landsat Level-2 Surface Reflectance products including Enhanced Vegetation Index (EVI) and Normalized Difference Vegetation Index (NDVI) derived from Landsat 5 (TM), Landsat 7 (ETM+) [63] and Landsat 8 (OLI) [64] scenes are also analyzed. EVI is calculated from red (R), near-infrared (NIR) and blue band (B) based on Equation (1) and incorporates parameters to adjust for canopy background, atmospheric resistance.  
\[EVI\ =\ 2.5*((NIR-R)/(NIR+6*R-7*B+1)) \tag{1}\]  
NDVI is calculated as a ratio between the R and NIR values in a traditional fashion (Equation (2)).  
\[NDVI\ =\ (NIR-R)/(NIR+R) \tag{2}\]  
Landsat 5 (TM), Landsat 7 (ETM+), and Landsat 8 (OLI) Level-2 surface reflectance products were generated, using the Land Surface Reflectance Code (LaSRC) and the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) algorithms [65; 66]. The LEDAPS and LaSRC surface reflectance algorithms correct for the temporally, spatially, and spectrally varying scattering and absorbing effects of atmospheric gases and aerosols, which is essential to derive the Earth's reflectance surface values. Landsat Spectral Indices are generated at 30 m spatial resolution on a Universal Transverse Mercator (UTM) mapping grid. The Spectral Indices could be obtained from the USGS Earth Resources Observation and Science (EROS) on-demand processing system.  
The GTOPO 30 1 km elevation map [67] in the China region with Landsat scene footprints at each site is presented in Figure S2. In each of the six sites, there are one to four sampling subplots. The location of each site with the sampling subplots is shown on the Landsat scenes with path and row information (Figure S3) and on the 90 m SARTM elevation map [68], respectively (Figure S4). There are in total 778 Landsat 5-7 and Landsat 8 scenes with good quality (collection category is Tier 1) and cloud cover smaller than 20% of the total scene area from the year of 2005 to 2017 for all of these forests. Because there is no quality assessment information for the Landsat spectral indices, we use the pixel quality assurance band (pixel_qa) instead to obtain quality information. Pixels with cloud shadow, aerosols, etc., are filtered out, and we keep pixels with the best quality.  
Landsat EVI and NDVI were extracted at the field sampling plots spatial level (20\(\sim\)100 m). Based on the Landsat scene quality information, we cleaned and kept just good quality Landsat pixel data. Then, daily EVI and NDVI were retrieved from the area-weighted averaged Landsat pixels inside the field sampling plot. We averaged the daily EVI and NDVI values by month. Finally, the monthly EVI and NDVI at the field sampling plots level were regressed with corresponding monthly in situ LAI values to calculate the statistics (R\({}^{2}\), etc.).  
We investigated whether forest cover fraction in different LAI products' spatial resolution levels influenced product bias and error. We use the 30 m resolution GlobeLand30 land cover product for 2010 [69] to extract the forest cover fraction inside the relevant pixels for each of the LAI products at each site. We classify the land cover for each 30 m pixel on different spatial scales (LAI products' scale: 300 m, 500 m, and 1 km) at each site and plot. We only kept the 30 m pixels located entirely inside the LAI pixels. The total number of 30 m pixels and pixels classified as forests were both counted and summed to determine the total area and the forest areas at different spatial scales. Then we can obtain the forests cover proportion at different spatial scales for each plot and site. Finally, the forests' cover proportion was compared with the LAI products' bias at each site.