A unified vegetation index for quantifying the terrestrial biosphere

See allHide authors and affiliations

Science Advances  26 Feb 2021:
Vol. 7, no. 9, eabc7447
DOI: 10.1126/sciadv.abc7447


Empirical vegetation indices derived from spectral reflectance data are widely used in remote sensing of the biosphere, as they represent robust proxies for canopy structure, leaf pigment content, and, subsequently, plant photosynthetic potential. Here, we generalize the broad family of commonly used vegetation indices by exploiting all higher-order relations between the spectral channels involved. This results in a higher sensitivity to vegetation biophysical and physiological parameters. The presented nonlinear generalization of the celebrated normalized difference vegetation index (NDVI) consistently improves accuracy in monitoring key parameters, such as leaf area index, gross primary productivity, and sun-induced chlorophyll fluorescence. Results suggest that the statistical approach maximally exploits the spectral information and addresses long-standing problems in satellite Earth Observation of the terrestrial biosphere. The nonlinear NDVI will allow more accurate measures of terrestrial carbon source/sink dynamics and potentials for stabilizing atmospheric CO2 and mitigating global climate change.


Quantifying vegetation cover, biochemistry, structure, and functioning from space is key to study and understand global change, biodiversity, and agriculture. In practice, remote sensing has relied vastly on the use (and abuse) of vegetation indices (VIs) derived from spectral reflectance owing to their generally decent performance. VIs are parametric transformations of a few spectral bands designed to maximizing their sensitivity to particular biophysical phenomena (e.g., greenness, water content, or photosynthetic activity) while minimizing their sensitivity to factors such as soil properties, solar illumination, atmospheric conditions, and sensor viewing geometry. A plethora of narrow-band indices has been proposed in the literature (1). Indices are designed for specific applications and conditions, and their parameters are fixed empirically.

The most widely used VI in Earth observation is undoubtedly the normalized difference vegetation index (NDVI) (2, 3). This index exploits the fact that green healthy vegetation shows contrasting behavior in how it reflects red and near-infrared (NIR) radiation. The more chlorophyll there is in a canopy, the more visible light (including the red) can potentially be absorbed to drive photosynthesis, and thus the higher the absorbed energy that can potentially be consumed in carbon fixation. On the other hand, as more living plant biomass is present, the vegetation will scatter and reflect more NIR radiation, which is unusable for photosynthesis. By calculating the difference between bands measuring red and NIR reflectances, NDVI accentuates the particular signature of green vegetation while attenuating undesired influences from nonvegetative elements. NDVI, and other similar indices, have proven effective in assessing chlorophyll content (4, 5), being a good proxy of vegetation density parameters, like the leaf area index (LAI) and the fractional vegetation cover (FVC) (68), as well as the fraction of absorbed photosynthetically active radiation (fAPAR). The success of NDVI relies on its ease of use and its availability over long observational records expanding more than three decades, notably thanks to the Advanced Very High Resolution Radiometer (AVHRR), Landsat optical sensors (Multi Spectral Scanner, Thematic Mapper, Enhanced Thematic Mapper, Operational Land Imager), and the Moderate Resolution Imaging Spectroradiometer (MODIS).

However, NDVI has two major limitations. First, the relationship between NDVI and green biomass is nonlinear and saturates. Some indices such as the enhanced vegetation index (EVI) (9) have tried to compensate for this using information from other bands, but the saturation problem remains. Other approaches have tried to improve NDVI heuristically to obtain a good proxy of both fAPAR and light-use efficiency, and hence suggested it for gross primary productivity (GPP) estimation (10). Actually, some authors have proposed NDVI2 (11) and other arbitrary exponentiations (12) to cope with the nonlinear issue. The second issue is that VIs, by construction, react to the presence of green leaves, but not to photosynthesis per se. GPP can thus decline without any leaf abscission (i.e., a reduction of LAI) or reduction in chlorophyll. A relatively new way to estimate GPP variability from satellite measurements to retrieve sun-induced chlorophyll fluorescence (SIF) (13). However, the relationship between canopy GPP and SIF retrieved from space is still not fully understood (14), and more importantly, this technique is still only available with an overly coarse spatial resolution and a very shallow temporal archive (15, 16).

