Land potential assessment and trend-analysis using 2000–2021 FAPAR monthly time-series at 250 m spatial resolution

General workflow

To model land potential we use the newly provided GLASS FAPAR V6 dataset as a baseline input (Ma, 2022). We estimated trends and gaps in FAPAR globally per pixel for the time period 2000–2021 using a spatio-temporal machine learning (ML) framework. An overview of the workflow is shown in Fig. 1. The whole processes can be represented with seven (7) main steps. First, we aggregated the GLASS FAPAR V6 8-day time series from 2000 to 2021 to monthly values of three percentiles (5th, 50th and 95th), representing their position in the value distribution. Second, we generate a sampling design and then overlay the points using the location of the points and the date to produce a regression matrix to model the actual FAPAR as a function of environmental predictors, including several variables related to human pressure. Third, we optimize and fit an ensemble machine learning model. Fourth, we use the fitted model to predict potential FAPAR, in this case by setting the value in the hypothetical state of no-human-impact, i.e., removing the impact of urbanization, intensive fertilization, and irrigation. Fifth, we calculated the gap between actual and potential FAPAR globally per pixel on a monthly basis for the year 2021. Sixth, we derived the long-term linear trend per pixel to supplement the gap analysis between potential and actual FAPAR, and finally (seventh) we ran a spatial overlay and match both the FAPAR long–term trend and the FAPAR gap with land cover maps to derive statistics per land cover type.

GLASS FAPAR V6

The GLASS FAPAR V6 data is a harmonized and consistent time series of FAPAR at 250 m resolution, which covers the time period from 26th February 2000 to 27th December 2021 in 8 day interval steps (Ma et al., 2022). It is an improved version of the GLASS FAPAR V5 in terms of time series trend stability. In contrast to V5, which makes use of radiative transfer functions, V6 is derived from a bidirectional long/short-term memory (Bi–LSTM) model to estimate FAPAR using Moderate Resolution Imaging Spectroradiometer (MODIS) reflectance and LAI data and several other FAPAR products (MODIS Collection 6, GLASS FAPAR V5, and PROBA-V1 FAPAR).

