Quantifying Urban Vegetation Dynamics from a Process Perspective Using Temporally Dense Landsat Imagery

Abstract

Urban vegetation can be highly dynamic due to the complexity of different anthropogenic drivers. Quantifying such dynamics is crucially important as it is a prerequisite to understanding its social and ecological consequences. Previous studies have mostly focused on the urban vegetation dynamics through monotonic trends analysis in certain intervals, but not considered the process which provides important insights for urban vegetation management. Here, we developed an approach that integrates trends with dynamic analysis to measure the vegetation dynamics from the process perspective based on the time-series Landsat imagery and applied it in Shenzhen, a coastal megacity in southern China, as an example. Our results indicated that Shenzhen was turning green from 2000–2020, even though a large-scale urban expansion occurred during this period. Approximately half of the city (49.5%) showed consistent trends in greening, most of which were located in the areas within the ecological protection baseline. We also found that 35.3% of the Shenzhen city experienced at least a one-time change in urban greenness that was mostly caused by changes in land cover types (e.g., vegetation to developed land). Interestingly, 61.5% of these lands showed trends in greening in the recent change period and most of them were distributed in build-up areas. Our approach that integrates trends analysis and dynamic process reveals information that cannot be discovered by monotonic trends analysis alone, and such information can provide insights for urban vegetation planning and management.

1. Introduction

Urban vegetation provides valuable ecosystem services, such as heat mitigation, air purification, and habitat preservation [1,2,3]. It also provides plenty of potential benefits for human health and thereby becomes of crucial importance for urban human wellbeing. For example, greenness exposure offers opportunities to reduce the risk of cardiovascular disease, metabolic disorder, and mental health [4,5,6]. Consequently, understanding vegetation dynamics through time is widely of concern for stakeholders of urban design and management.

Because of the advantages of temporally continuous and spatially explicit observations, satellite remote sensing data has been widely used and will be continually used to monitor the vegetation dynamics via measures of vegetation index [7,8,9,10]. Based on temporally dense data, such as AVHRR GIMMS, Terra MODIS, and SPOT VGT, the time-series vegetation indices, Normalized Difference Vegetation Index (NDVI) or Enhanced Vegetation Index (EVI), are largely used [11,12,13,14]. These studies mostly applied monotonic trends analysis that uses either the linear regression model or the Mann-Kendall test to investigate the trends of changes. An increase in NDVI/EVI time-series (generally defined as greening) indicates enhanced vegetation growth and a decreasing one (defined as browning) suggests reduced vegetation growth [15]. By integrating meteorological, environmental, and socio-economic data, the driving forces such as climate change, insect damage, and disturbance of fire, hurricane, and human activities can be inferred [16,17,18]. Such monotonic trends analysis has been successfully and extensively used to monitor the vegetation dynamics and evaluate the vegetation growth status in forested land and grassland [19,20].

However, the application of monotonic trend analysis for the vegetation dynamics in urban areas still remains challenging for reasons that are caused by the highly dynamic urban landscapes due to intensively anthropogenic activities [21,22]. First, the cover change with vegetation skewed the trends quantified by the monotonic trends analysis in a certain time interval [23,24,25]. Such total trends represented by monotonic greening or browning, without considering the dynamics, cannot accurately reflect the urban vegetation greenness. Second, the highly dynamic urban landscapes resulted in multiple processes in greenness, which cannot be reflected by monotonic trends analysis [26,27]. Numerous studies have demonstrated that the vegetation browning emerged when urban outward expansion converted the forested land, grassland, and cultivated land to the developed land [28,29]. Nevertheless, the complex socio-ecological activities, such as tree planting in new urban areas and lawns recreating in downtown areas, also led to vegetation recovery and thereby contributed to the vegetation greening [29,30,31,32,33]. In the context of massive urban expansion, this greening from the vegetation recovery will be hidden, without considering the change from a process perspective. Third, the information on when the vegetation change occurs and what status the vegetation (greening or browning) exhibits cannot be extracted by monotonic trends analysis, while the status of vegetation growth in each sub-period, particularly in the recent period that varied in the length of duration, offers insights for managers to evaluate policy effectiveness and create subsequent strategies [9,34]. Consequently, in addition to the monotonic trends in a certain time interval, the urban vegetation dynamics include the information on when the vegetation change occurs, how frequently the change happens, and what the changing process exhibits. An approach needs to be developed to measure the urban vegetation dynamics from the process perspective.