Using radiative transfer models, Sellers et al. (1719) noted early on that NIR reflectance is a better proxy for fAPAR than NDVI. The problem is then to disentangle the fraction of the NIR that is reflected from the vegetation from the remaining fraction of NIR reflected from nonvegetated elements within a mixed pixel. To address this issue, Badgley et al. (20) proposed considering NDVI as a proxy for vegetation coverage instead of a proxy for fAPAR, and thus multiply NDVI times NIR to calculate a new index, NIRv, which shows high correlations with SIF and GPP at specific temporal scales. Despite its wide reception in the community, NIRv also raises some intriguing questions. For example, given that fAPAR is estimated by both of its components (NIR and NDVI), how does this affect the interpretation of the index? Also, as NIRv linearly scales with the NIR reflectance, how does it deal with saturation? Last, NIRv still uses the same bands as NDVI, but it is not clear how the adopted approximations and assumptions affect NIRv or whether it exploits all available information in these spectral bands.

This paper introduces a methodology to generalize the broad family of VIs based on differences and ratios of spectral bands. Unlike previous approaches to improve indices based on principled (10, 20, 21) or heuristic parametric transformations (11, 12, 22, 23), here we adopt a machine learning standpoint using the theory of kernel methods, which has been widely used to derive nonlinear algorithms from linear ones, while still resorting to linear algebra operations (24, 25). Kernel methods map the involved spectral bands using a nonlinear feature map to a higher dimensional space where the index is defined. The calculation can be expressed in terms of the spectral channels by the definition of a kernel (similarity) function, so one does not need to define the feature map explicitly. The main property of kernel methods is that of linearizing the problem, which is what most of the indices seek either heuristically or based on first principles. Also, by using a particular kernel function, we have guarantees that all higher-order relations between the spectral channels are accounted for, not just the first-order ones. For example, when using differences between NIR and the red bands, the kernel function summarizes all monomials of the differences too, i.e., {NIR-red, (NIR-red)2, (NIR-red)3, …} in a single scalar. Although kernel methods can, in principle, be applied to any VI (see section S1.5 and table S1), the framework is illustrated here to generalize NDVI, largely because of the long history and wide utility of this index, most notably to perform global and long-term studies. We specifically define the NDVI in Hilbert spaces and adopt the radial basis function (RBF) reproducing kernel, k(NIR,red)=exp((NIRred)2/(2σ2)), where the σ parameter controls the notion of distance between the NIR and red bands. The presented kernel NDVI (kNDVI) reduces to computekNDVI=tanh ((NIRred2σ)2)where σ is a length-scale parameter to be specified in each particular application and represents the sensitivity of the index to sparsely/densely vegetated regions. A reasonable choice is taking the average value σ = 0.5(NIR + red) (see sections S1 and S2 for mathematical and ecophysiological justifications), which leads to a simplified operational index version expressed as kNDVI=tanh (NDVI2). The selection of the kernel function and prescription of its parameter allows the kNDVI to perform an automatic and pixel-wise adaptive stretching and guarantees that all moments of the relations between the NIR and red channels are taken into account. This also allows kNDVI to cope with saturation effects, complex phenological cycles, and seasonal variations, to deal with the mixed-pixel problem (20), and to propagate lower uncertainty than other indices (section S2.5). It can be shown that kNDVI actually generalizes NDVI and NIRv theoretically (see sections S1 and S2 and Properties S2.1 and S2.2), which ensures an improved performance. Last, the presented methodology, and the kNDVI in particular, are easy to implement and use in practice (section S10), which is of paramount relevance in operational studies.