We derived monthly FAPAR time series by extracting, per pixel, the p0.05, p0.50, p0.95 percentiles using time series of the 8–day GLASS FAPAR V6 data. The percentiles indicate the FAPAR values that encompass the specified percentage of data set values falling below them. The 50th percentile (p0.50), hence, represents the median of the aggregated values. We processed all data in MODIS Sinusoidal Tile Grid and using full parallelization as implemented in the joblib library version 1.1.0 (https://joblib.readthedocs.io) and the scikit-map library (https://github.com/scikit-map/scikit-map) in Python. For each month, we aggregated all composites within that month plus one composite each before and after that month, resulting in five to six composites as input for each monthly value. This helped to reduce potential missing dates in the time series and noise in the original FAPAR 8–day composites and allowed a better understanding of differences within each month. Finally, we merged the processed tiles into global Cloud-Optimized GeoTiffs, which resulted in a total of 792 composites throughout 22 years. These data (almost 1.4 TB in size) are publicly available via https://doi.org/10.5281/zenodo.8381409.

We explore the monthly aggregated FAPAR time series by spatially overlaying all training points used for modeling with monthly aggregated time series data of the Enhanced Vegetation Index (EVI) derived from MOD13Q1 (50th percentile) at 250 m spatial (Didan, 2021), monthly daytime land surface temperature derived from MOD11A2 (50th percentile) at 1 km spatial resolution (Wan, Hook & Hulley, 2015), monthly water vapor derived from MCD19A2 (mean value) at 1 km spatial resolution (Lyapustin & Wang, 2022; Parente, Simoes & Hengl, 2023), and monthly accumulated precipitation from CHELSA V2.1 (Karger et al., 2017, 2021) at 1 km spatial resolution. To aid interpretation of FAPAR patterns, we created density scatter plots to visualize FAPAR (all three percentiles p0.05, p0.50, p0.95) against each variable and calculated Pearson’s r coefficient.

We validated our monthly aggregated GLASS FAPAR V6 values against reference measurements provided by Ground Based Observation for Validation (GBOV). The “RM6 —Fraction of Intercepted Photosynthetically Active Radiation” data was downloaded from https://land.copernicus.eu/global/gbov/ (Brown et al., 2020). The stations included in this analysis are shown in Table 1. For sites where measurements from both upward and downward facing DHP (Digital Hemispherical Photography) images were available (FIPAR–up and FIPAR–down, respectively), we approximated total FIPAR following Eq. (1), as provided in Brown et al. (2020):

(1) $$FIPAR=FIPA{R}_{up}+(1-FIPA{R}_{up})\ast FIPA{R}_{down}.$$

Note that the difference between DHP-derived FIPAR and FAPAR is generally considered marginal (Brown et al., 2020).

We filtered out any values with a measurement error greater than 0.1 or any raised quality issue (quality flag greater than 1). Each site has differing time intervals between measurements. For a more coherent comparison with the monthly GLASS FAPAR V6 values, we aggregated the reference measurement values per site to monthly averages whenever more than one value per month and year was available. This resulted in 3,262 unique space-time points. Reference measurements, as provided by GBOV, represent one elementary sampling unit each. Depending on the ground station, the extent of an elementary sampling unit can range from 20 m × 20 m to 40 m × 40 m (see the Algorithm Theoretical Basis Document on RM6 provided by the GBOV service at https://land.copernicus.eu/global/gbov/products/). We spatially and temporally matched the stations with GLASS FAPAR V6 values, assuming a relatively homogeneous vegetation cover in the 250 m radius around each station and hence a neglectable effect of differing spatial resolutions between reference and observed data.

Predictor variables

A total of 68 spatially explicit and harmonized variables representing (bio)climatic, topographic, geographic, vegetation cover and human pressure factors were used for modeling purposes. All layers were resampled onto a standard grid covering latitudes between 87.37°N and 62.0°S and reprojected to the coordinate reference system EPSG: 4326. The majority of covariates, for example, variables from the Digital Terrain Model (DTM) variables, lights at night, snow probability and similar, were available at 250–500 m spatial resolution (Yamazaki et al., 2017). Other 1–km resolution variables were resampled to 250 m resolution using a cubic splines in GDAL (Warmerdam, 2008). Whenever time-series covariate data were unavailable for the years 2020 and/or 2021, we used the most recent available year (e.g., previous year). We also added the percentile value of FAPAR as a covariate so that the fitted model can be used to predict FAPAR for an arbitrary part of the distribution.

Training points

We generated training points using a stratified random sampling approach, excluding areas of permanent water or ice (Brus, 2022). To ensure the inclusion of a minimum number of samples from proportionally underrepresented feature-space areas of human impact and of areas with least human impact, we created a set of non-spatially overlapping strata focusing on:

• 1. areas with minimal human pressure: all four human pressure variables at zero in years at the start (2000 and 2001), and at the end of the time series (2020 and 2021),

• 2. areas with recent change in vegetation cover from forest/wetland to cropland/grassland and vice versa comparing years 2000/2001 to 2019/2020,

• 3. areas with a change in nightlight values above plus or below minus 1[ $\mathrm{n}\mathrm{W}.\mathrm{c}{\mathrm{m}}^{-2}.\mathrm{s}{\mathrm{r}}^{-1}$] comparing the average of the years 2000/2001 to 2020/2021.

We calculated the optimal sample size per stratum based on the total standard deviation of each stratum of four FAPAR images across the quarters (January, April, July and October) of the years 2001 and 2021, using Neyman’s allocation formula (Neyman, 1934; Wright, 2014):

(5) $$n_h=\frac{N_h{S}_h}{{\sum }_{i=1}^H{N}_i{S}_i}n,$$and calculating sample weights $W_h$ using

(6) $$W_h=\frac{N_h}{n_h},$$where $n_h$ is the sample size for stratum $h$, $N_h$ is the population size, $S_h$ is the standard deviation and $n$ is the total sample size.

In summary, a total of 12,500 unique spatial points were sampled using this sampling design. A minimum sample size of 200 was applied to each stratum, and the remaining samples were distributed proportionally according to the standard deviation of each stratum. A full overview of the strata and sample sizes can be found in Table 3. The stratification of points with simple random sampling within each stratum resulted in globally distributed samples (Fig. 2). Each training point was overlaid with the time series of each of the three FAPAR percentiles (p0.05, p0.50, p0.95) and of all covariate layers from 2000–2021. This resulted in a total of 9,825,000 unique space-time data points. In this first model test, the cross-validation of the ANN model exceeded the volatile memory available in the machine used for the computation (i.e., 1 TB), hence we finally decided to randomly subset the 10 millon points to about 3 million space-time points from the total as input for modeling. Note that the regression matrices are also large because we used a large number of covariate layers.

Potential FAPAR model building and evaluation

FAPAR gap and trend analysis of time-series data

We calculated the monthly gaps between the actual and the potential FAPAR for the 95th percentile per pixel for the year 2021 by subtracting the potential FAPAR predictions from the actual FAPAR values. We spatially matched the FAPAR gap with every combination of land cover between 2000 and 2020 according to the ESA CCI annual land cover maps, stable and transition areas (ESA, 2017) and visualized the distribution of gap values for each. We conducted the analysis in the aggregated land cover classes as shown in Table 5.

To perform trend analysis, we extracted the trend component of the FAPAR time series (2000–2021) by season–trend decomposition using locally estimated scatterplot smoothing (LOESS), Seasonal and Trend decomposition using Loess (STL) (Ben Khalfallah et al., 2021), and fitted an ordinary least squares (OLS) model on the trend component. For every FAPAR pixel we computed the OLS coefficients (alpha and beta), standard deviation, $p$ and $R^2$. The percentage of area with significant positive and negative trends ( $p<0.005$) was calculated for every land cover change combination (2000 and 2020) according to aggregated ESA CCI land cover classes.

Monthly aggregated FAPAR data exploration

The correlation between the monthly aggregated actual FAPAR in relation to EVI, precipitation, daytime land surface temperature (LST) and water vapor, together with the calculated Pearson’s r for each, is shown in Fig. 3. FAPAR shows a strong positive correlation with EVI (Pearson’s $r=0.89$). The correlation with monthly precipitation and water vapor is also positive but less strong. The density plots indicate a more complex distribution of these two parameters in relation to FAPAR. Water vapor and precipitation both show a similar range of values throughout the FAPAR range. The density plot of LST against FAPAR indicates a decrease of LST with increasing FAPAR. However, this is not reflected in the Pearson’s correlation value, which can be attributed to a wide range of values in the LST especially at a FAPAR value of 0, where values range from −50 °C to +50 °C. Generally, FAPAR values in the range 0.75–0.85 were underrepresented compared to the entire FAPAR range. In general, it can be said that FAPAR is strongly correlated with EVI, LST, precipitation, and water vapor, but the relationship is often non-linear and saturation can be observed for FAPAR values near to 0.0 and 0.9.

The validation of the monthly FAPAR shows a strong agreement and a low bias with the GBOV ground validation points ( $R^2=0.80$, Fig. 4). This matches the results of Ma et al. (2022).

FAPAR model performance and variable importance

The recursive feature elimination showed an optimal feature number of 52 (out of 68 variables provided). The spatial cross-validation results for modeling distribution of monthly FAPAR with the selected 52 variables using the EML model showed an R² of 0.89 and a CCC of 0.94. RMSE and MAE resulted in moderate values of 0.10 and 0.06, respectively. The selected hyperparameter sets for the “Extremely randomized trees” and “Gradient descended trees” model are shown in Table 4. The accuracy metrics for each submodel are only slightly lower, showing that each model performs similarly well (Table 6).

The bias for all tested models based on five–fold cross validation (CV) is close to zero. The evaluation of the EML model against independent validation set (500,000 spacetime points) resulted in similar values to the spatial CV results with an R² of 0.90, CCC of 0.95, RMSE of 0.09 and MAE of 0.06 (Fig. 5). The distribution of the density plot of reference against prediction for this independent point data set from GLASS FAPAR V6 shows an overall high agreement between the reference and predicted values.

The predicted actual FAPAR at ground stations also showed a high agreement with the GBOV points (R² = 0.75, Fig. 4). The results of both independent tests show that high FAPAR values (above 0.75) tend to be underestimated, while low FAPAR values (below 0.4) tend to be overestimated (Fig. 5).

Variable importance analysis showed that the three most important variables for modeling were: (1) the growing season length, (2) the forest cover indicator, and (3) the annual precipitation amount (BIO12) as shown in Fig. 6. In general, the climatic variables played a more significant role in explaining FAPAR than topographic and lithologic variables. All variables on human pressure (human footprint index, population count, cropland intensity, and nightlights) were retained in the recursive feature elimination, highlighting the additional information each of them provides in explaining FAPAR, despite appearing lower in the variable importance ranking than climatic variables. While both the “Extremely randomized trees” and the “Gradient descended trees” model showed a higher importance for the human footprint index than other human pressure variables, the “Gradient descended trees” gave cropland intensity, population density and nightlights a higher importance than “Extremely randomized trees”.

Visual examination of global maps of FAPAR gap actual vs. potential and long-term trend

We visually inspected the global maps of average annual actual and potential FAPAR produced for the year 2021 (Fig. 7). The global distribution of potential FAPAR indicates that Northern latitudes (Europe, North America, Northern Asia) show mostly negative gap values (actual FAPAR lower than potential), while highest values are found at latitudes close to the equator, especially in the Amazon region, the Congo River Basin and South-East Asia. The model deviance (standard deviation of the sub-models) of these regions is comparably low with other regions of the globe.

The gap of actual vs. potential FAPAR shows that the area of these same regions exhibits mostly little difference, but with patches of potential FAPAR lower than actual. Regions with a potential FAPAR lower than the actual relate mainly to locations that show a potential forest cover according to the potential biome maps. On the other hand, regions with a potential FAPAR higher than the actual relate to the distribution of potential grassland and shrub cover.

The model deviance seems to be higher for regions that show a potential FAPAR lower than the actual FAPAR than vice versa (Fig. 7). Predictions of potential FAPAR in deforestation areas in Brazil show a potential FAPAR close to the surrounding forested areas, although generally with slightly lower values. This is expected from the model results, which show a tendency to underestimate FAPAR in the range 0.75–1.0 (Fig. 5). We specifically looked at known areas of irrigation to identify the behavior of the model in these areas. The typical pivot irrigation areas in the U.S. and the region around the Nile in Egypt show a potential FAPAR lower than the actual FAPAR. This shows that the model detects human pressure variables (e.g., cropland intensity) as an important factor in predicting FAPAR.

The global FAPAR trend map for 2000–2021 shows that regions at latitudes close to the equator show large areas of negative trends (Fig. 8). Negative trends are also observed on a wide scale in Australia, North America, and Central Asia. Europe, on the other hand, shows mainly stable and positive trends. Both India and China show especially wide areas of positive FAPAR trends. While in India we observe mostly an actual FAPAR higher than its potential with a high model deviance of the potential FAPAR, in China we observe mostly an actual FAPAR lower than its potential with a low model deviance.

For three sample locations in Figs. 8–10, we provide a more detailed inspection of FAPAR trend, gap and time series. The location of the reforestation project at Las Gaviotas in Colombia shows a spatial overlap of positive FAPAR trend with reforestation and a decrease in LST (Fig. 8). Although the actual FAPAR shows lower values than its potential at the start of the time series, it increases over time to similar values as the potential and even grows beyond the potential values (Fig. 9). The FAPAR maps show that the potential FAPAR is similarly high across the region, while the actual FAPAR is lower in the vicinity of the reforestation area than the potential FAPAR (Fig. 10). The location close to the Iriri River in Brazil shows a spatial overlap of the negative trend of FAPAR with deforestation and an increase in LST. While the potential FAPAR remains stable over time, the actual FAPAR decreases from 2008 to 2016. Similarly, as for the site in Colombia, we can see a homogeneously high potential FAPAR in the area. The gap between actual vs. potential FAPAR is small in the stable forest area, while there is large gap in the deforestation part.

For the location of Kismaayo in Somalia, we observe a strong negative trend over the entire inspected area. Land cover remains stable over the period of time; however, an increase in LST observed, which is mostly observed in areas of “other vegetation” (shrubland in this case), while cropland shows a more stable LST over time. The centre pixel shows a decrease in the time series, whereas the potential FAPAR remains rather stable. The actual FAPAR is mostly lower than the potential; however, the model deviance is higher in this area than in the other two sample locations.

FAPAR gap and trend analysis with land cover

Our spatial analysis of the relationship between the average FAPAR gap for the year 2021 and land cover stable and change classes shows that the largest negative gaps (actual FAPAR lower than potential) are found in the stable classes ‘Lichens and mosses’ (LICMO) (90%), ‘Urban areas’ (URBAN) (86%), ‘Tree cover needleleaved deciduous’ (TREND) (67%), and ‘Shrub or herbaceous cover flooded’ (SHHFL) (61%); and in the change classes ‘Cropland rainfed-Urban areas’ (CRPRF-URBAN) (82%), ‘Tree cover needleleaved evergreen-Shrubland’ (TRENE-SHRUB) (72%), ‘Tree cover needleleaved deciduous-Mosaic tree and shrub or herbaceous cover’ (TREND-MTSHH) (72%), and ‘Tree cover needleleaved evergreen-Mosaic tree and shrub or herbaceous cover’ (TRENE-MTSHH) (71%) (see Fig. 11). The largest positive gaps (actual FAPAR higher than potential) are found in the stable classes ‘Tree cover broadleaved deciduous’ (TREBD) (64%), ‘Shrubland’ (SHRUB) (60%), and ‘Cropland, irrigated or post-flooding’ (CRPIF) (57%); and in the change classes ‘Tree cover broadleaved deciduous-Shrubland’ (TREBD-SHRUB) (87%), ‘Shrubland-Tree cover broadleaved deciduous’ (SHRUB-TREBD) (84%), ‘Shrubland-Cropland rainfed’ (SHRUB-CRPRF) (66%), and ‘Mosaic cropland or natural vegetation-Tree cover broadleaved deciduous’ (MCRNV-TREBD) (53%). Especially classes ‘Water bodies’ (WATER) and ‘Urban areas’ (URBAN) show a wider distribution of the gap values, while classes ‘Lichens and mosses’ (LICMO) and ‘Tree cover needleleaved deciduous’ (TREND) show a more narrow dispersion.

Regarding the relation between the FAPAR long–term trend 2000–2021 and land cover change grid, especially high positive values are found for stable classes ‘Cropland, irrigated or post-flooding’ (CRPIF) and ‘Tree cover mixed leaf type’ (TREMX); and change classes ‘Mosaic cropland or natural vegetation-Tree cover broadleaved evergreen’ (MCRNV-TREBE) and ‘Mosaic cropland or natural vegetation-Mosaic tree and shrub or herbaceous cover’ (MCRNV-MTSHH) (see Fig. 12). Strong negative trends are found for stable classes ‘Sparse vegetation’ (SPARV) and ‘Urban areas’ (URBAN); and change classes ‘Tree cover broadleaved evergreen-Mosaic cropland or natural vegetation’ (TREBE-MCRNV) and ‘Tree cover broadleaved evergreen-Cropland rainfed’ (TREBE-CRPRF). Land cover classes that show the largest percentage of positive trends within the class are ‘Cropland, irrigated or post-flooding’ (CRPIF) (85%), ‘Tree cover mixed leaf type’ (TREMX) (82%), ‘Tree cover broadleaved deciduous’ (TREBD) (80%), and ‘Tree cover needleleaved deciduous’ (TREND) (79%). Land cover classes that show the largest percentage of negative trends within the class are ‘Sparse vegetation’ (SPARV) (33%) and ‘Urban areas’ (URBAN) (31%). In total, we found that 19% of the global land area shows a decreasing trend, while 67% shows an increasing trend.

The key variables explaining spatial and temporal variation in FAPAR

The results of modeling actual FAPAR and variable importance assessment in our work clearly show that the key variables that explain the variation in the monthly FAPAR are the length of the growing season, the forest indicator and the annual precipitation. Our results are in line with the findings of Hengl et al. (2018) on modelling FAPAR, who applied a similar methodological approach to model FAPAR for the years 2014–2017 based on machine learning and environmental covariate layers. Especially the fact that precipitation is one of the key explanatory variables is confirmed in our study; however, we systematically extended the number of covariate layers and also achieved an overall higher model accuracy.

In our initial data exploration of the monthly aggregated FAPAR and its relationship with EVI, land surface temperature, and water vapor, we found expected relationships between FAPAR and the corresponding variables, such as a decrease in land surface temperature with increasing FAPAR. EVI has been widely used as a measure of vegetation density and an indicator of primary productivity (Brown et al., 2020; Jay et al., 2016), therefore a close relationship with FAPAR was expected. However, the wide range of present values in all density plots (see Fig. 3) indicates that the estimation of monthly FAPAR from other biophysical MODIS variables is complex. Both monthly FAPAR, EVI and LST seem to be unique biophysical variables and are all important for global modeling efforts.

Accuracy of the biophysical models explaining FAPAR

Our results further show that FAPAR can be modeled to a maximum R² of about 0.9 (based on strict cross-validation). This means that there is still a significant part of the variation (about 10%) that we could not map on this scale of work. Our validation of the aggregated monthly FAPAR against ground reference measurements shows agreement with the validation of GLASS FAPAR V6 by Ma et al. (2022) against ground measurements, who reported an R² of 0.8, but for their 8–day product. The validation of our FAPAR model predictions show high agreement with ground reference measurements as well ( $R^2=0.75$). The close match between the weekly and monthly model accuracies is possibly due to the relatively smooth variation of FAPAR.

Validation of FAPAR values in Ma et al. (2022) using ground reference measurements (GBOV) showed an overestimation of FAPAR at low values. This is also in agreement with Brown et al. (2020), who reported an overestimation of FAPAR derived from remote sensing of sparse vegetation. Additionally, the results of the validation of Ma et al. (2022) confirm the overestimation of FAPAR at forest sites, which may be related to the saturation of FAPAR in dense canopies.

The ensemble and the individual base models showed comparable model performance in our study. The ensemble approach, however, offers the main advantage of providing the model deviance, i.e., the standard deviation of the base models, per pixel. In addition, using an ensemble model could potentially increase the spatial consistency of uncertainty as some base learners may perform better or worse in specific feature spaces than others, as has been shown for land cover classes by Witjes et al. (2022). While the tested hyperparameter sets were limited, we achieved high accuracy. Nonetheless, testing a wider range of values could potentially improve results even further.

Visual inspection of the differences between potential and actual FAPARs (around the world) largely matches our expectations. This indicates that our models can be considered stable and could potentially also be used to predict future FAPAR values, e.g., given different climate change scenarios or similar (as in Bonannella et al. (2023)). Similarly, in future research, it will be interesting to model and understand the key drivers of the change in FAPAR in a complex Earth system with changing climate, changing policies, and shifting vegetation.

Mapping gaps between actual vs. potential FAPAR

Our results of comparing actual vs. potential FAPAR show that the biggest gaps between the two are found for the following land cover classes:

• Potential FAPAR bigger than actual (negative gap): urban, needle-leave deciduous trees, flooded shrub or herbaceous cover, and lichen and mosses.

• Potential FAPAR smaller than actual (positive gap): broad–leaf deciduous trees, shrubland, irrigated or post–flooded cropland and broad-leaf evergreen trees.

In the above section, we only list statistics based on the annual average of 2021, and monthly values might differ. Analysis of monthly values and mapping of gaps is complex. More robust visualization tools and analyzes are needed to further detect gaps and patterns over time in these complex data.

Also note that in this article we used the concept of “Potential Natural Vegetation”, which requires some discussion. With PNV concept we often assume that “natural” is the optimal state of the ecosystem. In reality, there will be many situations (especially where people intensively irrigate the land and make intensive fertilizer applications), where a human-managed or human-designed landscape can show much higher productivity than a “natural” one. Similarly, the actual FAPAR greater than its potential value does not necessarily refer to land recovery or improvement in land degradation status (Bai et al., 2008; Cherlet et al., 2018). We also observed, for example, that savannas often have a lower potential FAPAR than the actual one. While there is controversy about the term of what is considered “natural”, we believe that it is still a useful exercise to quantify the indirect effects of removing the impacts of urbanization, intensive fertilization, and irrigation on potential FAPAR.

Relationship between long-term trends and gaps in FAPAR in relation to land cover changes

In the last 20+ years, researchers have, without doubt, discovered that the planet is becoming greener. Especially China and India are showing a large overall increase in greenness due to extensive agricultural applications (Chen et al., 2019). However, Gao et al. (2023) have analyzed NDVI time-series data for 1982–2015 and have found notable increases in interannual variability (IAV) of NDVI over the northern high latitudes, Eastern Europe, and Central Australia, indicating that parts of the world are becoming more sensitive to climatic, environmental and anthropogenic changes. The greening of the Earth in the last two decades is likely driven by CO₂ fertilization (Chen et al., 2022), N deposition, climate change, and change in land cover (Zhu et al., 2016). However, croplands in China and India appear to have a decrease in the IAV of vegetation greenness, suggesting more controlled management practices through inter-cropping, irrigation and fertilization. Likewise, there is an increase in frequency, severity, duration, and scope of global droughts, so it is likely that the greening of the planet will not continue to compensate for the increased IAV and may even reverse in the coming years (Chen et al., 2022).

Our results also show that, for most of the planet (67%), there seems to be a positive increase in FAPAR over the time period 2000–2021. The results of matching the FAPAR trends (betas) with land cover change categories show that stable vegetation cover classes with the largest percentage of negative trends are: sparse vegetation, rain–fed cropland, and mosaic croplands. This could be expected, since these land cover types are potentially especially sensitive to changing climate, such as increased droughts. The most negative mean beta found was −0.14. The highest percentages of negative gaps are found in urban areas and in the vegetation classes: lichens and mosses, needle-leaf deciduous trees, and flooded shrub or herbaceous cover.

In summary, the planet over the 2000–2021 time period has on average been getting greener, but there are many local areas, especially in the tropics and areas of savanna and steppe, where there are distinct decreases in FAPAR over the last 20+ years. Attention needs to be paid to these areas to support land owners and governments in their efforts to restore, reverse and reduce degradation. Despite a general trend toward greening, we can see large gaps in the actual FAPAR around the world.

Limitations of this study

Although the FAPAR product by Ma et al. (2022) is to our knowledge the most comprehensive and longest period product at 250 m spatial resolution, the data set also comes with some limitations:

• 1. The resolution of images is relatively coarse, which means that many land use systems, especially in regions where crop fields are relatively small, are combined/mixed-in, making it difficult to untangle driving factors. Moving to using Landsat as the primary source of the FAPAR assessment could help produce data with a minimum number of “mixels” (pixels of mixed land use) and increase compatibility with the best global land cover products, e.g., Zhang et al. (2021).

• 2. It is well known that optical remote sensing sensors have limited accuracy in tropics due to clouds and saturation over dense canopies; hence there is possibly a significant difference in accuracy of all products we produced across the globe, with the tropics probably being the most difficult to track. To avoid effects of clouds and water vapor on estimates of FAPAR we used the p0.95 monthly values; however, in some tropical areas are >80% of time under cloud cover and, hence, probably not even the p0.95 monthly value helps.

• 3. GLASS FAPAR V6 is derived itself from a machine learning model (Bi–LSTM) based on remote sensing data, and hence introduces an inevitable uncertainty that will propagate to the prediction of potential FAPAR. However, comparison with GBOV ground reference measurements showed that our model predicts actual FAPAR with relatively high accuracy, despite the fact that we did not use surface reflectance values/indices.

Verger et al. (2023) have recently produced smoothed and gap-filled time series of global FAPAR at a spatially coarser resolution of 1 km. This data set overlaps a great deal with the work of Ma et al. (2022) and could potentially also be used to perform a similar analysis. However, we focused in our study on finer resolution data in order to attempt to capture small-scale changes, for example, of land use systems with individual land parcels that are relatively small in size and fractured. In future work, modeling potential FAPAR at even finer spatial resolutions of 30 m will be useful to match the recent highest resolution dynamic land cover products (Potapov et al., 2022; Zhang et al., 2021).

Regarding our modeling results, we highlight the following limitations of the methodology and interpretation of the results:

• 1. Our cross-validation results indicate that there is still a large component of variation (about 10%) in monthly FAPAR that we can not explain. This is probably a limitation of the moderate spatial resolution we worked with, including using covariates at coarser resolution (1 km) than the FAPAR product. In addition, variable importance can be related to the training sample distribution, possibly under-representing areas of some of the land form or lithologic variables, which are less present on the globe.

• 2. Although the EML model deviance can give an indication of the uncertainty of the prediction, the lack of ground-truth data on potential FAPAR requires the results to be interpreted cautiously. A more robust framework is needed to calculate the probability distributions of FAPAR per pixel (Andriuzzi et al., 2022).

• 3. For the prediction of potential FAPAR, we make use of a cross-walk table from the ESA CCI land cover towards the BIOME 6000 classification. This may introduce a bias in the prediction. Furthermore, it results in distinct potential FAPAR boundaries between adjacent biomes, which are unlikely to be so distinctly present in reality.

• 4. We refer here to the potential FAPAR as the FAPAR under no-human pressure by removing the impact of urbanization and intensive agriculture. This refers to the least human pressure under current conditions of the month we predict on. We again highlight the “shifting baseline syndrome” (Papworth et al., 2009), where despite using a global coverage of training points, the results may also be subjected to this syndrome. We cannot preempt that there is human pressure present in areas that are not indicated by our variables, and there may be legacy effect of preceding human pressure (since likely most of the globe has been subjected to human influence at some point).

• 5. Climate data sets we used are only an approximation of true conditions; climate depends also on present vegetation type and will change if there is a forest where there has been none before, and vice versa; therefore, separation of climate and vegetation is not trivial. Here we assumed that climate is an independent variable for FAPAR modeling, but this is likely an oversimplification, especially at finer resolutions.

• 6. FAPAR does not just depend on the current month climate/precipitation, but also on preceding conditions, e.g., the cumulative precipitation of 2–3 months before. We have so far ignored this aspect in this study, as it would increase the computational complexity and would require additional testing.

• 7. The most important variables in our results were not human pressure variables, which overall have small influence on the prediction, although clearly shown in irrigated areas. The forest cover indicator, on the other hand, had a high importance for the prediction of potential FAPAR.

• 8. Vegetation cover variables based on ESA CCI land cover change maps are at 300 m resolution. We noticed that GLASS FAPAR V6 picked up a deforestation signal at some sites in Brazil (Fig. 9), but ESA CCI land cover did not show a change, which also affects the model training.

Finally, there are two global vegetation modeling paradigms in the literature: (1) the Dynamic Global Vegetation Models (DGVMs), and (2) data science models based on fitting ML models to training data. DGVMs such as the Lund–Potsdam–Jena-DGVM, have been used to simulate potential primary productivity by excluding the land-use component from the process-based model (Haberl et al., 2007). DGVMs are generally restricted by the simplified mechanistic rules set in the model and the definitions of vegetation (Scheiter, Langan & Higgins, 2013). Data-driven ML approaches can be used to estimate the possible potential state of vegetation based on correlations with biophysical variables (Hengl et al., 2018). In this article, we used a data-driven approach, but this comes at the cost that it is purely correlative, e.g., based on an extrapolation of observed relationships and not on modeling actual processes.

Implications of this study and next steps

With the FAPAR analysis produced and a global monthly map of the potential FAPAR, we have tried to contribute to international programmes such as Land Degradation Neutrality (UNCCD) and the Sustainable Goals of the United Nations (UN SDGs). We have released our data as open, and our code (used to produce models, visualizations, and predictions) via GitHub (https://github.com/Open-Earth-Monitor/Global_FAPAR_250m). We hope that this will inspire other groups to review our results and expand our research in new directions. We especially hope that national and local agencies will find these data as a good reference for spatial planning and that many small medium enterprises will see an opportunity for practical applications and products.

Although this work has taken significant time and computing resources, many sections of modeling and assessment could be further improved. Some potential extensions of this research include:

• 1. For a more thorough understanding on the impact of modifying the human pressure variables for the potential FAPAR prediction results, a sensitivity analysis of the model could be performed, i.e., analysing the effect each of them has individually but also in synergy with each other. The interaction between human pressure variables and climate variables could be investigated to better understand the magnitude of uncertainty in potential FAPAR.

• 2. In future studies, it would be interesting to investigate ways to quantify the “shifting baseline syndrome”, and to include these estimates in the modeling of potential FAPAR. As an example approach, we could think of finding ways to model FAPAR for historic times and use these values as a reference (Enquist et al., 2020). In addition, it would be interesting to compare our results with an approach not reliant on any hard class distinction of biomes or land cover types to account for the spatial gradient change of vegetation cover. One way could be to include fractions of plant functional types, for example recently made available by ESA CCI (Harper et al., 2023).

• 3. Due to limited computing resources, we had to subset the globe to a set of training points as input for modeling. It would be interesting to increase the size of the training points and compare the results with our current model. In our approach we made use of a linear regressor as meta-learner. The deployment of other meta-learners, such as logistic regressors, could potentially improve results further. In addition, to get a more robust way to estimate uncertainty of the EML predictions, it would be interesting to test the conformal prediction (Johansson et al., 2014) for this purpose.

• 4. For the land cover change analysis, we compared the years 2000 and 2020 to derive change classes. Some areas might experience multiple changes within this time series, which could be further investigated in following studies.

• 5. To derive robust conclusions on the degree of land degradation, additional information from environmental variables should be included to supplement the data on land productivity, such as land cover change or soil organic carbon. Shao et al. (2024), for example, made use of specific variables, such as wind speed, known to drive land degradation in their local study area. Our analysis of land productivity could be used as one of the variables in such frameworks that make use of the “convergence of evidence” to indicate potential degradation areas (Gianoli, Weynants & Cherlet, 2023).

• 6. Although our results give estimates of the gap between actual and potential FAPAR for current times, this approach could also be applied to future climate simulations. With estimates of future potential FAPAR, we could include this information already today in restoration, mitigation, and adaptation planning.