Direct and indirect impacts of urbanization on vegetation growth across the world's cities.

Urban environments, regarded as "harbingers" of future global change, may exert positive or negative impacts on urban vegetation growth. Because of limited ground-based experiments, the responses of vegetation to urbanization and its associated controlling factors at the global scale remain poorly understood. Here, we use satellite observations from 2001 to 2018 to quantify direct and indirect impacts of urbanization on vegetation growth in 672 worldwide cities. After controlling for the negative direct impact of urbanization on vegetation growth, we find a widespread positive indirect effect that has been increasing over time. These indirect effects depend on urban development intensity, population density, and background climate, with more pronounced positive effects in cities with cold and arid environments. We further show that vegetation responses to urbanization are modulated by a cities' developmental status. Our findings have important implications for understanding urbanization-induced impacts on vegetation and future sustainable urban development.

INTRODUCTION

Human society has entered an era of increasing urbanization. More than half of the world’s population now lives in cities, and approximately 68% (i.e., 6.7 billion people) will live in urban areas by 2050 (). Urbanization and human activities radically modify landscapes and their ecology (). Rapid urban development has led to widespread conversion of vegetated areas to impervious surfaces, profoundly changing the atmospheric and climatic conditions of urban areas [e.g., urban heat island (UHI) effect, increased CO2 concentrations, and air pollution] (, ). Urban environments are thus often considered as “harbingers” of future global change (, ), exert multiple effects on vegetation growth, and serve as “natural laboratories” for vegetation growth studies (, ). As a result, an increasing number of studies has focused on the ecology of populated landscapes over the past decades (, , ). Identifying the effects of urban environments on vegetation growth across the globe along with the mechanisms driving it can help improve our understanding of vegetation growth responses to global environmental change (, ).

The impacts of urbanization on vegetation growth can be decomposed into direct and indirect effects (). The direct effect is generally negative and refers to the transformation of land cover from natural surfaces to impervious ones, reducing vegetation cover and growth (, ). The indirect urbanization effect on vegetation growth is usually caused by human management practices in cities and higher air temperature compared to surrounding natural areas (). This indirect impact on vegetation has been reported in previous studies (–), but the findings often differed in direction and extent. While some early horticultural studies showed that urban stressors have a negative impact on tree growth (, ), recent studies, based on manipulative experiments (, ) or regional remote sensing observations (–), suggested that vegetation growth in cities is stimulated by the altered urban environment (e.g., more productive or longer growing seasons). However, the studied cities largely differed in their background climate and urban development levels, which may explain the ambiguous and controversial findings. The different observation scopes and methods used in previous studies make it difficult to infer general patterns.

Quantifying the indirect impact of urban environments on vegetation growth is challenging. Satellite-based remote sensing has been widely used to monitor vegetation growth, productivity, and phenology by using vegetation indices (VIs) such as the normalized difference VI (NDVI) or the enhanced VI (EVI). Remote sensing is therefore a powerful tool to characterize spatiotemporal changes of land surfaces, including urban vegetation. Zhao et al. () proposed a conceptual framework to explicitly quantify the direct (ωd) and indirect (ωi) urbanization effects on vegetation growth by characterizing the relationship between urbanization intensity (UI) (β; represented by the proportion of impervious surface area) and satellite-based VIs (Fig. 1). This approach makes it possible to compare the theoretical linear decrease of VI with increasing β along urban-rural gradients (i.e., the zero-impact line) with the observed nonlinear decrease (i.e., the fitted VI~β curve). The indirect impact can then be quantified by calculating the difference between the fitted VI~β curve and the zero-impact line (see Materials and Methods for details). This method for quantifying the indirect impact of urbanization on vegetation growth was first applied to 32 major cities in China (), and several studies have since confirmed the applicability of the methodology (, , ). All studies showed a prevalent enhancement of vegetation growth in urban environments but were restricted to the national or regional scale in temperate zones. This severely limits generalizations at the global scale, and the spatiotemporal patterns of urbanization-induced indirect effects on vegetation and their dependence on climatic background conditions have yet to be evaluated. Furthermore, the factors and mechanisms driving vegetation responses to urban environments across spatial gradients are still poorly understood.

