This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution license (http://creativecommons.org/licenses/by/3.0/).

Long-term trends in photosynthetic capacity measured with the satellite-derived Normalized Difference Vegetation Index (NDVI) are usually associated with climate change. Human impacts on the global land surface are typically not accounted for. Here, we provide the first global analysis quantifying the effect of the earth’s human footprint on NDVI trends. Globally, more than 20% of the variability in NDVI trends was explained by anthropogenic factors such as land use, nitrogen fertilization, and irrigation. Intensely used land classes, such as villages, showed the greatest rates of increase in NDVI, more than twice than those of forests. These findings reveal that factors beyond climate influence global long-term trends in NDVI and suggest that global climate change models and analyses of primary productivity should incorporate land use effects.

Climate change research often focuses on long-term productivity trends in vegetation at regional [1–3] to continental and global scales [4,5]. A key variable in these studies of vegetation productivity is the satellite-derived Normalized Difference Vegetation Index (NDVI), that is either used directly as a measure of photosynthetic capacity [6], or serves as an input variable for climate change models predicting terrestrial net primary production (NPP [5,7,8]). On a global scale, these studies often ignore that NDVI may not be driven by climate alone. Because the majority of the Earth’s ice-free land surface has been altered by human action [9–12] and because these same areas account for ~90% of terrestrial net primary production [13], trends in NDVI may also be linked to human land-use practices, such as land conversion, irrigation, and nitrogen deposition. Although several regional and local analyses demonstrate significant impacts on NDVI trends from changes in human land-use practices [14–16], global-scale analyses relating trends of NDVI to these direct anthropogenic effects have received little attention [5].

Here, we provide one of the first global analyses of direct anthropogenic effects on long-term trends in NDVI from 1981 to 2010 (however see [17] for an analysis that also include land use but on a different scale). First, we quantified variation in NDVI trends across “anthromes” (anthropogenic biomes, after [13]). Anthromes were introduced by Ellis and Ramankutty [13] and classify the earths land surface based on a combination of data on population density, land-use, and land cover [13]. They characterize land surface based on global patterns of sustained, direct human interaction with ecosystems. Comparisons of trends in NDVI among anthromes thus represent a useful probe of direct anthropogenic effects on long-term trends in NDVI. We also examined the effect of human-population density, as a surrogate for the degree of anthropogenic effects. Lastly we tested for specific mechanistic linkages that may underlie spatial variation in trends of NDVI among ecoregions and investigated effects of human land conversion (measured as percentage village, urban areas and croplands per ecoregion), nitrogen deposition, and irrigation.

2.Results and Discussion

Globally, anthromes with intensive land use (dense settlements, villages, and croplands) all had significant (p < 0.05) and steep increases in annual mean NDVI over a 29-year period compared to other anthromes. Forests exhibited marginally significant increases (β = 0.0006, p = 0.06), whereas we found no significant trends for wildlands or rangelands (Figure 1A,C,E,G,I,K). We emphasize that wildlands in the anthropogenic biomes classification used here [13] simply lump all land surfaces that are free or almost free of any human footprint. They thus comprise extremely varied ecosystems across the globe where some may show positive and other may show negative trends. Our statement only refers to the average across all these diverse ecosystems.

We also calculated robust linear regression Theil-Sen (TS) slope estimators for each spatial location in the global NDVI data (using the GIMMS3g NDVI data). Across all anthromes, the medians of the per-pixel TS estimators agreed closely with the slope estimates of the annual mean (compare slope estimates in Figure 1(left) with medians in Figure 1(right)). Based on the average NDVI value in the respective anthromes, the TS estimators showed that NDVI of villages changed the most over the 29 year study period (5.9% increase), followed by dense settlements (4.3%) and croplands (3.7%). With a 2.5% increase, global net NDVI of forests increased less than half as much as villages. Global net NDVI increases in wildlands and rangelands were close to zero.

We show that in years where global NDVI decreased, the global average NDVI of land types such as villages or croplands appear to have dampened minima compared to wildlands or forests (compare minima in Figure 1A,C,E,G,I,K). Variation induced by seasonal changes in NDVI, such as those associated with the interplay between human activities and monsoon dynamics [18], complicate interpretation of this pattern.

In our second analysis, we tested for the effect of human-population density, as a surrogate for the degree of anthropogenic effects, and found strong linkages between trends in NDVI and human population density. Here, we used ecoregions as spatial units [19]. Ecoregions are distinct with respect to natural communities and resident species and thus ideal to test for spatial variation in anthropogenic effects on NDVI. In ecoregions with low human population densities, NDVI trends were small and close to zero, whereas in ecoregions with high population densities, NDVI increased significantly over time (Figure 2). For example, densely populated ecoregions with an average population density of ~500 people/km2 showed an average increase of ~0.025 NDVI units over the 29 year study period (a 6% increase compared to the global mean). Except for Australia, average annual increases in NDVI were also greater in more densely populated ecoregions when data were analyzed at the continental scale. Europe, which has the greatest average population density across ecoregions, showed the strongest relationship between trends in NDVI and population density (Figure 3).

Lastly, we tested for specific mechanistic linkages that may underlie spatial variation in trends of NDVI among ecoregions. In particular, we tested for human land conversion (measured as percentage village, urban areas, and croplands per ecoregion), nitrogen deposition, and irrigation and found that all three of these factors had significant effects. The percentage converted lands per ecoregion, the percent irrigated area per ecoregion, and the amount of nitrogen deposition all were positively related to increasing trends of NDVI (Table 1). Together these three factors explained 21% of the variation in trends in NDVI.

These results were robust to several methodological concerns and alternative hypotheses. For example, analyses excluding the moist tropical biome yielded similar results, suggesting the uncertainty of remote sensing estimates in the moist tropics [20] did not drive our findings (Table 2). Results were also robust to the hypothesis that humans have settled preferentially in areas with high NDVI (alternatively, that areas with high NDVI had greater absolute increases in NDVI) because standardizing the per-pixel TS estimators by the mean NDVI did not significantly change model outcomes (Table 2). Finally, results were robust to the spatial unit of analyses, because substituting the biologically defined ecoregions with square grid cells of length 500 km also did not substantially change model outcomes (Table 2).

Ecoregions with the fastest increases in NDVI were predominately located in temperate Europe, the Mediterranean, the Sudano-Sahelian belt in Africa, western India, and northern China (Figure 4A). Ecoregions with the fastest decreases were predominately located in boreal forests of North America, as well as dry forests in South America (mostly Paraguay) and drylands of northeastern Australia. When comparing ecoregions with the fastest NDVI increases (the 5% of land surface with the fastest increases in NDVI, n = 73 ecoregions) against those showing the largest NDVI decreases (the 5% of land surface with the fastest decreases in NDVI, n = 38 ecoregions), we found that the ecoregions exhibiting the fastest NDVI increases had a greater average population density (median of 94 vs. 7 persons per km2, Figure 4), greater percentage of converted lands (88% vs. 9%, Figure 4B), greater percentage of irrigated lands (49% vs. 4%, Figure 4B), and greater deposition of nitrogen (573 vs. 209 mg N/m2/year, Figure 4B).

Ecoregions with the fastest changes in NDVI corresponded well with decreases and increases in NDVI shown in other global and regional studies [5,21]. Decreasing trends in NDVI have been reported for ecoregions with the fastest decreases, like the Dry Chaco region in South America [14], northeastern Australia [22], and boreal Canada and Alaska [23,24]. Likewise, increasing trends in NDVI have been reported for ecoregions with the fastest increases across Europe [25,26], northern China [27], central India [2], and the Sudano-Sahelian belt in Africa [28,29].

Although trends in NDVI in some of these ecoregions, especially those with negative trends in remote areas like boreal Canada and Alaska or northeastern Australia, have been predominately attributed to climatic effects [22–24], there is growing evidence for direct anthropogenic effects in almost all regions with significant human populations. For example, increasing trends in NDVI in India and northern China have been related to increasing rates of fertilization and irrigation [2,28]. Likewise, while previous studies have chiefly related the Sahelian greening trend to climate, e.g., [7], new investigations point to land conversion to croplands as one of the main causes for changes in NDVI [29]. In Europe, nitrogen deposition may be one of the main drivers of carbon sequestration in forests and—together with reforestation and climate changes—may well explain the large increases in NDVI observed across the continent [25,26,30]. We explain 35% of the variability in trends in NDVI by population density alone (Figure 3D), supporting the hypothesis that anthropogenic drivers are contributing to NDVI increases in Europe. There are also examples for direct anthropogenic effects in ecoregions with rapidly decreasing NDVI. Results of a recent study on deforestation undertaken by the UN-FAO [31] show that forest loss in South America accounts for more than 50% of the global forest area loss in the past two decades. One region particularly affected is the Dry Chaco in Paraguay, a region with rapid land conversion from forest to soybean plantation [32], identified by our analyses as an extremum for decreasing NDVI.

3.Experimental Section3.1.NDVI Data and Pre-Processing

We used NDVI data from the Advanced Very High Resolution Radiometer (AVHRR) instrument developed into the GIMMS3g dataset (Global Inventory Modeling and Mapping Studies) [33,34]. Previous versions of the GIMMS data set have been used extensively for testing temporal trend in NDVI [35], have been tested against ground controls points [36] and are thought to be the best long term data set for analyzing long term trends in vegetation dynamics [37]. The current GIMMS3g dataset represents an improvement over previous AVHRR NDVI data sets because it was coprocessed with SeaWiFS and MODIS Terra data [20]. The SeaWiFS and MODIS instruments have spectral bands specifically designed for vegetation monitoring, better atmospheric correction algorithms, reduced geometric distortions, and improved radiometric sensitivity [38]. MODIS NDVI has been collected over the past 13 years from the Terra platform and the past 11 years from MODIS on the Aqua platform. For the respective MODIS overlap periods, the per-pixel trends in NDVI from GIMMS3g agree well with those from Terra and Aqua MODIS, except in the moist tropics, which are difficult for remote sensing due to humidity and cloud cover [20]. MODIS Terra and Aqua NDVI data do not compare well for the tropics either. Using the longer-term GIMMS3g dataset, we conducted the analyses discussed below twice, once with all terrestrial pixels and once excluding the moist tropics. Since our analyses used yearly means (see below) and the studies cited above did not compare yearly means of NDVI between MODIS and GIMMS3g, we performed additional analyses comparing annual means between MODIS terra and Aqua, and the GIMMS3g data set. We sampled 34 1-degree areas across the globe and found a high linear correlation between the MODIS and GIMMS data sets. A Gaussian mixed model between annual means of MODIS Terra and GIMMS3g accounting for sampling area with a random effect showed very close relationships between both data sets (β = 0.74, t = 24.57, N = 298, coefficient det. = 0.99) as did the same comparison with MODIS Aqua (β = 0.75, t = 25.79, N = 247, coefficient det. = 0.99). If the data was standardized by the mean NDVI across years of each area to make areas more comparable, the relationship was still strong. MODIS Terra vs. GIMMS3g: (β = 0.68, t = 21.25, N = 298, coefficient det. = 0.65) and MODIS Aqua vs. GIMMS3g: (β = 0.69, t = 22.62, N = 247, coefficient det. = 0.65).

The GIMMS3g dataset of bi-monthly NDVI composites provides global coverage at a spatial resolution of 8 km × 8 km. For our analyses, we used the 672 composites for the 29 year period 1982–2010. We rescaled the data to the original floating-point NDVI scaling of −1 to 1, we identified and removed falsely included water pixels with a correction layer based on vector layers of continental coastline, lakes and reservoirs [39], and we applied a minimum threshold that set all pixels with NDVI values less than 0.05 to 0.05. NDVI values below 0.05 indicate non-vegetated areas like bare soil, water, snow, and ice. This 0.05 threshold is based on an analysis of stable pixels in the Sahara Desert that had an average NDVI of 0.094 and has been used previously in studies examining trends in NDVI [40].

3.2.Trend Analysis & Anthropogenic Effects

We estimated linear trends of annual means of NDVI for the 29 point time series from 1982 to 2010 for each spatial location in the data set (n = 1,975,140 pixels) using the Theil–Sen (TS) median estimator (function mblm in package mblm in R). The TS estimator represents a robust linear regression method that calculates the median slope among all pairs of observations. Commonly used for studies in vegetation trends based on NDVI time-series [20,41,42], the TS estimator is relatively insensitive to outliers and is accurate even for skewed and heteroscedastic data.

To relate these estimates to anthropogenic effects we performed three different analyses. First, we evaluated trends in NDVI for each of the six groups of anthropogenic biomes identified by Ellis and Ramankutti [13]; wildlands, rangelands, forests, croplands, villages, and dense settlements. We calculated histograms and medians for the TS estimators for each group of anthropogenic biomes and the global trend for each group based on the yearly average of NDVI for all pixels of a group using generalized least squares (GLS). We incorporated the temporal correlation structure with a 1st order autoregressive moving average. For all GLS analyses, we calculated the coefficient of determination based on the improvement of the log-likelihood from the fitted null model (intercept only) to the full fitted model (function r.squaredLR in R package MuMIn, [43]).

Our second and third analyses were carried out on the scale of ecoregions. We used the World Wildlife Fund’s classification scheme, which distinguishes ecoregions [19] that are distinct with respect to natural communities and resident species. We first tested whether human-population density, as a surrogate for the degree of anthropogenic effects, was related to NDVI trends. We estimated human population density for each ecoregion as the mean of the 1990 and 2010 population estimates [44]. We next tested for the effects of land conversion, nitrogen deposition, and irrigation on NDVI trends by performing multiple regressions that included these factors as predictors. As a proxy for converted lands, we used the percent area covered by dense settlements, villages, and croplands per ecoregion [13]. For nitrogen deposition we used the mean nitrogen deposition per ecoregion in 1993 ([45], an average from the 1980s to the 2010s would have been preferable but was not available). For irrigation we calculated the percentage of irrigated area per ecoregion based on the UN-FAO irrigation map [46].

For both analyses we used the ecoregional means of the per-pixel TS estimators as the response variable. We used GLS accounting for spatial autocorrelation with an exponential spatial correlation function and accounting for differences in sizes of ecoregions using the inverse of ecoregion size as a weight in the variance function. For the multiple regression, we tested predictors for collinearity and used a cutoff correlation coefficient of 0.7 [47].

For robustness, we first repeated the analyses excluding all pixels located in the moist tropical broadleaf biome, because of unresolved and ongoing debates about the validity of remote sensing estimates in moist tropical regions [20,48] and the poor agreement between (1) GIMMS3g and MODIS Terra NDVI and (2) MODIS Terra and Aqua NDVI data in the moist tropics found by Fensholt and Proud [20]. Second, to ensure that our results do not stem from humans choosing areas with greater NDVI or from areas of greater NDVI exhibiting greater absolute increases of NDVI, we repeated our analyses by standardizing the per-pixel TS estimators by the mean NDVI value of each pixel. Third, even though ecoregions represent ecologically distinct units and as such are appropriate geographic units for examining the effects of covariates on trends in NDVI, we also repeated our analyses substituting squared grid cells of 500 km by 500 km for ecoregions. Lastly, to examine whether our results hold true also within continents and not just on the global scale, we performed the analyses for population density separately for each continent. Because some continents have only relatively few ecoregions, we did not perform the by-continent analyses for the multiple regression analyses.

4.Conclusions

In conclusion, our study provides one of the first global synopses quantifying the effect of the Earth’s human footprint on trends in photosynthetic capacity. Globally, more than 20% of the variability in NDVI trends can be associated with human land-use practices. Consequently, global climate change models and analyses of primary productivity should incorporate direct anthropogenic effects. Land use practices in densely populated areas, and other areas featuring major human impacts through land conversion, irrigation, and nitrogen deposition, may contribute more to global changes in photosynthetic capacity than previously anticipated.

Acknowledgments

WFF and TM were supported by NSF ABI award 1062411 and TM was supported by the Robert Bosch Foundation. This work was also supported by the National Socio-Environmental Synthesis Center (SESYNC) under funding received from the National Science Foundation DBI-1052875.

Conflicts of Interest

The authors declare no conflict of interest.

Author Contributions

Thomas Mueller, Gunnar Dressler, Compton J. Tucker and William F. Fagan conceptualized the study, Thomas Mueller and Gunnar Dressler analyzed the data, and all authors contributed to the writing of the manuscript.

References1.AnyambaA.C.TuckerJ.Analysis of Sahelian vegetation dynamics using NOAA-AVHRR NDVI data from 1981–2003J. Arid Environ2007635966142.JeyaseelanT.RoyP.S.YoungS.S.Persistent changes in NDVI between 1982 and 2003 over India using AVHRR GIMMS (Global Inventory Modeling and Mapping Studies) dataInt. J. Remote Sens2007284927494610.1080/014311607012532793.OlssonL.EklundhL.ArdöJ.A recent greening of the Sahel—Trends, patterns and potential causesJ. Arid Environ20056355656610.1016/j.jaridenv.2005.03.0084.MyneniR.B.NemaniR.R.RunningS.W.Estimation of global leaf area index and absorbed par using radiative transfer modelsIEEE Trans. Geosci. Remote Sens1997351380139310.1109/36.6497885.NemaniR.R.KeelingC.D.HashimotoH.JollyW.M.PiperS.C.TuckerC.J.MyneniR.B.RunningS.W.Climate-driven increases in global terrestrial net primary production from 1982 to 1999Science20033001560156310.1126/science.1082750127919906.TuckerC.J.SellersP.J.Satellite remote-sensing of primary productionInt. J. Remote Sens198671395141610.1080/014311686089489447.HicklerT.EklundhL.SeaquistJ.W.SmithB.ArdoJ.OlssonL.SykesM.T.SjostromM.Precipitation controls Sahel greening trendGeophys. Res. Lett20053210.1029/2005GL0243708.CaoM.K.PrinceS.D.SmallJ.GoetzS.J.Remotely sensed interannual variations and trends in terrestrial net primary productivity 1981–2000Ecosystems200472332429.SandersonE.W.JaitehM.LevyM.A.RedfordK.H.WanneboA.V.WoolmerG.The human footprint and the last of the wildBioscience20025289190410.1641/0006-3568(2002)052[0891:THFATL]2.0.CO;210.MittermeierR.A.MittermeierC.G.BrooksT.M.PilgrimJ.D.KonstantW.R.da FonsecaG.A.B.KormosC.Wilderness and biodiversity conservationProc. Natl. Acad. Sci. USA2003100103091031310.1073/pnas.17324581001293089811.FoleyJ.A.DeFriesR.AsnerG.P.BarfordC.BonanG.CarpenterS.R.ChapinF.S.CoeM.T.DailyG.C.GibbsH.K.Global consequences of land useScience200530957057410.1126/science.11117721604069812.HurttG.C.ChiniL.P.FrolkingS.BettsR.A.FeddemaJ.FischerG.FiskJ.P.HibbardK.HoughtonR.A.JanetosA.Harmonization of land-use scenarios for the period 1500–2100: 600 years of global gridded annual land-use transitions, wood harvest, and resulting secondary landsClim. Chang201110911716110.1007/s10584-011-0153-213.EllisE.C.RamankuttyN.Putting people in the map: Anthropogenic biomes of the worldFront. Ecol. Environ2008643944710.1890/07006214.BaldiG.NosettoM.D.AragonR.AversaF.ParueloJ.M.JobbagyE.G.Long-term satellite NDVI data sets: Evaluating their ability to detect ecosystem functional changes in South AmericaSensors200885397542510.3390/s809539715.BradleyB.A.MustardJ.F.Comparison of phenology trends by land cover class: A case study in the Great Basin, USAGlob. Chang. Biol20081433434616.NeighC.S.R.TuckerC.J.TownshendJ.R.G.North American vegetation dynamics observed with multi-resolution satellite dataRemote Sens. Environ20081121749177210.1016/j.rse.2007.08.01817.EastmanJ.R.SangermanoF.MachadoE.A.RoganJ.AnyambaA.Global trends in seasonality of Normalized Difference Vegetation Index (NDVI), 1982–2011Remote Sens201354799481810.3390/rs510479918.NiyogiD.KishtawalC.TripathiS.GovindarajuR.S.Observational evidence that agricultural intensification and land use change may be reducing the Indian summer monsoon rainfallWater Resour. Res201010.1029/2008WR00708219.OlsonD.M.DinersteinE.WikramanayakeE.D.BurgessN.D.PowellG.V.N.UnderwoodE.C.D’AmicoJ.A.ItouaI.StrandH.E.MorrisonJ.C.Terrestrial ecoregions of the worlds: A new map of life on earthBioscience20015193393810.1641/0006-3568(2001)051[0933:TEOTWA]2.0.CO;220.FensholtR.ProudS.R.Evaluation of earth observation based global long term vegetation trends—Comparing GIMMS and MODIS global NDVI time seriesRemote Sens. Environ201211913114710.1016/j.rse.2011.12.01521.De JongR.de BruinS.de WitA.SchaepmanM.E.DentD.L.Analysis of monotonic greening and browning trends from global NDVI time-seriesRemote Sens. Environ201111569270210.1016/j.rse.2010.10.01122.DonohueR.J.McVicarT.R.RoderickM.L.Climate-related trends in Australian vegetation cover as inferred from satellite observations, 1981–2006Glob. Chang. Biol2009151025103910.1111/j.1365-2486.2008.01746.x23.VerbylaD.Browning boreal forests of western North AmericaEnviron. Res. Lett2011610.1088/1748-9326/6/4/04100324.BeckP.S.A.JudayG.P.AlixC.BarberV.A.WinslowS.E.SousaE.E.HeiserP.HerrigesJ.D.GoetzS.J.Changes in forest productivity across Alaska consistent with biome shiftEcol. Lett20111437337910.1111/j.1461-0248.2011.01598.x2133290125.JulienY.SobrinoJ.A.VerhoefW.Changes in land surface temperatures and NDVI values over Europe between 1982 and 1999Remote Sens. Environ2006103435510.1016/j.rse.2006.03.01126.StockliR.VidaleP.L.European plant phenology and climate as seen in a 20-year AVHRR land-surface parameter datasetInt. J. Remote Sens2004253303333010.1080/0143116031000161814927.PiaoS.L.FangJ.Y.ZhouL.M.GuoQ.H.HendersonM.JiW.LiY.TaoS.Interannual variations of monthly and seasonal Normalized Difference Vegetation index (NDVI) in China from 1982 to 1999J. Geophys. Res.: Atmos200310.1029/2002JD00284828.HeumannB.W.SeaquistJ.W.EklundhL.JonssonP.AVHRR derived phenological change in the Sahel and Soudan, Africa, 1982–2005Remote Sens. Environ200710838539210.1016/j.rse.2006.11.02529.BegueA.VintrouE.RuellandD.CladenM.DessayN.Can a 25-year trend in Soudano-Sahelian vegetation dynamics be interpreted in terms of land use change? A remote sensing approachGlob. Environ. Chang20112141342010.1016/j.gloenvcha.2011.02.00230.De VriesW.PoschM.Modelling the impact of nitrogen deposition, climate change and nutrient limitations on tree carbon sequestration in Europe for the period 1900–2050Environ. Pollut20111592289229910.1016/j.envpol.2010.11.0232116356131.The State of the World’s ForestsFood and Agriculture Organization of the United NationsRome, Italy201132.ClarkM.L.AideT.M.GrauH.R.RinerG.A scalable approach to mapping annual land cover at 250 m using MODIS time series data: A case study in the Dry Chaco ecoregion of South AmericaRemote Sens. Environ20101142816283210.1016/j.rse.2010.07.00133.PinzonJ.E.TuckerC.J.A non-stationary 1981–2012 AVHRR NDVI3g time seriesRemote Sens2014submitted34.TuckerC.J.PinzonJ.E.BrownM.E.SlaybackD.A.PakE.W.MahoneyR.VermoteE.F.El SaleousN.An extended AVHRR 8-km NDVI dataset compatible with MODIS and SPOT vegetation NDVI dataInt. J. Remote Sens2005264485449810.1080/0143116050016868635.ZhouL.M.TuckerC.J.KaufmannR.K.SlaybackD.ShabanovN.V.MyneniR.B.Variations in northern vegetation activity inferred from satellite data of vegetation index during 1981 to 1999J. Geophys. Res.: Atmos2001106200692008310.1029/2000JD00011536.SobrinoJ.A.JulienY.Global trends in NDVI-derived parameters obtained from GIMMS dataInt. J. Remote Sens2011324267427910.1080/01431161.2010.48641437.BeckH.E.McVicarT.R.van DijkA.SchellekensJ.de JeuR.A.M.BruijnzeelL.A.Global evaluation of four AVHRR-NDVI data sets: Intercomparison and assessment against Landsat imageryRemote Sens. Environ20111152547256310.1016/j.rse.2011.05.01238.JusticeC.O.TownshendJ.R.G.VermoteE.F.MasuokaE.WolfeR.E.SaleousN.RoyD.P.MorisetteJ.T.An overview of MODIS Land data processing and product statusRemote Sens. Environ20028331510.1016/S0034-4257(02)00084-639.Natural Earth Data SetsAvailable online: http://www.naturalearthdata.com/(accessed on 30 May 2014)40.SlaybackD.A.PinzonJ.E.LosS.O.TuckerC.J.Northern hemisphere photosynthetic trends 1982–99Glob. Chang. Biol2003911510.1046/j.1365-2486.2003.00507.x41.De BeursK.M.HenebryG.M.A statistical framework for the analysis of long image time seriesInt. J. Remote Sens2005261551157310.1080/0143116051233132665742.Alcaraz-SeguraD.ChuviecoE.EpsteinH.E.KasischkeE.S.TrishchenkoA.Debating the greening vs. Browning of the North American boreal forest: Differences between satellite datasetsGlob. Chang. Biol20101676077010.1111/j.1365-2486.2009.01956.x43.NagelkerkeN.J.D.A note on a general definition of the coefficient of determinationBiometrika19917869169210.1093/biomet/78.3.69144.Socioeconomic Data and Applications Center, Center for International Earth Science Information Network (CIESIN), Columbia University and Centro Internacional de Agricultura TropicalGridded Population of the World, Version 3, SEDAC2005Available online: http://sedac.ciesin.columbia.edu/gpw(accessed on 30 May 2014)45.Global Maps of Atmospheric Nitrogen Deposition, 1860, 1993 and 2050Available online: http://daac.ornl.gov//CLIMATE/guides/global_N_deposition_maps.html(accessed on 30 May 2014)46.Global Map of Irrigation AreasAvailable online: http://www.fao.org/nr/water/aquastat/irrigationmap/index10.stm(accessed on 30 May 2014)47.DormannC.F.ElithJ.BacherS.BuchmannC.CarlG.CarréG.MarquézJ.R.G.GruberB.LafourcadeB.LeitãoP.J.Collinearity: A review of methods to deal with it and a simulation study evaluating their performanceEcography201236274648.SaleskaS.R.DidanK.HueteA.R.da RochaH.R.Amazon forests green-up during 2005 droughtScience200710.1126/science.1146663Figures and TablesFigure 1.

Annual change of NDVI of ecoregions, based on ecoregional means of per-pixel Theil-Sen estimators) vs. human population density. Circle diameters are proportional to the size (area) of each ecoregion. Red line indicates trend based on generalized least squares model with exponential spatial autocorrelation structure and accounting for differences in sizes of ecoregions in the variance function (β: 0.00032, p < 0.0001, Δ AIC: 105, coef. determination: 0.13). For corresponding analyses by continent, see Figure 3.