We show that kNDVI exhibits consistently stronger correlations than NDVI and NIRv in key independent products [GPP at flux tower estimates and SIF from Global Ozone Monitoring Experiment–2 (GOME-2)]. In general, the proposed index performs better than NDVI and NIRv in all applications, biomes, and climatic zones. The kNDVI is more resistant to saturation, bias, and complex phenological cycles and shows enhanced robustness to noise and stability across spatial and temporal scales (sections S6.2 and S6.3). Additional results for approximating MODIS LAI (section S4), correlation with other related parameters (like fAPAR and FVC) acquired in situ (section S7), crop yield estimation (section S8), and kNDVI’s use for image change detection (section S9) further confirm the validity of the approach. All these properties and performance are achieved without adopting any specific assumption, just exploiting all higher order statistical relations between the involved reflectances.

Accurate proxy to GPP

We evaluated and compared the performance of kNDVI with NDVI and NIRv as a GPP proxy using flux tower GPP estimates from the FLUXNET database (section S5). The proposed kNDVI provides correlations with GPP similar to or better than the other indices over all considered biomes and across all the 169 flux tower sites (Table 1). The weakest relationships are observed for evergreen broad-leaved forests, which can be expected because of the stronger saturation effect in such ecosystem (similarly clear when using the index for LAI estimation, see section S4). The kNDVI excels in each biome individually, confirming its adaptive nature, and globally shows a clear gain (Fig. 1). Although photosynthesis is driven by the amount of vegetation photosynthetic mass within a pixel, solar irradiation and environmental constraints also play a critical role. The latter is not accounted for by the spectral information provided by NIR and red bands. This explains why all indices present lower correlation with GPP and SIF than with LAI for all biomes (Table 1, cf. section S4). The correlation is, however, higher for the kNDVI in almost all cases. Alternative measures of nonlinear association between GPP and the indices, such as Spearman’s correlation (26), mutual information (27), and distance correlation (28), yielded similar results and conclusions (see sections S5 and S6), thus confirming the good capabilities of kNDVI to implicitly linearize the problem.

Table 1 Temporal correlation coefficient between the VIs and the parameters GPP and SIF per biome.

Only vegetation biomes are considered and classes in IGBP were grouped as indicated in parentheses: C1 = NF=Needle-leaf Forest (1 + 3), C2 = EBF = Evergreen Broadleaf Forest (2), C3 = DBF=Decidious Broadleaf Forest (4), C4 = MF = Mixed forest (5), C5 = SH=Shrublands (6 + 7), C6 = SAV=Savannas (8 + 9), C7 = GRA = Herbaceous (10), C8 = CRO=Cultivated (12). Best results per biome indicated in bold and darker green indicates higher correlation.

Embedded Image
View this table:
Fig. 1 Correlation between the indices and parameters.

Histogram of the correlation coefficient between the VIs and the parameters: for GPP (left) correlation computed over 169 FLUXNET sites, and for SIF (right) averaged over all 506 global images.