To characterize the urban vegetation dynamics from the process perspective, it needs to accurately detect the timing of change occurrence and track the time-varying trends at the same time, which is not easy because urban landscapes are highly heterogeneous and temporally dynamic, and the change trajectories of land cover in urban are various [21,35,36,37]. With the opening of the Landsat archive and the free use of the Sentinel data, several time-series-based approaches have been developed to help overcome these challenges [26,38]. The approaches using the yearly imagery derived from the compositing technique, such as Trends in Disturbance and Recovery (LandTrendr) [25] and Vegetation Change Tracker (VCT) [39], have abilities to capture the timing information on when the land cover abruptly changes [40,41]. In recent years, the approaches, including the Continuous Change Detection and Classification (CCDC) [42] and COntinuous monitoring of Land Disturbance (COLD) [43], move forward to use all the available data to monitor the inter-, intra- and seasonal change, and capture the timing of the change. Benefit from the Google Earth Engine (GEE), a free cloud-based platform consists of multi-petabyte publicly available geospatial datasets that have been widely used for the geospatial analysis on a variety of high-impact at Earth-scale [44], the issue of frequent cloud coverage that commonly affects the detection accuracy of vegetation dynamic when using single-image approaches can be largely addressed by optimizing the input of image composition, such as seasonal composition strategies, statistical operators, and band composition [45]. Additionally, several applications have also been developed on the GEE for improving change detection, such as changes based on the time-series urban land cover with high-resolution (15 m) using pan-sharpened Landsat imagery [46].

Building upon the vegetation dynamics research, this study aims to further contribute to this field by (1) developing an approach that integrates trends with dynamic analysis to quantify the urban vegetation dynamics from the process perspective using temporally dense Landsat imagery and (2) applying this approach in the megacity of Shenzhen, as an example, to investigate the vegetation dynamics from 2000–2020, including the information on when and how frequent the change occurs, and what the changing process exhibits.

2. Materials and Methods

2.1. Study Area

Shenzhen covers an administrative area of 1997 km², with a population of 13.43 million and GDP of 2692.71 billion RMB in 2019 [47]. Forty years ago, it used to be a village, with only a population of 0.31 million and GDP of 0.2 billion RMB in 1979. Because of the dramatic urban expansion, the size in 2017 was approximately 35 times to that in 1978 and the cover percentage of developed land in 2017 had increased to 48%, which was larger than the together of forest, grass, and cropland, 44% [48]. In order to keep the urban sustainable, the local government released an urban land-use management policy in 2005, among which an ecological protection control baseline is delineated to limit the urban expansion and restore the ecosystem function (Figure 1) [49]. The vegetation coverage of the entire Shenzhen area increased from 48% to 62% in 2000–2018. The forest-cover percentage, in particular, decreased firstly in 2000–2003 and then increased in 2003–2018 [50]. In this study, we defined the areas located outside the ecological protection control baseline as urban built-up areas.

2.2. Aquisition of Time-Series Landsat Imagery

We selected the EVI, a commonly used index and generated from Near-Infrared, Red, and Blue bands of Landsat imagery [51], to quantify the vegetation dynamics. First, the EVI optimized the vegetation signal with improved sensitivity in high biomass regions, and thus can successfully address the saturated signals problem over high biomass conditions that are common for other ratio-based indices [51]. Second, the EVI decouples the canopy background signal and reduces the atmosphere influences by using the atmosphere-sensitive blue band to correct the red band. The calculation of EVI can be found in Equation (1).

$$\mathrm{EVI}=G\times \frac{{\rho }_{NIR}+{\rho }_{red}}{{\rho }_{NIR}+C_1{\rho }_{red}-C_2{\rho }_{blue}+L}$$

(1)

where ${\rho }_{NIR}$, ${\rho }_{red}$ and ${\rho }_{blue}$ are the surface reflectance of the Near-Infrared, Red, and Blue bands which are atmospherically or partially atmospherically corrected. The coefficients C₁ and C₂ are the aerosol resistance term and L is the adjustment for vegetation canopy background. The values of G, C₁, C₂, and L adopted in the EVI are 2.5, 6, 7.5, and 1, respectively [52,53].