Fig. 1.

The conceptual framework for quantifying the direct and indirect effects of urbanization on vegetation based on satellite-derived observations.

The color changes from green to gray represent the gradient from rural to urban areas with a reduction of vegetation cover and an increase in impervious surface cover. The satellite-derived low-resolution blocks of VI and high-resolution blocks of impervious surface area (ISA) are overlaid to characterize the relationship between UI (β) and VI. The blue points and the solid blue line represent the observed VI (Vobs) values and their regression line, respectively. Vv and Vnv are the VI values of fully vegetated rural surfaces and unvegetated urban surfaces, respectively. The dashed blue line represents the VI change along the urbanization gradient assuming no indirect impact (Vzi; the zero-impact line). A positive or negative indirect impact of urbanization is then inferred if the VI values fall above or below the zero-impact line, respectively. The computational framework reflected in this figure is derived from the work of Zhao et al. ().

Here, we provide a global analysis of the impact of urban environments on vegetation growth from 2001 to 2018 across the world’s 672 largest urban areas, including tropical (99 cities), temperate (286 cities), cold (187 cities), and arid environments (100 cities) (fig. S1). We used EVI as an indicator of vegetation productivity (), extracted from the Moderate Resolution Imaging Spectroradiometer (MODIS) data product. The impervious surface data from high-resolution Landsat imagery were used to calculate UI (see Materials and Methods for details about data sources and data processing). Continuous relationship curves between VI and β were generated for each city. We then examined geospatial and temporal patterns of vegetation growth responses to urbanization across global cities within different climate zones to answer the following questions: Does the indirect impact of urbanization on vegetation depend on regional climatic conditions? How is this indirect impact driven by the changing environmental and anthropogenic factors in urban areas and how does it vary across climate zones? And does the indirect effect of urbanization on vegetation growth and its drivers depend on the levels of development in cities? Answering these questions is important for understanding urban vegetation growth dynamics and their underlying mechanisms, thereby providing information that can aid to the sustainable development and management of cities.

RESULTS

VI~β relationship and ωi across world cities

The relationship between the VI and UI (β) was fitted for each city. The VI~β curves and their corresponding zero-impact lines [the straight line connecting the VI observed for pixels with 100% vegetation (β = 0) to the VI observed for pixels with 100% impervious surface (β = 1)] for 12 representative cities are shown in fig. S2. Across all cities, VI decreased with increasing β. This demonstrates the negative direct effect of urbanization, whereby the coverage of vegetated surface decreases with increasing surface sealing (fig. S3). However, in most cases, the actual change of observed VI values along β did not fully match the straight zero-impact line, which can be attributed to indirect effects (ωi). The observed VIs were commonly higher than the zero-impact line predictions (fig. S2), indicating a widespread positive response of vegetation growth to urban environments (Fig. 2A).

Fig. 2.

Spatial and temporal patterns of the urbanization-induced indirect impact (ωi) on vegetation.

(A and C) The spatial distribution of the magnitude (A) and the trend (C) of ωi in the world’s major cities are overlaid on the map of climate zones. The temporal trends and statistical significance were estimated using the Mann-Kendall test. Points with green or brown color represent cities in which ωi significantly (P < 0.05) increased or declined over time, respectively. Points with no filling had no significant trends. (B and D) Changes with latitude in the magnitude (B) and the trend (D) of ωi. Solid lines and shaded areas represent the mean and SD. (E and F) Averaged values of ωi (E) and ωi trends (F) for varying climatic backgrounds (mean annual air temperature and precipitation).

Overall, 655 (97%) of the 672 global cities exhibited an indirect vegetation enhancement due to urbanization (ωi > 0), with a global average enhancement of ~26%. The cities with the largest vegetation enhancement were located in the east of China. The few cities with negative ωi were located in South America and Asia. The magnitude of the indirect impact of urbanization was larger in the Northern Hemisphere than in the Southern Hemisphere (Fig. 2B). In both the Northern and the Southern Hemispheres, ωi increased from low to high latitudes, with areas near the equator showing the lowest average indirect impact.