We studied the robustness of the indices across sites. Figure 2 shows the density and boxplots of the slopes (scaled between 0 and 1) for all 169 flux tower sites. The NIRv index shows a mean closer to 0.5, but the spread is higher than for the kNDVI. Both NDVI and NIRv show very wide whiskers (and hence pathological behaviors and high sensitivity to outliers), while kNDVI shows higher robustness and stability across sites. A simple analysis over all the towers shows that kNDVI outperformed in 84 of the towers (50%), NIRv in 59 (35%), and NDVI in 26 (15%). The kNDVI gains are more noticeable in deciduous and evergreen forests, which confirms the good adaptation to varying photosynthetic phenology of different biomes, primarily forests. This is confirmed when looking at the seasonal patterns of stand photosynthesis for some illustrative sites in Fig. 3, expressed as monthly GPP. For example, the CA-TP4 (Ontario–Turkey Point 1939 Plantation White Pine site) is a region dominated by densely covered woody vegetation and displays green foliage all year round. Unlike NDVI that shows relatively too much and too little sensitivity, respectively, to seasonally changing GPP, the kNDVI follows much better the temporal shape and captures the higher and lower GPP values too. This might be due to the subtle pigment shifts that are largely invisible to NDVI, but may be more detectable by kNDVI, as it was recently shown with NIRv (29). For grasslands, like the CH-Oe1 (Oensingen, Switzerland), neither NDVI nor NIRv can disentangle the phenological cycle of the vegetation from the background noise, while the kNDVI returns acceptable results with larger dynamic range. Here, the tree and shrub cover is less than 10% and a permanent mixture of water and herbaceous or woody vegetation is observed, inducing a strong mixed-pixel problem aggravated by complex topography. The IT-Ro1 (Roccarespampani-1 near Viterbo site) is a deciduous broad-leaved forest consisting of broadleaf tree communities with a clear annual cycle of long leaf-on and leaf-off periods, which are followed faithfully by the kNDVI index. NIRv and kNDVI reveal very similar characteristics. An interesting case is that of closed shrublands. The mixed shrub foliage in the Kennedy Space Center site CSH US-KS2, which can be either evergreen or deciduous, is efficiently handled by kNDVI (R = 0.72) over NIRv (R = 0.68) and NDVI (R = 0.57). Here, unlike NIRv, the proposed kNDVI does not over- and underestimate GPP. Overall, we observed that the kNDVI closely tracked the seasonal dynamics of photosynthesis, presenting a better agreement with GPP. This is achieved by adaptively stretching the dynamic range to better capture time-series extremes (e.g., sparsely and densely vegetated, as well as cold and dry regions). The proposed kNDVI seems to largely correct for “background effects” (important in sparse vegetation or snow) and saturation and may be more sensitive to subtle greenness shifts (e.g., evergreens) based on pigments rather than structure per se.

Fig. 2 Goodness of fit between the indices and GPP.

Distribution of slopes of site-level linear regressions (normalized between 0 and 1) between the indices and biweekly GPP from 169 FLUXNET sites.

Fig. 3 Monitoring GPP at tower level.

Illustrative results over four flux towers covering evergreen needle-leaved forests (CA-TP4), grasslands (CH-Oe1), deciduous broadleaf forest (IT-Ro1), and closed shrublands (US-KS2).

Closer monitoring of photosynthetic activity of ecosystems

Recent studies have linked SIF and VIs, such as NDVI and NIRv (20), as a pragmatic alternative to more sophisticated machine learning approaches (30). We here evaluate the kNDVI computed from the MODIS reflectance bands to approximate globally gridded GOME-2 SIF at 16-day temporal resolution. Despite the fact that GOME-2 can measure both SIF and the NIR and red bands simultaneously, we intentionally estimated all indices independently from coincident MODIS data (see processing details in section S6). We computed the correlation between time series. The kNDVI outperforms the other indices in general (Fig. 1) and in all biomes individually (Table 1), especially in DBF, GRA, and CRO: 5 to 11% gain in correlation over NIRv and 20 to 35% over NDVI.

To confirm the robustness to capture extreme SIF values, we studied the spatial maps of temporal correlation coefficients. kNDVI dominates in all regions (Fig. 4, top) and correlates better with SIF in 69.69% of the pixels compared to NIRv and in 91.32% of cases compared to NDVI (see Fig. 4, bottom). Results suggest that the kNDVI clearly outperforms the other indices in densely vegetated tropical (e.g., Amazonia and Indonesia) and arid regions (e.g., Australia and Mediterranean). As for the case of GPP, other measures of correlation yielded identical conclusions (table S7). Further analysis confirmed the dominant performance of kNDVI in all latitudes, especially in higher and lower ones (table S10), as well as in all climatic zones, especially in the arid and cold regions (table S9).

Fig. 4 Temporal correlation between indices and SIF globally.