To maintain the quality and consistency for the time-series observation, the EVI was calculated based on the Surface Reflectance (SR) product derived from the Landsat Level-1 data because the SR product improves comparisons between multiple images over the same region by accounting for atmospheric effects such as aerosol scattering and thin clouds [43]. We selected the Landsat 5 and 7 Surface Reflectance Tier 1 that atmospherically corrected using the LEDAPS and Landsat 8 Surface Reflectance Tier 1 that atmospherically corrected using the LaSRC from GEE. The cloud cover and cloud shadow pixels were excluded using the Quality Assessment (QA) band that integrated into the Tier 1 collection. Because of the relatively low temporal frequency (16-day revisit capability), the number of available clear observations, referring to pixels that are free of cloud, cloud shadow, and snow, is reduced. Such a data preprocessing scheme significantly reduces the reliably representative annual values generated by the image compositing method, such as annual median EVI and monthly maximum EVI, for every year [54]. As a result, we used all the clear observations in 2000–2020, instead of the best-pixel selection, to construct the temporally dense EVI time-series. In this study, a total of 1702 images acquired in 2000–2020 were selected from Landsat 5, 7, and 8. The spatiotemporal variation revealed by the percentage of the clear observations can be found in Figure A1.

Based on the EVI time-series, the unique phenological characteristic of vegetation, represented as a harmonic function with periodical change, can be observed [42]. When the vegetation change occurs, the EVI time-series will be out of the original shape. If the land cover with vegetation converts to other types, the EVI abruptly decreases and the periodical property disappears on the time-series. If the vegetation suffers from diseases, pests, and typhoons, or the vegetation in a pixel partially changes, the inter-annual change reflected by the variation in EVI trends occurs. If the vegetation species changes or the climate condition impacts the vegetation phenology, the periodical property on the EVI time-series changes.

2.3. Change Detection for Urban Vegetation

To investigate the vegetation dynamics from the process perspective, we first detected the timing and frequency of vegetation change by the approach of Continuous Change Detection and Classification (CCDC) [23,42] based on the EVI time-series constructed from the temporally dense Landsat imagery in 2000–2020. The CCDC is a robust approach for monitoring land cover/land-use change and generating maps for any given time using dense time-series satellite imagery. This algorithm contains three components: seasonality, trend, and break. The seasonality corresponds to the seasonal change, mostly caused by the vegetation phenology that is led by seasonal patterns of environmental factors, such as precipitation and temperature. The trend focuses on gradual change, driven by vegetation growth and degradation, climate change, and pests. The break refers to abrupt change, generally caused by land cover change. The study from Zhu et al. (2014) that first presented this approach by input all the Landsat bands to detect the various types of land cover change. Recently, many studies also demonstrated that the CCDC is an adaptive approach in that the single band, such as the NDVI or EVI, can be used to detect the land cover change for a single type [55]. The model of CCDC is shown in Equation (2):

$$\widehat{\rho }\left(i,x\right)=a_{0,i}+a_{1,i}\cos \left(\frac{2\pi }{T}x\right)+b_{1,i}\sin \left(\frac{2\pi }{T}x\right)+c_{1,i}x$$

(2)

where $\widehat{\rho }\left(i,x\right)$ is the predicted value for the i-th band of Landsat at Julian x, $a_{0,i}$ is the coefficient for overall value for the i-th Landsat band, $a_{1,i}$ and $b_{1,i}$ are the coefficients for the intra-annual change for the ith band, $c_{1,i}$ is the coefficients for inter-annual change for i-th band, and T is the number of days in one year, which equals to 365.