The trend of the indirect impact from 2001 to 2018 is shown in Fig. 2C, with 43% (288 of 672) of cities showing a significant increase of ωi (slope > 0, P < 0.05). These cities are mainly distributed in mid-latitude areas in the Northern Hemisphere. Only 3% (22 of 672) of cities experienced significant declines in ωi (slope < 0, P < 0.05), and these are mainly located in low-latitude areas. The latitudinal pattern of the temporal trend in ωi was similar to that of the average ωi (Fig. 2D). Detailed quantitative results are shown in table S1.

Climatic and anthropogenic factors driving the indirect urbanization effect on vegetation growth

The latitudinal pattern of the magnitude and the trend of the indirect effect on vegetation growth in urban environments suggests that the spatial distribution of ωi might be determined by changes in temperature. Our analysis shows that ωi depends on a city’s climate background (Fig. 2, A to C), whereby the magnitude of ωi and its temporal trend increased toward dry and cold regions, i.e., areas with low annual precipitation and high mean annual air temperature (Fig. 2, E and F). This indicates that vegetation growth exhibits a stronger positive indirect response to urban environments when water and thermal supply are limited.

Figure 3A shows the relationship between VI and β aggregating the city-level estimates within each climate zone (tropical, temperate, cold, and arid regions). The quantified ωi along β in these regions is shown in Fig. 3B. The most pronounced indirect vegetation enhancement occurred in arid regions (average ωi of ~39%), followed by cold (~32%) and temperate regions (~22%). The lowest vegetation enhancement with an average ωi below 20% was found for tropical regions. Furthermore, while the average direct effect of urbanization did not change between the 2001 to 2009 and 2010 to 2018 periods (except for slight increases in arid regions; Fig. 3C), the positive indirect effect was higher from 2010 to 2018 than from 2001 to 2009 (Fig. 3D), indicating an increasing urban vegetation enhancement over time. This increase was most pronounced in cold and temperate regions, in which the majority of the world’s largest cities are located.

Fig. 3.

Relationships between VI and urban intensity (β) under different climatic backgrounds.

(A) Relationships between the normalized EVI and β, globally and across four climate zones. Points are the averaged VI values within each β bin with an interval of 0.02 (n = 50). The solid and dashed lines are cubic regression lines and zero-impact lines. Data from 2001 to 2009 and from 2010 to 2018 are represented in blue and orange colors, respectively. (B) Relationships between the indirect impact of urbanization on vegetation (ωi) and β, globally and across four climate zones. Points are the averaged ωi values within each β bin. The blue, orange, and gray dashed lines are means of ωi from 2001 to 2009, means of ωi from 2010 to 2018, and zero lines, respectively. Mean values of direct impact (ωd) (C) and ωi (D) of cities in two time periods, globally and across four climate zones. Error bars represent SEM.