Top: Color composite of indices-to-SIF correlation, (R, G, B) = (NIRv, NDVI, kNDVI). Bluish means kNDVI outperforms the rest, which generally happens [in 91.32% of the pixels over NDVI (left) and 69.69% of the cases over NIRv (right)] and particularly in the extreme (low and high) vegetation covers or in cold and dry regions. Bottom: Differences of correlation-with-SIF between the proposed index kNDVI and NDVI (left) and NIRv (right), both globally and for extreme regions. Red colors indicate a higher correlation for kNDVI, and blue indicates a lower correlation for kNDVI (relative to the other indices).

The study areas in Fig. 4 showed the biggest differences between the kNDVI and NDVI and NIRv, and are further scrutinized in Fig. 5. The kNDVI provides improved fit scores in all cases, larger excursions in general, and more resistance to noise and saturation. The higher accuracy by kNDVI (e.g., in California, +19% in R over NIRv) comes mainly from the better behavior in the presence of sharp phenological cycles. In the Iberian peninsula, kNDVI and NIRv perform similarly in quantitative terms, but the proposed kNDVI appears less affected by high-frequency components and covers the whole dynamic range nicely. In Australia, the favorable numerical gain in R (+25%) and the much lower scatter highlight that kNDVI better approximates SIF and closely follows the cycles (especially in March-April-May periods). Despite the big challenges in the Amazon for SIF estimation with GOME-2, the kNDVI can be a more convenient choice compared to other indices, as it deals better with noise and background effects (e.g., soil, standing water, or snow). All in all, the proposed kNDVI seems better qualified to cope with noise, saturation, and complex phenologies.

Fig. 5 Temporal analysis over selected study areas.

Scatterplots of the different indices versus SIF (left), and the average time series over the study areas (right). Axes limits were optimized to improve visualization of all indices.

Similar conclusions were obtained when we studied spatial correlations through time: The proposed index achieves noticeable improvements over NDVI and NIRv, especially between August and November, thus improving autumn phenology owing to its adaptive stretching (see fig. S10 in section S6). The kNDVI is more competitive at finer temporal resolutions (native biweekly) with a noticeable advantage over NDVI (+15%) and NIRv (+4%), but the gain over NIRv disappears at bimonthly scales, since the temporal aggregation induces a “more linear” problem. Likewise, a broader spatial aggregation (from 0.5 up to 2) yielded improved results of all indices, but kNDVI still outperformed the others independently of the spatial scale (section S6 and fig. S11).

We lastly studied the capabilities of kNDVI to deal with the mixed-pixel problem (Fig. 6). Both kNDVI and NIRv scale with the total NIR, NIRT, unlike NDVI that clearly saturates. The kNDVI strongly correlates with SIF over highly vegetated pixels, but the correlation decreases with lower vegetated fractions (Fig. 6). The difference between kNDVI and NDVI stands out, and kNDVI is slightly higher correlated with SIF than NIRv, thus suggesting that the index can reliably isolate the proportion of reflectance attributable to vegetation as well. These properties emerge directly from the NIR-red relations since no assumption is made in designing the index. Accounting for all NIR-red relations allows us to optimally disentangle the mixed-pixel problem efficiently, especially in the densely vegetated areas (e.g., LAI and GPP phenology of crops in section S4 and Fig. 3).

Fig. 6 Correlation with vegetated fraction.

Correlation coefficient between the indices and SIF increases with vegetated fraction (computed from NDVI percentiles). We include the total NIR, NIRT, as a reference. The lower bounds of the NDVI quartiles are as follows: 0, 0.25, 0.50, and 0.75.

The study of natural and agricultural systems should greatly benefit from the kNDVI proposed here because of its solid theoretical foundation combined with its ease of calculation and application. The high correlation with GPP and SIF across all biomes, especially in grasslands, croplands, and mixed forests as well as in arid regions, suggests that the index can efficiently cope with both the saturation and the mixed-pixel problems encountered with traditional indices. The proposed kNDVI explains a large fraction of the variance of GPP at flux tower level, showed good robustness capabilities to noise and saturation, and enhanced stability across space. The kNDVI also highly correlates with SIF derived from an independent sensor, paving the way toward improving our quantification and understanding of photosynthesis at the global scale. Its application and usefulness goes beyond vegetation monitoring and embraces change and extreme detection, phenological and greening studies, upscaling parameters, and all applications where VIs in general and NDVI in particular have previously demonstrated their utility. Our results demonstrate that an agnostic statistical approach is sufficient to explain most of the observed signal.