Figure 3.

Annual change of NDVI of ecoregions (Africa (A), Asia (B), Australia (C), Europe (D), North America (E), South America (F), based on ecoregional means of per-pixel Theil-Sen estimators) vs. log10 of human population density. Circle diameters are proportional to the size (area) of each ecoregion. Red line indicates trend based on generalized least squares model with exponential spatial autocorrelation structure and accounting for differences in sizes of ecoregions in the variance function. β: coefficient of log10 of human population density, d: coefficient of determination, significance codes: *** p < 0.0001, ns: not significant.

Figure 4.

(A) Trends in NDVI across the globe, 1981–2010. Ecoregional [19] extremes for NDVI increase (defined as the 5% of land surface with the fastest increases in NDVI, n = 73 ecoregions) are in red, whereas ecoregional extremes for NDVI decrease (defined as the 5% of land surface with the fastest decreases in NDVI, n = 38 ecoregions) are in blue; (B) Boxplots contrast the ecoregional extremes in A for increases (n = 73) and decreases (n = 38) in terms of NDVI trends, population density, percent converted lands [13], percent irrigated lands, and nitrogen deposition.

Additional models relating trends in NDVI to (A) population density and (B) to nitrogen deposition, irrigation, and converted lands in order to ensure robustness of results against potential methodological concerns (see main text).