In addition to the climatic background, anthropogenic factors may also play an important role in affecting the response of vegetation to urban environments. We further analyzed relationships of ωi with the UHI effect (∆T), urban greenness (UG), averaged UI, and population density (POP) (fig. S4) across the globe and within climate zones by using partial correlation analysis (Fig. 4, A to E). Globally, ωi was negatively correlated (P < 0.05) with air temperature (Tair) and ∆T and positively correlated (P < 0.01) with UI and POP. However, our analyses also show that the effect magnitudes and directions differ across climate zones. In tropical regions, mean annual temperature and precipitation barely affected ωi, and, instead, UI was most correlated with ωi in this region. In the temperate region, the extent of the indirect impact on urban vegetation was dominated by anthropogenic factors, with UI and POP both exhibiting a positive effect on ωi. The insensitivity of ωi to mean annual temperature and precipitation may be attributed to the fact that human management alleviates the limiting effects of temperature (extreme heat or cold conditions) and water availability on vegetation. ωi was negatively correlated with ∆T in temperate regions, contradicting the positive response of vegetation to the UHI effect observed in two previous studies (, ). Instead, these results support a negative response of vegetation growth to increased temperature in warm areas (), which might be partly explained by declining phenological sensitivity in these regions (, ). In cold regions, mean annual temperature and precipitation within cities had a positive effect on ωi, most likely because the limiting effect of low temperature and water supply on vegetation is more pronounced in cold environments (). The negative relationship between UG and ωi in cold regions suggests that positive indirect effects of urbanization on vegetation growth are most pronounced in cities with only little vegetation (, ). In arid regions, cities with low temperature and green coverage as well as high UI and POP showed the strongest ωi increase, indicating that positive indirect effects of urbanization on vegetation growth are particularly pronounced in cold cities with a high degree of urban management (). Overall, these analyses show that climatic and anthropogenic factors jointly influence the urbanization-induced indirect impact on vegetation, with anthropogenic factors explaining a higher share of the variation in ωi than climatic factors across most regions.

Fig. 4.

Effects of climatic and anthropogenic factors on the indirect urbanization impact on vegetation (ωi).

(A to J) Partial correlations between ωi (A to E) or the temporal trend in ωi (F to J) and six explanatory variables in global (A and F), tropical (B and G), temperate (C and H), cold (D and I), and arid (E and J) cities. Tair, mean annual air temperature; PRE, mean annual precipitation; ∆T, urban-rural temperature difference; UG, UG as the mean EVI in the urban area; UI, mean of UI; POP, mean of urban POP. All variables with superscript represent the trend of those variables estimated using the Mann-Kendall test. *P < 0.05 and **P < 0.01.

To study the drivers of temporal trends in ωi, we calculated the trends of the six variables representing climatic and anthropogenic factors (Tair, PRE, ∆T, UG, UI, and POP) and ran partial correlation analyses. The variables affecting ωi trends were consistent across regions. Temporal trends in ωi were best explained by changes in UI (UI) with a significant positive relationship (P < 0.01) at the global scale and within climate zones (Fig. 4, F to J). However, increases in POP were associated with decreases in ωi in temperate and tropical regions, indicating that fast population growth hinders effective management of urban vegetation. In temperate and arid regions, increasing air temperature and a more pronounced UHI effect over time weakened the increasing trend of ωi, likely because excessive heat waves in those cities pose an increasing threat to vegetation growth. Conversely, in cold regions, cities with strong increases in mean annual temperature showed the strongest increases in ωi, confirming that urban vegetation in cold environments responds positively to climate warming (, ). In tropical and temperate regions, increases in POP within cities were negatively correlated with changes in ωi, probably because increases in population lead to deforestation () and greater investment in hard infrastructure, such as buildings and streets (), reducing the positive response of vegetation to urban environments.

Differences in indirect impact between cities at different development levels

While the indirect enhancement of vegetation growth in urban environments generally increased with UI, the effects of urbanization on vegetation growth tended to differ between cities with different levels of economic development. For example, the fitted VI~β curve was below the zero-impact line within some areas with low urban intensities in Beijing, while this was not evident in Paris (fig. S5). We further compared the magnitudes of the indirect urbanization impact on vegetation (ωi) in city areas with low urbanization (β from 0 to 0.5) and high urbanization (β from 0.5 to 1) between cities with different levels of development (Fig. 5A). The average indirect impact (~31%) in developed cities was higher than that (~24%) in developing cities, with particularly pronounced differences in low urbanization areas. Notably, in developing cities, the average indirect impact was near zero or even negative in low urbanization areas, while in the majority of developed cities, a positive indirect effect was observed in low urbanization areas. This difference may relate to the stage of urban development, with rapid urban expansion in developing cities exerting a negative indirect impact on the vegetation surrounding the city. By contrast, the average ωi in highly urbanized areas was similar for both types of cities, although the variance of ωi in developing cities was larger than that in developed cities. Moreover, within highly urbanized city areas, developing cities showed stronger temporal increases in ωi than developed cities (Fig. 5B). These results suggest that the management of urban vegetation in developing cities is concentrated in highly urbanized core areas, while the management of urban vegetation in developed cities (at a relatively mature urbanization stage) is more balanced between highly urbanized and suburban areas.