The kernel methods framework allowed us to generalize all VIs, but we focused on the NDVI case only. Kernel methods, in general, and the kNDVI, in particular, implement the original operation (e.g., NDVI) in a high-dimensional feature space where spectral bands are mapped to. The solution of kNDVI is thus a nonlinear version of NDVI. The framework allows us to accomplish the ever-sought linearization operation implicitly. This means that no ad hoc parametric transformations are needed, just the kernel operation. This also implies that virtually no gain should be obtained over other indices when the relation between the bands and the parameter of interest is linear, such as for instance when an appropriate PAR normalization is applied (see sections S5.2 and S6.4) or whenever one averages over larger spatial or temporal scales (see section S6.3). Our results, however, suggested that the kNDVI instantiation improved results in all problems, even when the domain was previously linearized. This makes the index a very powerful and practical default choice. We anticipate a wide use and development of the proposed index in particular, and of the family of nonlinear VIs in general, to derive informative indicators for operational Earth monitoring and the quantification of the terrestrial biosphere vital signs.


Datasets and processing

GPP and FLUXNET data. The GPP data were obtained from FLUXNET, which is a collection of sites from multiple regional networks (31). This network provides a compilation of in situ observations to measure the exchanges of carbon dioxide, water vapor, and energy between the biosphere and atmosphere (32). To calculate the GPP, the carbon dioxide flux, i.e., net ecosystem exchange, is measured by means of the eddy covariance method. This flux is further partitioned into ecosystem respiration and GPP [gC m−2 day−1] using the daytime (33) or nighttime (34) partitioning methods. For our analyses, we used GPP estimates from the freely available Tier 1 dataset that were obtained with the daytime partitioning method. Of all available sites (212), we selected a subset of 169 sites corresponding to natural vegetation having less than 50% of missing remotely sensed data due to cloud contamination. In addition, we only considered sites where we had more than 4 months of available flux data.

SIF data from GOME-2. We generated GOME-2 0.5 fluorescence at 740 nm and reflectance at 670 and 780 nm from level 2 data obtained from measurements of the GOME-2 sensor flying onboard MetOp-A. The retrieval algorithm of SIF [mW/m2/sr/nm] proposed in (35) uses the filling-in of Fraunhofer lines caused by the plants’ chlorophyll fluorescence. Data were gridded to 16 day and 0.5° resolutions from the individual soundings and cover 11 years (2007–2017). No spatial smoothing or temporal averaging was performed before computing or averaging results. High sun zenith angle (SZA) observations (SZA > 70°) were removed from the analysis as well as cloudy scenes with a cloud fraction over 50% and observations taken between 2 p.m. and 8 a.m. local time. The illumination corrected SIF/cos(SZA) was considered, cf. section S6.

MODIS BRDF-corrected reflectances. MODIS reflectance data were derived from the MCD43A4.006 bidirectional reflectance distribution function (BRDF)-Adjusted Reflectance 16-Day L3 Global 500m product (36). They are disseminated from the Land Processes Distributed Active Archive Center (LP DAAC) also available at Google Earth Engine (GEE). The MCD43A2 MODIS product, which contains ancillary quality information for the corresponding MCD43A4 product, was also used for avoiding low-quality BRDF estimates. We computed the indices and conducted the analysis at 16-day temporal and 500-m spatial scales over the 11 years of SIF data.