For each pixel, the CCDC estimated the seasonal, gradual, and abrupt changes for vegetation and recorded the positions of the change occurrence on the EVI time-series. Based on these positions, the CCDC identified the year of vegetation change occurrence and counted the number of the position for quantifying the vegetation change frequency. It has been confirmed that the CCDC is a computationally expensive method and needs huge data storage because it requires a high temporal frequency of clear observation and it updates the time-series model once the new observations are new available [42]. However, benefitting from the powerful computational infrastructure from the GEE, these issues can be addressed and the CCDC algorithm, tabbed as ee.Algorithms.TemporalSegmentation.Ccdc under the Earth Engine algorithms documentation, has been officially released (https://developers.google.com/earth-engine/apidocs/ee-algorithms-temporalsegmentation-ccdc?hl=en, accessed on 10 June 2021). Here, we implemented the CCDC algorithm to detect the vegetation dynamics and lay out the frequency and the latest change year.

2.4. Trends Analysis and Changing Processes Characterization

After detecting the year of change for each pixel, the EVI time-series was segmented into several fragments based on the change year, and the trends in each fragment (also called change period) were calculated by the Ordinary Linear Square (OLS) regression model (Equation (3)). As the change frequency of each pixel is different and the year of change occurrence is time-varying, the duration of each fragment varied in length.

$${\mathrm{EVI}}_{\mathrm{t}}=ax+b$$

(3)

where a is slope, b is constant, t represents the observation time, and the EVI_t is the model estimated EVI value for a pixel at time t.

Based on the trends in each fragment, three kinds of trends, nochange, browning, and greening, respectively assigned as 1, 2, and 3, were identified and the processes of vegetation dynamics can be exhibited from the matrix. For example, if there were two changes, three intervals were segmented and the trends in each fragment were filled in the matrix as shown in

Table 1. Matrix for processes of vegetation dynamics based on the trends in each fragment with different duration.

		Status in the Second Fragment	Status in the Third Fragment
			1	2	3
Status in the first fragment	1	1	111	112	113
		2	121	122	123
		3	131	132	133
	2	1	211	212	213
		2	221	222	223
		3	231	232	233
	3	1	311	312	313
		2	321	322	323
		3	331	332	333

. The class of 123 in the matrix means the process of the vegetation dynamics contains three stages, which change from nochange in the first fragment to browning in the second fragment, and then turning green in the recent period. Here, although all of the pixels had the same investigated period, 2000–2020, each pixel had different numbers of fragments and each fragment had a different duration. As a result, the sub-periods in the process of vegetation dynamics were not fixed. The area proportion of each sub-period can be filled in the matrix to evaluate the how the urban vegetation changes and what the recent status presents. The overview of the integrated approach and its application can be found in Figure 2.

3. Results

3.1. Non-Monotonous Trends of Vegetation Greenness in the Urban Area

We summarized the vegetation dynamics along with the urban outward expansion and internal redevelopment based on all of the available EVI values calculated from the time-series Landsat imagery in 2000–2020 (Figure 3). The EVI trajectories showed that the vegetation greenness followed the non-monotonous trends and the processes varied from different areas. Such non-monotonous trends were sometimes misleading from the statistical significance and the positive or negative values.

As shown in Figure 3a, the pixel experienced a changing process of urban expansion, in which the EVI abruptly decreased in 2012 and continuously increased in 2014–2020. Although a browning trend was found on Line AB with a slope of −0.002 (p = 0.156), two hidden greening trends, reflected from a positive slope of 0.013 (p = 0.000) in 2000–2012 and another positive slope of 0.031 (p = 0.000) in 2014–2020, were revealed from the process-based analysis.

The changing process of urban redevelopment is summarized in Figure 3b–d, although the EVI values are from different residential areas. When no land cover change occurs and the vegetation continuously grows, a monotonous trend can be simply evaluated. Figure 3b shows a significant increase of the EVI trend on Line AB with a slope of 0.008 (p = 0.000). When a residential/commercial/industrial area is dismantled, an abrupt change occurs on the EVI time-series. Figure 3c shows an insignificant trend presented by Line AB (0.001, p = 0.203) during 2000–2020, while a significantly decreased trend in 2018–2020 was found on Line CD, with a negative slope of −0.049 (p = 0.000). When the residential/commercial/industrial area is redeveloped, the trend will turn green along with the tree replanting. Figure 3d shows a vegetation recovery process, due to which the EVI experienced fluctuant changes of browning and greening, and finally turned green from 2015 with a positive slope of 0.012 (p = 0.008).

3.2. Spatial Pattern of the Vegetation Dynamics from a Process Perspective

The results of vegetation dynamics from the maps and statistical analysis showed 35.3% of the entire Shenzhen saw at least one-time change and the frequent changes were mostly distributed in the built-up areas, located from the central to the northeast urban areas (Figure 4a,b). Comparing the spatial extents of each level on change frequency, the cover percentages, in order of the change occurrence number, were 22.0%, 9.2%, 3.0%, and 1.1%.

Integrating the trends with dynamics, the statistical results showed 71.2% of all of Shenzhen was greening, among which 49.5% showed consistent trends in greening during the period of 2000–2020 and 21.7% revealed its status of turning green in the recent change period, which began from the latest year of change occurrence detected by the CCDC algorithm to the year of 2020. The results also showed only 13.7% was browning, including 5.2% of persistent browning and 8.5% turning brown. Here, we only display the matrix for the processes of vegetation dynamics with persistence, one-time, and two-time changes caused by land cover, seasonal, and trend changes because their area proportion was over 95%. After summarizing the matrices with one- and two-time changes, three major processes included browning-to-turning green assigned as 23, 123, 223, 233, and 323; no change-to-turning green assigned as 13,113, and 133; and probably turning green assigned as 131 and 231. The areas were 99.0 km², 108.3 km², and 20 km², respectively (Figure 5b,c). It was also found that 106 km² exhibited green loss, including types of 31, 32, 311, 312, 321, 322, 331, and 332.

3.3. Trends in Recent Change Period for Vegetation Greenness

Focusing on the recent change period, the latest shifts in the trends mostly appeared before 2012, with 66.4% of Shenzhen city with change, indicating a large number of long-term changes occurred in 2000–2020 (Figure 6b). Except for the persistent change, the relatively long-term changes were primarily located in the built-up area, while the short-term changes varied greatly from different areas. The spatial pattern of the latest year of change occurrence showed those short-term changes clustered in four regions, including the western coastal area, central but besides the ecological protection control baseline, the northeast build-up area, and the zone around the northern reservoirs (Figure 6a).

According to the duration of the latest period that starts from the latest year of change occurrence to the year 2020 and the corresponding trends in the latest period, the long-term changes mostly had positive values of EVI slopes and the area proportion of greening in Shenzhen decreased from 2.3% in 2001–2020 to 0.2% in 2019–2020, while the short-term changes showed more insignificant EVI trends with an increased percentage of 0.4% in 2012–2020 to 1.8% in 2019–2020 (Figure 6b). The statistical results of time-varying trends from the boxplot showed those long-term observations had slower change rates. Interestingly, the data distribution of the slopes with positive and negative values were similar, although the amount of them was quite different (Figure 6e,f).

4. Discussion

4.1. The Application of the Temporally Dense Imagery for the Process Identification in Urban Vegetation Greenness

Relatively few studies investigated the urban vegetation dynamics from the process perspective. Our results showed the non-monotonous trends of the vegetation dynamics in urban areas were more complicated because the frequent changes generated several turning points in the EVI time-series, such that the greenness experienced multiple processes (Figure 3 and Figure 4). These multiple processes skewed the trends derived from the linear regression model, misleading the identification of greening and browning via statistical significance and the slope values of trends (Figure 3 and Figure 4). For example, the EVI slope in Figure 3a calculated by the linear regression model exhibited an insignificant decreased trend, with a slope of −0.002. While the process-based analysis showed the vegetation dynamics experienced three stages, including greening in the first stage, insignificant change in the second stage, and greening again in the latest stage, with the slope values of 0.013 (p = 0.000), −0.010 (p = 0.571), and 0.031 (p = 0.000) (Figure 3a). Such effects led by the highly dynamic landscapes, however, can only be improved by the process-based analysis using the temporally dense imagery because the time-series approaches, such as CCDC in this study, segment the EVI time-series into fragments and then evaluate the trends in each sub-period from the piecewise perspective [9,23].

The trends analysis using the temporally dense Landsat data allows more subtle understandings of the urban vegetation change, compared to the composition approach applied to Landsat data that follows a best-pixel selection [56]. As the revisit period of Landsat 5, 7, and 8 was 16 days, the strategy of best-pixel selection, such as annual median EVI or monthly median EVI, will reduce the reliably representative annual value for every year by the Landsat data after removing the cloud cover and cloud shadow pixels. The advantage of the temporally dense data accelerates the detection of structural changes, such as vegetation phenology [57,58] and lake/reservoir shoreline change [59]. The outliers resulting from the simple criteria such as the maximum NDVI or mean NDVI in the growing season based on the time-series can be largely eliminated, which is especially applicable to the regions where the data is frequently absent due to the cloudy and rainy conditions. Additionally, except for the abrupt changes caused by conversions among different land cover types, the inter- and intra-annual changes reflected from the time-series spectral characteristics can only be detected based on the temporally dense imagery [34].

4.2. Process and Recent Trends for the Urban Vegetation Greenness

The information provided from the process analysis, including when the change occurs, how frequently the change happens, and what the changing process experiences, help to better evaluate the vegetation quality and ecosystem function. Compared to most concerns of the vegetation dynamics, including the trends before and after disturbances in forest and grassland, urban ecologists prefer the multiple changes, such as repeated conversions from barren land to grassland in the urban redevelopment process, because the frequent change potentially influences the urban ecosystem properties (e.g. biodiversity and provision of ecosystem services) [26,60,61,62]. The process matrix showed 71.2% (versus 66.7% evaluated from the linear regression model, Figure 3) of Shenzhen was greening, including not only 49.5% of consistent green but also 21.7% of turning green in the latest period. This suggests that the temporal decomposition reveals more information which is frequently ignored and hidden in the total trends quantified by the linear regression model (Figure 5 and Figure 6d). The finding also showed the processes varied greatly even for the greening was found in the latest period (Figure 5). For example, 99 km² had experienced the process of browning first and then turning green after a period of insignificant change, while 108 km² had experienced the process of insignificant change first and then turning green. These variations display different stages of urban development, which are beneficial to infer the possible drivers.

The temporally variant trends in the latest change period that starts from the latest change year detected by the CCDC algorithm and ends in the year 2020 exhibit the recent status of vegetation growth and provide opportunities to accurately evaluate the provision of ecosystem services, although the durations varied from different areas (Figure 6). For example, the pixels in the northeast built-up area showed greening in the recent period, suggesting the dwelling environment was improved, although those pixels used to be forest and the urban expansion abruptly decreased the EVI [49]. Such time-varying trends are especially appropriate for the comparison of the ecosystem service provision between the old and new urban areas because the effects of urban outward expansion can be exclusive and can be used to better project the vegetation change.

More importantly, the vegetation dynamics investigated from the process perspective will help government agencies better assess the effectiveness of urban management policies and then implement subsequent planning strategies. The first outcome, change frequency, from the dynamic analysis can be used to evaluate the ecosystem stability. The spatial pattern of change frequency showed the vegetation growth in the eastern area was more stable than that in the western area, although both of them were located inside the ecological protection baseline. When combining the changing process, timing information, and the trends in the recent change period, a large number of dynamics were mostly caused by vegetation recovery after the year 2005 but before 2012, suggesting the previous protection policy, delineating the baseline in 2005, was effective [49]. This forest recovery led to fast greening in a short time interval (Figure 6); the ecosystem continues to function, however, relies on the stability that copes with the obstacles arising from nature and changing environment [63,64]. As a result, how to improve the ecosystem stability, such as increasing the vegetation diversity and making appropriate management schemes, becomes a further consideration for urban planning to achieve the goal of sustainability.

4.3. Limitation and Future Work

In this study, we implemented the CCDC algorithm to detect the vegetation dynamics only based on the EVI time-series constructed from the Landsat SR imagery collection. However, previous studies showed the change detection accuracy will be improved when all of the spectral bands are used as input [42]. Comparing to the Landsat image, the Sentinel-2 data had finer spatial resolution (10 m and 20 m) and temporal resolution (5-day). Additionally, the wavelengths of Sentinel data are more sensitive to the chlorophyll content and phenological states [45,65,66]. It will be more promising when using the Sentinel data to quantify the urban vegetation dynamic. Our results from the process-based analysis showed the Shenzhen city is greening but varies greatly in different areas (Figure 6c,d). Several studies suggested the local climate and the economic input, including irrigation and fertilization, are related to urban greening [31,67], while the economic levels are varied from different regions, which may impact the vegetation greenness. Therefore, more studies need to be performed across many other cities to test and explore the drivers behind urban greenness.

5. Conclusions

This study presented a method that integrates the trends with dynamic analysis and then applied it to investigate the vegetation dynamic from a process perspective in the megacity of Shenzhen by using the temporally dense Landsat imagery. Compared to monotonic trend analysis, this integrated approach has the ability to provide more detailed dynamic information on when and how frequently vegetation change occurs, and what the changing process is. Using the integrated approach, we found that Shenzhen was turning green in 2000–2020, even though a large-scale urban expansion occurred during this period. Nearly half of the city (49.5%) showed consistent trends in greening, most of which were located inside the ecological protection baseline. We also found that 35.3% of Shenzhen experienced at least one abrupt change in vegetation greenness, among which 61.5% of these lands showed trends in turning green in the recent change period. Such information on vegetation greenness, frequently hidden by the trends analysis alone, can be revealed by our approach that integrates the trends with a dynamic process and provides insights for urban vegetation management.