Fig. 5.

The indirect impact (ωi) of urbanization on vegetation varying with different levels of urban development (developing or developed cities) and different ranges of UI (β).

ωi (A) and ωi trend (B) in developed and developing cities within low (β from 0 to 0.5) and high (β from 0.5 to 1) β ranges. Box plots show the interquartile range (IQR) (box), the median (horizontal line in box), and 1.5× IQR (whiskers); outliers are omitted for clarity.

DISCUSSION

Limited availability of ground-based observations and incomplete integration of climatic and human-related factors have hampered our understanding of the patterns and causes of urbanization-induced impacts on vegetation growth. Here, we used remote sensing data with high spatial and temporal resolution to characterize the impact of urban environments on vegetation growth across the world’s major cities. By separating direct and indirect effects of urbanization on vegetation growth along urban-rural gradients (), we show that, across the majority of cities, positive indirect urbanization effects exist in parallel with negative direct effects. This indirect enhancement of vegetation growth might offset the direct loss of vegetation due to land transformation in urban areas to some extent and matches previous studies that found longer growing seasons in urban areas compared to their surrounding rural areas (, , ). We further show that, at a global scale, climatic and anthropogenic factors jointly control the indirect urbanization effect on vegetation. Positive indirect urbanization effects were most pronounced in highly urbanized cities located in cold and arid regions, indicating that urban vegetation in these areas is strongly affected by human management (e.g., irrigation) and thus is less limited by water, temperature, and/or nutrient supply than the natural vegetation outside urban areas (). This suggests that the climatic background of cities mediates the link between urbanization and vegetation growth and that city-specific human activities (adaptation to local climate conditions) act as an important driver of indirect urbanization effects on vegetation. The observed enhancement of vegetation greenness in cities may also have implications for the effects of future climate change on vegetation growth (, ), predicting increasing vegetation productivity under future environmental conditions. Yet, the complex feedbacks between the environment and global vegetation lead to great uncertainty in future projections (, ).

The positive indirect effect of urbanization on vegetation greenness may also be partly explained by changes in plant biomass allocation (, , ). High UI usually means more built-up areas and less open space, reducing the light available to vegetation (). To avoid light limitation, plants may thus increase their relative biomass allocation to shoots (, ). Yet, belowground resources might also become more limiting in urban environments. While human management of urban vegetation might, to some extent, alleviate belowground resource supply, for instance, through supplement of water and nutrients (, ), plant root health is often compromised by soil contamination, water shortage, and reductions in litterfall and decomposition because of surface sealing and soil compaction (, ). A healthy balance between the crown and root systems will be critical to sustain urban vegetation growth under more extreme climate conditions in the future.

The negative correlation between ωi and the UHI effect in temperate regions (Fig. 4C) indicates that urban heat does not necessarily enhance vegetation growth and might even hinder vegetation growth in warm biomes (, ). Global remote sensing observations have instead shown enhanced vegetation growth under rising temperatures (, ), pointing toward different temperature sensitivities of vegetation in natural and urban environments (, ). Enhanced urban vegetation growth may also increase the cooling effect of vegetation, mitigating excessive urban heat (, ). Furthermore, the insignificant correlation between ωi and the UHI effect in other climate zones indicates that the spatial variability in temperature caused by the UHI effect has only limited influence on urban vegetation growth outside temperate areas (, , ), underscoring the importance of other site-specific factors, such as air and light pollution, soil modifications, and biotic interactions, that merit further investigation ().