General rationale. In all our experiments, we used reflectance values from MODIS, yet radiances or digital counts could also be used. The flux tower GPP estimates in our experiments come from the site-level data in (31). The SIF product comes from GOME-2, so the product is fully independent of MODIS reflectances. GPP and SIF correlations are computed in the time domain, while for SIF, we additionally compute correlations in space and then average results over time (results shown in section S6).

In all cases, we compute correlations between indices (NDVI, NIRv, and kNDVI) and the considered product only in meaningful vegetation classes: Needleleaf Forest, Evergreen Broadleaf Forest, Decidious Broadleaf Forest, Mixed forest, Shrublands, Savannas, Herbaceous, and Cultivated. These resulted from a meaningful grouping of International Geosphere-Biosphere Programme (IGBP) classes (see section S3). Analysis of the SIF results also considered aggregated climatic zones (Tropical, Arid, Temperate, Cold, and Polar), monthly means, and latitude averages (see section S6).

kNDVI calculation. The kNDVI index is defined askNDVI=k(n,n)k(n,r)k(n,n)+k(n,r)(1)where n, r ∈ ℝ refer to the reflectances in the NIR and red channels, respectively, and the kernel function k measures the similarity between these two bands. We used in all cases the RBF kernel, k(a, b) = exp (− (ab)2/((2σ2)), where the σ parameter controls the notion of distance between the NIR and red bands. This kernel function induces an important simplificationkNDVI1k(n,r)1+k(n,r)=tanh ((nr2σ)2)(2)

Other kernel functions are possible, but the RBF kernel is the most widely used one because of its theoretical and practical advantages (see sections S1 and S2) (24, 25). We calculated the kNDVI fixing the length-scale parameter σ equal to the mean distance between the NIR and red bands, σ = 0.5(n + r), which is a standard heuristic in the kernel methods literature, makes the index adaptive to each pixel, and worked very well in practice. Note that this simplification further reduces the index tokNDVI=tanh (NDVI2)(3)

Further optimization of σ per biome was done, but results did not improve substantially (results not shown).

Reproducibility: Open-source software and data. All calculations, visualization, and analyses were performed using the MATLAB programming language. We stored and processed netCDF files and tabular data. The kNDVI can be easily coded and applied. We give implementations in five standard programming languages (MATLAB, R, Python, Julia, and IDL) and in the GEE in section S10.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank J.L. Rojo-Álvarez and M. Martínez-Ramón for the many discussions on kernel methods, and L. Alonso who provided feedback on our logic. This work used eddy covariance data acquired and shared by the FLUXNET community, including the following networks: AmeriFlux, AfriFlux, AsiaFlux, CarboAfrica, CarboEuropeIP, CarboItaly, CarboMont, ChinaFlux, Fluxnet-Canada, GreenGrass, ICOS, KoFlux, LBA, NECC, OzFlux-TERN, TCOS-Siberia, and USCCC. Funding: G.C.-V. was supported by the European Research Council (ERC) under the ERC Consolidator Grant 2014 project SEDAL (647423). M.C.-T. and F.J.G.-H. were supported by the EUMETSAT Satellite Application Facility on Land Surface Analysis (LSA-SAF). SR research was financially supported by the NASA Earth Observing System MODIS project (grant NNX08AG87A). J.A.G. acknowledges the support of NASA ABoVE award number NNX15AT78A. S.W. acknowledges funding from the Emmy Noether Programme (GlobFluo project) of the German Research Foundation (GU 1276/1-1) as well as funding from the European Union’s Horizon 2020 research and innovation program under grant agreement 776186 (CHE project) and agreement 776810 (VERIFY project). Author contributions: G.C.-V. and M.C.-T. conceived and developed the methodology. S.W. provided harmonized SIF data, M.J. provided harmonized GPP data. M.C.-T. collected and harmonized the LAI dataset. G.C.-V., A.M.-M., and M.C.-T. did the experiments and performed the analysis with help from G.D., S.W., and L.G. J.M.-M. implemented the index in several languages, including the GEE platform. All authors analyzed the results and contributed to the writing of the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article