Our results also demonstrate differences in the indirect effects of urbanization between cities at different development levels (Fig. 5). The larger positive indirect urbanization effect in highly developed cities relative to developing cities might be explained by intensified management of vegetation in both city centers and rural areas. The large variation in the indirect urbanization impact within cities in developing countries is likely linked to their large variation in urban development levels and management practices (, ). The lower indirect impact in low urbanization areas in cities of developing countries compared with that in developed countries indicates that developed cities can invest more resources in the ecological governance in peri-urban areas, while increasing population pressure and excessive industrialization make it difficult for developing cities to balance economic development with ecological management. Yet, the temporal increases of ωi in highly urbanized areas in developing cities may reflect environmental improvements in the urban centers of cities characterized by recent, rapid urbanization.

Our results, based on satellite-derived observations of urban plant activity and artificial impervious area, show general patterns of vegetation responses to urban environments across global cities. A limitation of our approach is that the quantification of spatial variation in urban density, compactness, and sprawling shapes within cities () using impervious surface data is challenging. Different urban structures and development schemes may lead to complex variation in the indirect effects on vegetation. In addition, other environmental factors, such as tree species composition, air pollutants, and nitrogen deposition, may also affect vegetation growth in cities. To deepen our understanding of the urbanization effects on vegetation, it will thus be important to study the effects of other potential factors once high-resolution data on those factors are available at a global scale. Our findings are mostly in agreement with experimental studies on urban vegetation, highlighting the combination of long-term manipulative experiments with satellite observations as a promising future avenue for untangling the mechanisms underlying vegetation responses to urbanization.

Urban areas are hot spots that drive environmental changes, and global environmental changes are linked to changes in local urban environments (). Urban habitats represent an important but underused data source for studying the effects of environmental and anthropogenic factors on vegetation dynamics (). By providing insight into the responses of plant activity to climatic and anthropogenic stressors in urban areas, our quantitative results expand our understanding of the vegetation response to urbanization and the responses of plant growth to altered environmental conditions over space and time. Given the major challenges posed to cities by global climate change and population growth, knowledge of the uneven effects and diverse drivers of urbanization on vegetation can provide a scientific reference for improving urban governance and working toward sustainable and resilient cities of the future.

MATERIALS AND METHODS

Datasets

Urban boundary data

The urban areas in our study were extracted from the multitemporal global urban boundary (GUB) dataset (). Comparisons with other products about global urban areas suggested that delineated urban boundaries in GUB data can well capture the geometries of urban extent around urban fringe areas (). We selected the cities with areas larger than 100 km2, which ensured that each city contains a sufficient number of available remote sensing pixels for further calculations. To comprehensively cover changes in vegetation along the urban-rural gradient for each city, we created a buffer outward from each urban boundary. To ensure a balance of size between urban and rural areas, we determined the buffer distance as followswhere S is the area of a city. This function can determine the buffer distance by making the size of rural area approximately the same with the size of urban area under an assumption of a circular shape of city. This approach avoids the drawback of adopting a fixed threshold for cities with different sizes. We excluded the pixels representing water bodies or crops by using the MODIS yearly land cover product (MCD12Q1) with the classification system of the International Geosphere-Biosphere Programme (). Croplands were removed because they are strongly affected by human management. Moreover, to avoid the impact of elevation difference on vegetation growth, we excluded pixels in rural areas that had elevations 50 m higher or lower than the average elevation of urban pixels according to the GTOPO30 digital elevation model data ().

Satellite-sensed VI data

The EVI from the MOD13A1 version 6 product (16-day composite) with a spatial resolution of 500 m was used as the VI to characterize the spatial variability of vegetation growth in cities (). The EVI data have greater sensitivity in high-biomass regions and partially eliminate the effect of canopy background compared with NDVI data (). All MODIS EVI data from 2001 to 2018 were aggregated into six time periods with a 3-year interval (i.e., 2001–2003, 2004–2006, 2007–2009, 2010–2012, 2013–2015, and 2016–2018) to analyze the vegetation growth in cities. The pixels with low confidence (i.e., pixels covered with cloud, snow, or aerosol) were removed on the basis of the quality assessment layer in the MODIS product. The EVI values for each pixel over each time period were averaged for representing the average state of urban vegetation for the corresponding period.

Imperious surface area data for determining urban intensity

We adopted annual maps of global artificial impervious area (GAIA) () to quantify urban intensity (β). This dataset has a high-resolution of 30 m with long-term records from 1985 to 2018 using the full archive of Landsat images. The nighttime light data and the Sentinel-1 Synthetic Aperture Radar data were also used as the ancillary data for improving the accuracy of the dataset in arid areas. The overall accuracy of the GAIA dataset is higher than 90%, and the uncertainty of the data was lower in recent years (after 2000) because of the higher availability of Landsat data (). The urban intensity of each MODIS pixel was determined as the fraction of pixels falling into the impervious surfaces within the MODIS pixel (, ). Maps of urban intensity were generated for all six time periods, and values ranged from 0 (fully vegetated surfaces) to 1 (fully built-up surfaces).

Climate data

The global map of the Köppen-Geiger climate classification at a 1-km resolution was used to define the four climate zones (). We reclassified the original climate zones into the tropical, temperate, cold, and arid areas (fig. S1). To acquire the temperature and precipitation for cities, we used the atmospheric reanalysis product of the global climate based on the ERA5 monthly averaged data, which is produced by the European Centre for Medium-Range Weather Forecasts (). We extracted the mean annual air temperature at 2-m height and the annual precipitation of each city from this dataset, with a spatial resolution of 0.25°. We used the MOD11A2 V6 product (8-day composite) with a spatial resolution of 1 km to calculate the urban island intensity (), which was quantified by the difference of average land surface temperature (LST) between urban and rural areas ().

Data about urban development

We used three datasets to describe the urban development of population, economy, and UG. The Gridded Population of World, version 4 (GPWv4) product (obtained from the Socioeconomic Data and Applications Center) with a spatial resolution of 30 arc sec (approximately 1 km) was used to capture the human POP of cities (). The POP values of all pixels within each city were aggregated and averaged in each time period.

The economic development level of a country in which a city is located is also an important indicator of urban development. We collected the income levels of countries including high income, upper middle income, lower middle income, and low income as classified by the World Bank dataset (fig. S6) (). We reclassified the cities into developed (>$12,375; in countries that are high income) and developing (<$12,375; in countries that are not high income) types according to the gross national income per capita in 2018.

In addition, the UG is also an important indicator reflecting the sustainable management of vegetation in urban areas. The green spaces of a city such as urban parks or green roofs are usually developed with city expansion and improved urban environments (, ). To avoid interference from vegetation phenology across different regions around the world, we used the maximum greenness of vegetation (maximum EVI value at each pixel) in each time period to represent the best state of UG (). Thus, UG was determined as the mean of all maximum EVI values in an urban area. The details on the aforementioned datasets can be found in table S2.

Calculation of the direct and indirect impacts of urbanization on vegetation growth

A conceptual framework proposed in () was adopted to quantify the impacts of urbanization on vegetation growth, which can be separated into direct and indirect ones (Fig. 1). The direct impact is caused by land cover change due to the increase in impervious surface coverage in a city. Generally, the VI (represented by EVI) decreases with the gradient of urban intensity β (represented by the proportion of impervious surfaces within a VI pixel) from rural to urban area. Therefore, if only the direct impact is considered, then a linear relationship between VI and β would be present. This linear relationship was also defined as a theoretical zero-impact straight line without the indirect impactwhere Vzi was the theoretical VI of a 500-m resolution pixel, Vv was the mean of VI at the pixels completely filled by vegetated surfaces (β = 0), and Vnv was the minimum of VI at the pixels completely filled by built-up surfaces (β = 1). Thus, the relative direct urbanization impact (ωd) on vegetation growth can be calculated as

When the actual observed VI (Vobs) values against β of all pixels along the gradient of urban intensity are plotted, the distribution of points may not be completely consistent with the zero-impact straight line, suggesting the existence of indirect impact. This indirect impact means that vegetation growth may be enhanced or abated due to the difference between urban and natural environments. The points (observed VI values) above the zero-impact line indicate a positive indirect impact of urbanization on vegetation growth, and the points below the straight line indicate the existence of a negative impact. Therefore, the relative indirect impact (ωi) of urban environments on vegetation growth can be measured by the relative change of the observed VI values from the zero-impact VI line

To characterize the VI~β relationship, we derived the mean value of VI within each urban intensity β bin with an interval of 0.02 for each city. We used a cubic polynomial model to fit the VI~β curve: y = a0 + a1x + a2x2 + a3x3, where y was the observed EVI and x was β. This model was proved to be able to capture the VI~β relationship (, ). The coefficient of determination (R2) values of fitted functions were all larger than 0.6, and most (77%) were larger than 0.8 (fig. S7). We used the mean of all ωd and ωi values along urban intensities as the average direct and indirect impact values for each city. In this study, we focused on quantifying the indirect impact derived from complex altered urban environments and analyzed the spatial and temporal patterns of the indirect impact across the world’s cities.

The Mann-Kendall method (, ) was used to examine the trend of ωi from 2001 to 2018 for each city, as this method does not assume a specific distribution of the input data, and it is a nonparametric test for monotonic trends. The Theil-Sen method () was used to estimate the slope of the Mann-Kendall trend.

Analysis of controlling factors of the indirect impact on urban vegetation

The impact of urban environments on vegetation growth and its changing trends over time are probably controlled by the climatic and anthropogenic factors of a city. The spatial distributions of ωi and the ωi trend based on latitude, mean annual air temperature, and precipitation were analyzed (Fig. 2). We also aggregated all pixel data in each climate zone to determine the VI~β relationships of cities in different climate zones (Fig. 3). To investigate possible effects of climatic conditions and anthropogenic factors on ωi, we selected six variables as the potential driving factors, including mean annual air temperature (Tair), annual precipitation (PRE), UHI intensity (∆T; expressed by the urban-rural temperature difference), UG, averaged UI, and mean of POP. The mean of all pixel values for each potential driving factor within each city was calculated for the following analysis. Then, we used the partial linear regression model to calculate the partial correlations of ωi with these variables. The partial correlations of ωi trend with trends of air temperature (Tair) and precipitation (PRE) and trends of anthropogenic factors (∆T, UG, UI, and POP) were also determined. The partial correlation (two tailed) of each factor was calculated while controlling other factors. To remove the extreme cases from results of partial correlation analysis, we ran bootstrapped partial linear regression models by randomly subsampling 33% samples with replacement out of the entire data 500 times. For each subsample, a partial regression model was fitted for each covariate, and the corresponding slope (coefficient) was calculated. Then, we aggregated the results and calculated the mean coefficient across all model runs.

To test whether the selected six variables can be effectively used as the predictors of ωi, we used the random forest (RF) model to build the relationship between ωi and the covariates and adopted fivefold cross-validation to evaluate the performance of the fitted RF model. The mean R2 values (R2cv) of the fivefolds based on the data in cities at the global scale and in the four climate zones were calculated, respectively. All results showed an acceptable fitting effect (fig. S8), indicating that the selected six variables could be represented as the potential drivers of the magnitude of urbanization-induced indirect impact on vegetation. Furthermore, we used the fitted RF models to quantify the relative importance of each of the six variables. The importance metrics were determined by calculating the mean decrease in Gini (MDG), which is a measure of how each variable contributes to the homogeneity of the nodes and leaves in the resulting RF. The higher value of MDG means the higher improvement in the model performance if a variable was selected into the tree-building process in RF and the higher importance of this variable in the model. From the results, the obtained variable importance of the six variables was generally consistent with that of the partial correlation results (fig. S9), such as the general higher importance of UI and POP compared to other factors, which can also be found in Fig. 4.

Moreover, we distinguished the different patterns of VI~β and ωi~β relationships between cities in developed and developing countries. The differences of ωi in low (β from 0 to 0.5) and high (β from 0.5 to 1) urbanized areas of these two types of cities were analyzed.