Correlations between Urbanization and Vegetation Degradation across the World’s Metropolises Using DMSP/OLS Nighttime Light Data

Abstract

Changes in biodiversity owing to vegetation degradation resulting from widespread urbanization demands serious attention. However, the connection between vegetation degradation and urbanization appears to be complex and nonlinear, and deserves a series of long-term observations. On the basis of the Normalized Difference Vegetation Index (NDVI) and the image’s digital number (DN) in nighttime stable light data (NTL), we delineated the spatiotemporal relations between urbanization and vegetation degradation of different metropolises by using a simplified NTL calibration method and Theil-Sen regression. The results showed clear and noticeable spatiotemporal differences. On spatial relations, rapidly urbanized cities were found to have a high probability of vegetation degradation, but in reality, not all of them experience sharp vegetation degradation. On temporal characteristics, the degradation degree was found to vary during different periods, which may depend on different stages of urbanization and climate history. These results verify that under the scenario of a vegetation restoration effort combined with increasing demand for a high-quality urban environment, the urbanization process will not necessarily result in vegetation degradation on a large scale. The positive effects of urban vegetation restoration should be emphasized since there has been an increase in demand for improved urban environmental quality. However, slight vegetation degradation is still observed when NDVI in an urbanized area is compared with NDVI in the outside buffer. It is worthwhile to pay attention to landscape sustainability and reduce the negative urbanization effects by urban landscape planning.

1. Introduction

In the past 50 years, while cities have experienced population growth and land consumption, the pressures on vital ecosystem functions have escalated rapidly [1,2]. Urbanization, which means a population shift from rural to urban areas, is often accompanied by urban expansion and land use change. To interpret the land consumption pressures, the impact of urbanization has been documented in the growing literature on the urban–rural gradient, which shows consistent changes in species richness and species composition [3,4,5]. With worldwide land cover change, biogeochemical cycles, hydrologic systems, and climate and biodiversity change driven by urbanization, increasing numbers of ecologists have accepted that urban areas are hot spots that drive environmental change on multiple scales [6]. These hot spots are estimated not only to be current threats for ecosystems but also will probably last for a long period in developing nations. The global urban population will reach 5 billion by 2030 [7] and will increase by 2.7 billion, nearly doubling today’s urban population of 3.4 billion by 2050 [8], which indicates that the dramatic urbanization phenomenon will continue. Therefore, the threat of changes in biodiversity with an increase in global urbanization is a concern that needs to be brought to the foreground [9].

The total vegetation in urban and suburban areas is an important indicator of urbanization pressure on biodiversity, because plants can be lost during either the initial habitat transformation or the landscape fragmentation processes [5]. Vegetation degradation, which means a reduction in the available biomass, often represents as a decline in the vegetative ground cover. Vegetation degradation triggered by a change of land cover to impervious surfaces in urban areas may result in eco-environmental threats with net primary production reduction and surface temperature variation [10,11]. However, some opposite perspectives have recently put forward a different relation, which suggest some urbanization factors can enhance urban vegetation activity [12]. The plant growth in urban areas might be promoted by warmer temperature and greater tropospheric CO₂ concentration [13,14]. Because the value of urban ecosystem services has been deeply rooted in landscape management [15], green infrastructure investments are supported and some negative urbanization effects are mitigated by landscape planning [16,17]. On balance, the interconnection between vegetation degradation and urbanization seems to be a complex and nonlinear system [12].

Accordingly, a hypothesis can be proposed that the spatiotemporal relations between urbanization and vegetation degradation are diversified and may relate to the stage of urbanization and geographical location. However, a systematic evaluation of the spatiotemporal trends of vegetation activities across multiple cities over large areas is still lacking [12]. The Normalized Difference Vegetation Index (NDVI) has been adopted as an effective index for describing the dynamics of urban vegetation [18,19]. As a remote sensing index, NDVI can indirectly estimate gross and net primary productivity, biomass, and green leaf area in a variety of grassland and forest ecosystems [20,21,22,23,24]. If a regular distribution of NDVI variation rate on the urban–rural gradient exists, the landscape is likely to experience a negative effect related to urbanization on vegetation. Therefore, in this article, the NDVI variation rate is adopted as a vegetation degradation indicator that can be extracted across multiple cities over large areas.

The Defense Meteorological Satellite Program (DMSP)/Operational Line-Scan System (OLS) nighttime stable light data (NTL) s been demonstrated to be an effective data source for mapping urban expansion and urbanization dynamics at a large spatial scale [25,26,27,28,29,30]. NTL not only can be applied to urban area extraction, but also widely reflect population density, economic activity, and energy use in an urban area [31,32,33,34,35]. Therefore, the image’s digital number (DN) may take on a growing trend in rapidly urbanized areas of developing countries. However, the DN may seem relatively stable in the core urban areas of developed countries, because of saturation at the highest value. To find the differentiation of nighttime light change among various cities, urbanization trends derived from DN were calculated and compared with the NDVI variation rates in the same geographical extent.

In light of the hypothesis that the spatiotemporal relations between urbanization and vegetation degradation should show differences, the study presented here contains four major parts. (a) To develop a simplified calibration method on NTL series to extract the urbanization area of the world’s metropolises; (b) To detect the variation trends of vegetation and urbanization by Theil-Sen regression [36,37]; (c) To classify the metropolises and differentiate between different relations in the variation trends; (d) To discuss the impact of human factors on vegetation changes in urbanization.

2. Data

2.1. Data Sources

A long-term NDVI series was obtained from the Vegetation Index and Phenology (VIP) Research Lab at the University of Arizona, with a spatial resolution of 0.05° (approximately 5.6 km). For this data set, multisensor data sets include the Advanced Very High Resolution Radiometer (AVHRR), the Moderate Resolution Imaging Spectroradiometer (MODIS), and the Visible Infrared Imager Radiometer Suite (VIIRS), and daily surface reflectance data were fused to build a homogeneous vegetation cluster model [38]. Then, monthly MODIS NDVI composed by MOD09GA as the MODND1M data set was downloaded from the Geospatial Data Cloud, Chinese Academy of Sciences, as a cross-reference to VIP NDVI. The maximum value composite (MVC) method was then adopted to derive the annual maximum NDVI in each year [39,40]. Version 4 of the NTL data was obtained from National Geophysical Data Center at the National Oceanic and Atmospheric Administration , with the background noise in the composite defined as zero and the final data values ranging from 1 to 63 [41]. The NTL data set included data acquired by six different DMSP satellites: F10, F12, F14, F15, F16, and F18 with 30-s grids, spanning −180° to 180° longitude and −65° to 75° latitude [42]. Furthermore, the global land cover map, GlobCover 2009, was downloaded from European Space Agency (ESA) as a cross reference to NTL data. The data specifications are listed in Table 1.

Table 1. Data specifications of the four data sets adopted in this study.

Products and Sensors	Time Period	Resolution	Data Source	Processing
DMSP/OLS NTL	1992 to 2012	Yearly/30-s grids	National Geophysical Data Center at the National Oceanic and Atmospheric Administration	Averaging the pixel values of each city’s urbanization area to derive the annual DN
VIP NDVI	June 1981 to December 2010	Monthly/5.6 km	Vegetation Index and Phenology Research Lab at the University of Arizona	MVC method to derive the annual maximum NDVI
MODIS MODND1M NDVI	February 2000 to December 2010	Monthly/1 km	Geospatial Data Cloud, Chinese Academy of Sciences	MVC method to derive the annual maximum NDVI
GlobCover 2009	2009	Single year/300 m	European Space Agency (ESA)	Visually compare with NTL

Table 1. Data specifications of the four data sets adopted in this study.

Products and Sensors	Time Period	Resolution	Data Source	Processing
DMSP/OLS NTL	1992 to 2012	Yearly/30-s grids	National Geophysical Data Center at the National Oceanic and Atmospheric Administration	Averaging the pixel values of each city’s urbanization area to derive the annual DN
VIP NDVI	June 1981 to December 2010	Monthly/5.6 km	Vegetation Index and Phenology Research Lab at the University of Arizona	MVC method to derive the annual maximum NDVI
MODIS MODND1M NDVI	February 2000 to December 2010	Monthly/1 km	Geospatial Data Cloud, Chinese Academy of Sciences	MVC method to derive the annual maximum NDVI
GlobCover 2009	2009	Single year/300 m	European Space Agency (ESA)	Visually compare with NTL

2.2. Study Areas

We selected 50 large metropolises with wide areas to represent the most clearly urbanized places around the world. To diversify the distribution, we dispersed the metropolises across the continents (Table 2). The city selection was based on a comprehensive qualitative recognition for population, economy, climate, country, and distance from each other. If a number of big cities were in accordance with the climate condition or economic level, only the best known city was selected. For example, Mumbai, Dhaka, Karachi, and New Delhi all are similar metropolitan areas, but New Delhi is best known in South Asia because it is the capital of India. Therefore, representativeness was an important factor in our selection. Finally, no more than three metropolises were chosen from the same country, except for the United States, which includes 14 metropolitan areas with large tracts of built-up land and high nighttime light dispersed around the locale.

Table 2. The distribution of the selected metropolises.

Continent	Metropolis
Asia	Beijing, Shanghai, Guangzhou, Taipei, Singapore, Bangkok, Dubai, New Delhi, Tehran, Tokyo, Kyoto, Seoul
Europe	London, Liverpool, Berlin, Athens, Lisbon, Madrid, Barcelona, Rome, Milan, Paris, Brussels, Stockholm, Moscow, Saint Petersburg, Istanbul
North America	Mexico City, New York, Miami, Houston, Dallas, Phoenix, Atlanta, Los Angeles, St. Louis, Washington D.C., Cleveland, Detroit, Boston, Chicago, Minneapolis, Toronto, Montreal
South America	Buenos Aires, Sao Paulo, Rio de Janeiro
Africa	Johannesburg, Cairo
Oceania	Melbourne

Table 2. The distribution of the selected metropolises.

Continent	Metropolis
Asia	Beijing, Shanghai, Guangzhou, Taipei, Singapore, Bangkok, Dubai, New Delhi, Tehran, Tokyo, Kyoto, Seoul
Europe	London, Liverpool, Berlin, Athens, Lisbon, Madrid, Barcelona, Rome, Milan, Paris, Brussels, Stockholm, Moscow, Saint Petersburg, Istanbul
North America	Mexico City, New York, Miami, Houston, Dallas, Phoenix, Atlanta, Los Angeles, St. Louis, Washington D.C., Cleveland, Detroit, Boston, Chicago, Minneapolis, Toronto, Montreal
South America	Buenos Aires, Sao Paulo, Rio de Janeiro
Africa	Johannesburg, Cairo
Oceania	Melbourne

The spatial extent of urbanization is usually dependent on data source and research goals. Two items should be noted that DMSP/OLS cannot precisely detect impervious surfaces in urban environments, and the urbanization processes can occur outside the built-up land as well. In this article, for the sake of vegetation restoration analysis, the boundary should be a little broader than the size of the impervious surface. We set the DN in the F182012 descend from 63 at 1 interval and compared the area with land cover images from GlobCover 2009 or Google Earth. The stepwise comparison showed that when the DN descended to 50, the area contained all the continuous impervious surface area in the 50 metropolises. So in this article, DN > 50 in F182012 was set as the urbanization boundary. It was a little larger than the impervious surface area, with relatively high population/economic density. By this extraction, some metropolises comprise more than one large city and virtually perform as holistic urban clusters. For example, New Delhi, Shanghai, Guangzhou, and Taipei, respectively, represent New Delhi and Delhi, Shanghai and Suzhou, Guangzhou and Shenzhen, and Taipei and Xinzhu. That is to say, we only take the name of the largest city to name the metropolitan areas for easy recognition in this article (Figure 1). Then the annual DN series of the 50 metropolises were built by averaging the pixel values in each city’s urbanization area where DN > 50 in F182012.

Figure 1. The boundary of three metropolises: (a) DMPS/OLS NTL in 2012; (b) GlobCover 2009 land cover map. The legends of the land cover types were merged for the display; (c) Image from Google Earth in 2014.

Figure 1. The boundary of three metropolises: (a) DMPS/OLS NTL in 2012; (b) GlobCover 2009 land cover map. The legends of the land cover types were merged for the display; (c) Image from Google Earth in 2014.

3. Methods

NTL data cannot be used directly to extract the dynamics of urbanization, because the DN value of NTL data is the average value from archived data without any on-board calibration [30,43]. To confirm the continuity and comparability of NTL data, former research was consulted, and the intercalibration method was improved to correct the data systematically for urbanization detection [30,42,43,44]. In this article, a simplified intercalibration method was used before the DN yearly series was generated. In the trend judgment part, the former research indicated that the unbiased predictor of the Theil-Sen approach was suggested as a potential replacement for ordinary least squares for linear regression in remote sensing applications [45]. Therefore, the Theil-Sen median trend is used in this article to quantify the variation of NDVI and DN. Finally, in the overlay analysis part, after the trend maps are overlaid and the different urbanization types are classified, we can determine whether the urbanization process definitely results in vegetation degradation. This procedure can be simplified as in Figure 2.

Figure 2. The flowchart for this research.

Figure 2. The flowchart for this research.

3.1. Calibration and Composition of NTL

DN values from the same satellite for different years or from different satellites for the same year often show discrepancies or abnormal fluctuations, which are partly attributable to unstable spill light [30]. As Figure 3 shows, the post-urbanization stage for New York and swift urbanization stage for Beijing are evident. However, in Melbourne and Rio de Janeiro, the trend is in fluctuation and is confused because increasing and decreasing trends may appear during the same period from different satellites. To make the result credible, intercalibration is needed to detect the brightness changes across the time series [46]. Faced with this obstacle, Elvidge et al. (2009, 2014) developed a second-order regression model, with F121999 used as the reference composite and with Sicily chosen as the reference region [43,46]. The form of the calculation is Y = C₀ + C₁X + C₂X², and calculated values that run beyond 63 are truncated at 63. Although some issues are debatable, it is generally reasonable to be focused on the whole statistical value rather than particular pixels [33].

After comparing the 50 metropolitan areas, we found the DN in New York most stable. Because the spatial extent of urbanization has been extracted by F182012, and most of the pixel values in this extent are higher than 30 in the former years, F182012 is used as the reference composite, and DN > 30 in New York is chosen as the reference region. Also, because we only deal with data ranging from 30 to 63 instead of 0–63, a second-order regression model may be insufficient. The inflection point in the quadratic curve may not appear in this half of the data range. In fact, we find no significant difference between the R² from the second-order regression model and the first-order regression model. Thus we simplified the calculation as Y = aX + b, where X means the DN in the NTL data set and Y means the DN after intercalibration. The intercalibration applied on the basis of the offsets and coefficients is listed in Table 3.

Figure 3. Examples of DN values in the urbanization area from different satellites before calibration during 1992–2012.

Figure 3. Examples of DN values in the urbanization area from different satellites before calibration during 1992–2012.

Table 3. Coefficients of the simplified linear regression models for NTL from 1992 to 2012.

Satellite	Year	a	b	R2	Satellite	Year	a	b	R2
F10	1992	0.50	26.61	0.74	F15	2000	0.57	23.28	0.83
	1993	0.48	28.80	0.75		2001	0.64	19.89	0.84
	1994	0.54	23.70	0.72		2002	0.70	15.93	0.87
F12	1994	0.52	26.46	0.82		2003	0.52	27.78	0.88
	1995	0.60	21.64	0.79		2004	0.50	29.38	0.89
	1996	0.63	19.66	0.81		2005	0.51	29.03	0.88
	1997	0.52	26.30	0.77		2006	0.46	32.07	0.92
	1998	0.53	26.38	0.79		2007	0.44	33.15	0.88
	1999	0.73	13.71	0.82	F16	2004	0.51	27.41	0.82
F14	1997	0.44	32.16	0.84		2005	0.50	29.57	0.90
	1998	0.47	31.24	0.86		2006	0.53	27.54	0.89
	1999	0.49	29.30	0.87		2007	0.71	16.16	0.89
	2000	0.84	7.65	0.88		2008	0.55	25.79	0.89
	2001	0.62	20.32	0.83		2009	0.65	18.75	0.84
	2002	0.51	27.24	0.79	F18	2010	1.21	−15.91	0.91
	2003	0.62	21.21	0.87		2011	0.56	23.71	0.76

Table 3. Coefficients of the simplified linear regression models for NTL from 1992 to 2012.

Satellite	Year	a	b	R2	Satellite	Year	a	b	R2
F10	1992	0.50	26.61	0.74	F15	2000	0.57	23.28	0.83
	1993	0.48	28.80	0.75		2001	0.64	19.89	0.84
	1994	0.54	23.70	0.72		2002	0.70	15.93	0.87
F12	1994	0.52	26.46	0.82		2003	0.52	27.78	0.88
	1995	0.60	21.64	0.79		2004	0.50	29.38	0.89
	1996	0.63	19.66	0.81		2005	0.51	29.03	0.88
	1997	0.52	26.30	0.77		2006	0.46	32.07	0.92
	1998	0.53	26.38	0.79		2007	0.44	33.15	0.88
	1999	0.73	13.71	0.82	F16	2004	0.51	27.41	0.82
F14	1997	0.44	32.16	0.84		2005	0.50	29.57	0.90
	1998	0.47	31.24	0.86		2006	0.53	27.54	0.89
	1999	0.49	29.30	0.87		2007	0.71	16.16	0.89
	2000	0.84	7.65	0.88		2008	0.55	25.79	0.89
	2001	0.62	20.32	0.83		2009	0.65	18.75	0.84
	2002	0.51	27.24	0.79	F18	2010	1.21	−15.91	0.91
	2003	0.62	21.21	0.87		2011	0.56	23.71	0.76

After the intercalibration, intra-annual composition and inter-annual series corrections are common methods used to remove any unstable or inconsistently lit pixels [30]. Because these pixels may disturb the stability of DN series, we treat them as noisy pixels to be eliminated. We define that the conditional statement DN_n-1 − DN_n ≥ 3 and DN_n+1−DN_n ≥ 3 (n ∈ [1993, 2011]) is probably the noise. There is a low probability that the DN from last year and next year are close, and the intervals between the current year and these before and after years change sharply. This finding is based on the concept of urbanization as a successive process, and the present urbanization stage is related to the past and future. The moving average method is adopted to smooth the series, with three years in succession adopted by an inverted sequence that started from 2011. The noise pixel would be replaced by the mean value of the other DN on the same position in the before and after years. The not-noise pixel would get the mean value of all DN on the same position in these three years. Finally, some examples of the smooth result are shown in Figure 4. By calibration, the tendency may seem more reasonable than the initial NTL data.

Figure 4. Examples of DN values of the urbanization area from different satellites after calibration during 1992–2012. F10–F18 express the original DN value, and the “average” line was the final value after calibration.

Figure 4. Examples of DN values of the urbanization area from different satellites after calibration during 1992–2012. F10–F18 express the original DN value, and the “average” line was the final value after calibration.

3.2. Variation Trend Judgement

We calculated trends on the basis of annual DN and annual maximum NDVI by Theil-Sen slope and used the Mann-Kendall method to test the significance. The Theil-Sen slope estimator is the median of the slopes calculated between observation values at all pairwise time steps [47]. It is calculated between observations X_j and X_i at pairwise time steps t_j and t_i [48]:

$$slope=Median\left(\frac{X_j-X_i}{t_j-t_i}\right)$$

(1)

where for NDVI variation, slope > 0 means restoration and slope < 0 means degradation.

Similar to the Theil-Sen procedure, the Mann-Kendall test examines the slopes between all pairwise combinations of samples, where each data point is treated as the reference for the data points in successive time periods [47]. The indicator of Kendall’s S is defined as [49]

$$S={\displaystyle \sum }_{i=1}^{n-1}{\displaystyle \sum }_{j=i+1}^{n}sign\left(x_i-x_j\right)$$

(2)

$$sign\left(x_i-x_j\right)=\{\begin{array}{ccc}1& if& x_i-x_j<0\\ 0& if& x_i-x_j=0\\ -1& if& x_i-x_j>0\end{array}$$

(3)

where n is the length of the time series, and x_i and x_j are observations at time i and j, respectively. When there is independent and identical distribution between data values, then the variance is given by [50,51]:

$$\text{Var}\left(S\right)=\frac{n\left(n-1\right)\left(2n+5\right)}{18}={\text{σ}}^2$$

(4)

where σ is the standard deviation. Then the equation for Mann-Kendall significance (Z) is computed by

$$Z=\{\begin{array}{ccc}\frac{S-1}{\sqrt{Var\left(S\right)}}& for& S>0\\ 0& for& S=0\\ \frac{S+1}{\sqrt{Var\left(S\right)}}& for& S<0\end{array}$$

(5)

where $\left|Z\right|\ge 1.96$ (equivalent to p ≤ 0.05) is judged as significant.

4. Results

4.1. Vegetation Variation Trend

Evaluating global long-term vegetation trends on the basis of NDVI is not a new topic in Earth observation studies [52,53]. However, when the trend is calculated in different periods, the results may vary significantly. The vegetation degradation and restoration trends during five periods are calculated by the Theil-Sen slope estimator and the Mann-Kendall test (Figure 5). Although the Earth has experienced wide vegetation variation during 1981 to 2010, non-significant variations appeared to be the main process in the later stage from 1996 to 2010. The significant degradation is remarkable in Outback Australia and Southwest Asia. However, because desert is the main cover type in these regions, the absolute quantity of greenness does not change a great deal. No strong evidence indicates that the degradation in these regions have direct relations to urbanization. The east of America, west of Europe, and east of China, where high-density areas consisting of cities are distributed among large volumes of cultivated land, experienced vegetation restoration during 1981–1995 (Figure 5a,d). Many regions have experienced a vegetation restoration process despite a large number of highly urbanized cities in these regions. However, it is unclear whether this restoration has some relation to urbanization.

Figure 5. Global land vegetation degradation and restoration trends indicated by NDVI variation: for (a) 1981–2010; (b) 1986–2010; (c) 1991–2010; (d) 1996–2010; and (e) 2001–2010.

Figure 5. Global land vegetation degradation and restoration trends indicated by NDVI variation: for (a) 1981–2010; (b) 1986–2010; (c) 1991–2010; (d) 1996–2010; and (e) 2001–2010.

4.2. Light Variation Trend

One the basis of the annual DN series from 1992 to 2012, the urbanization trend of every metropolis can be calculated by the Theil-Sen slope (Table 4). To delineate the urbanization stages, we divided them into five classifications by urbanization speed: Rapid, Relatively Fast, Moderate Speed, Relatively Slow, and Sluggish. After ranking the slope, we established that each class has 10 cities in a high-low slope sequence. Most of the metropolitan areas in the Rapid classification are in Asia, whereas a large number of North America and Western Europe metropolises have Relatively Slow or Sluggish classifications. As Figure 6 shows, the urbanization recognized by remote sensing clearly varies. In the post-urbanization stage of London and Boston, the urban sprawl process is weak, which indicates that their economic growth does not rely on urban expansion. However, the swift urbanization stage in East Asia indicates that the economic growth in these cities is accompanied by a substantial increase of impervious surfaces in suburban areas. As the land cover type changes severely, vegetation degradation may be difficult to avoid in this urbanization stage.

Figure 6. Urban expansion in several samples from 1992 to 2012. The definition of urbanization threshold is DN ≥ 50 after calibration, and the extent in the layout is only for the core area of the metropolis or the urban clusters. Guangzhou signifies the Pearl River Delta urban clusters.

Figure 6. Urban expansion in several samples from 1992 to 2012. The definition of urbanization threshold is DN ≥ 50 after calibration, and the extent in the layout is only for the core area of the metropolis or the urban clusters. Guangzhou signifies the Pearl River Delta urban clusters.

Table 4. Classification of the urbanization trend of metropolises during 1992–2012. Each class has 10 cities in a high-low slope sequence.

Classification	Metropolis or the Urban Clusters	Range of Slope
Rapid (R)	Shanghai, Dubai, Beijing, Cairo, Guangzhou, New Delhi, Bangkok, Tehran, Phoenix, Istanbul	0.76–1.58 **
Relatively Fast (RF)	Lisbon, Singapore, Seoul, Moscow, Taipei, Saint Petersburg, Madrid, Atlanta, Dallas, Sao Paulo	0.51–0.74 **
Moderate Speed (MS)	Athens, Milan, Johannesburg, Rome, Houston, Melbourne, Mexico City, Buenos Aires, Barcelona, Rio de Janeiro	0.25–0.51 **
Relatively Slow (RS)	Berlin, Los Angeles, Washington D.C., St. Louis, Chicago, Minneapolis, Kyoto, Miami, Tokyo, Paris	0.13–0.21 **
Sluggish (S)	Toronto, Detroit, Cleveland	0.08–0.13 **
	Brussels, Liverpool, London, Montreal, Boston, Stockholm	0–0.08

** Significance at the p < 0.05 level. The slope was calculated by the Theil-Sen method, and in the Mann-Kendall method, 1.96 was set as the threshold of significance.

Table 4. Classification of the urbanization trend of metropolises during 1992–2012. Each class has 10 cities in a high-low slope sequence.

Classification	Metropolis or the Urban Clusters	Range of Slope
Rapid (R)	Shanghai, Dubai, Beijing, Cairo, Guangzhou, New Delhi, Bangkok, Tehran, Phoenix, Istanbul	0.76–1.58 **
Relatively Fast (RF)	Lisbon, Singapore, Seoul, Moscow, Taipei, Saint Petersburg, Madrid, Atlanta, Dallas, Sao Paulo	0.51–0.74 **
Moderate Speed (MS)	Athens, Milan, Johannesburg, Rome, Houston, Melbourne, Mexico City, Buenos Aires, Barcelona, Rio de Janeiro	0.25–0.51 **
Relatively Slow (RS)	Berlin, Los Angeles, Washington D.C., St. Louis, Chicago, Minneapolis, Kyoto, Miami, Tokyo, Paris	0.13–0.21 **
Sluggish (S)	Toronto, Detroit, Cleveland	0.08–0.13 **
	Brussels, Liverpool, London, Montreal, Boston, Stockholm	0–0.08

** Significance at the p < 0.05 level. The slope was calculated by the Theil-Sen method, and in the Mann-Kendall method, 1.96 was set as the threshold of significance.

4.3. Correlation between Light and Vegetation

The mean value of NDVI for all the pixels in each urbanization area was calculated during 1992–2010, and the relations between mean NDVI and mean DN were diverse. Four samples in Northeast Asia show that when light increases, the variation trend of mean NDVI is uncertain (Figure 7). For these four cities, Tokyo’s classification of urbanization trend is Relatively Slow; Seoul’s classification is Relatively Fast, with the increase mainly during 1992–2002; and the classification of Beijing and Shanghai, which experienced high-speed urban growth, is Rapid. As the scatterplot shows, a positive correlation exists in Tokyo, but a negative correlation exists in Shanghai, whereas in Seoul and Beijing, the two indicators are nearly uncorrelated. The correlations between mean NDVI and mean DN among the world’s metropolises in each year are drawn in Figure 8. In the early years, because many cities did not have such high DN, high positive correlations appear, which mean the higher the light, the more the vegetation biomass. This finding maybe attributed to the location of high DN cities in warm regions with the climate benefits vegetation growth, whereas more evidence is needed to support this conjecture. Nevertheless, the positive correlation declined and negative correlation appeared in later periods, when cities have relatively higher DN value. This phenomenon reflects that the relation between urbanization and vegetation may vary not only in different regions but also in different periods.

Figure 7. The relation between mean NDVI and mean light in the urbanization areas during 1992–2010 in some Northeast Asia cities. Shanghai signifies the Yangtze River Delta urban clusters.

Figure 7. The relation between mean NDVI and mean light in the urbanization areas during 1992–2010 in some Northeast Asia cities. Shanghai signifies the Yangtze River Delta urban clusters.

Figure 8. The correlation between mean NDVI and mean DN in the urbanization areas during 1992–2010 among all the 50 metropolises.

Figure 8. The correlation between mean NDVI and mean DN in the urbanization areas during 1992–2010 among all the 50 metropolises.

5. Discussion

5.1. Delineating the Correlations by Statistical Methods

To further detect the temporal correlation between vegetation degradation and urbanization, three periods are separated: 1981–2010, 1991–2010, and 2001–2010. When the percent of vegetation degradation pixels is calculated in the period 1991–2010 for the defined urbanization area, the degradation percent will drop when the urbanization area is enlarged, which reveals that some large cities may experience vegetation restoration rather than degradation. However, the correlation is relatively weak (Figure 9a). Furthermore, the correlation between the percent of vegetation degradation pixels and the slope of light variation trend is also weak, which means a higher urbanization speed possibly causes more vegetation degradation on a macroscopic scale (Figure 9b). This finding verifies that the urbanization process is not necessarily resulting in vegetation degradation to a great extent. The climate history and urbanization stage may jointly contribute to the uncertain correlation.

Figure 9. The relations among urbanization area, percent of degradation pixel, and slope of light in all 50 metropolises.

Figure 9. The relations among urbanization area, percent of degradation pixel, and slope of light in all 50 metropolises.

To detect the background vegetation variation trend, a 100-km buffer was created from the boundary of every metropolis. Then, the percent of vegetation degradation/restoration between the urbanization area and the outside buffer was contrasted, and the top 10 vegetation variation percentages during different periods were listed (Table 5 and Table 6). During the period 1981–2010 (Table 5), some metropolitan areas underwent obvious vegetation degradation, but others experienced a vegetation restoration process. The degradation may happen in the areas of rapid or relatively fast urbanization, whereas restoration likely happens in areas of moderate or slow urbanization. However, some exceptions exist. For example, Melbourne and Kyoto underwent obvious vegetation degradation, whereas Moscow and Atlanta experienced vegetation restoration processes. These instances maybe attributed more to global climate change than to landscape transformation [54,55]. In addition, during the period 1981–2010, the degradation percentages in outside buffers are much less than the percentage in the urbanization areas, whereas in most of the high restoration cities, the restorations in the urbanization areas are slower than those in outside buffers. This phenomenon may be explained in two ways. On the one hand, the characteristics of vegetation variation may vary significantly at different urbanization speeds. The fast urbanization cities have a higher probability of vegetation degradation. On the other hand, some fast urbanization cities show high vegetation restoration. This observation reveals that prevention of vegetation degradation in fast urbanization processes is possible.

Table 5. The top 10 cities with the highest variation percentages in the urbanization area and the comparison with outside buffers from 1981 to 2010.

Name	Degradation Inside	Degradation Outside	Slope (Class)	Name	Restoration Inside	Restoration Outside	Slope (Class)
Guangzhou	92.7%	6.9%	0.99 (R)	Berlin	98.0%	98.8%	0.21 (RS)
Istanbul	90.2%	16.5%	0.76 (R)	New York	89.1%	98.4%	0.10 (S)
Tehran	89.5%	69.1%	0.85 (R)	Moscow	86.8%	99.4%	0.58 (RF)
Madrid	88.2%	55.9%	0.54 (RF)	Brussels	86.2%	95.6%	0.08 (S)
Phoenix	87.2%	84.2%	0.77 (R)	Atlanta	84.9%	99.8%	0.53 (RF)
Singapore	86.1%	9.5%	0.71 (RF)	London	84.3%	94.9%	0.02 (S)
Melbourne	82.4%	4.6%	0.38 (MS)	Washington D.C.	82.1%	98.6%	0.20 (RS)
Shanghai	79.6%	14.5%	1.58 (R)	Boston	81.5%	97.2%	0.01 (S)
Kyoto	70.5%	13.8%	0.15 (RS)	Stockholm	78.8%	91.9%	0.00 (S)
Bangkok	70.0%	10.3%	0.85 (R)	Minneapolis	75.8%	99.5%	0.17 (RS)

Table 5. The top 10 cities with the highest variation percentages in the urbanization area and the comparison with outside buffers from 1981 to 2010.

Name	Degradation Inside	Degradation Outside	Slope (Class)	Name	Restoration Inside	Restoration Outside	Slope (Class)
Guangzhou	92.7%	6.9%	0.99 (R)	Berlin	98.0%	98.8%	0.21 (RS)
Istanbul	90.2%	16.5%	0.76 (R)	New York	89.1%	98.4%	0.10 (S)
Tehran	89.5%	69.1%	0.85 (R)	Moscow	86.8%	99.4%	0.58 (RF)
Madrid	88.2%	55.9%	0.54 (RF)	Brussels	86.2%	95.6%	0.08 (S)
Phoenix	87.2%	84.2%	0.77 (R)	Atlanta	84.9%	99.8%	0.53 (RF)
Singapore	86.1%	9.5%	0.71 (RF)	London	84.3%	94.9%	0.02 (S)
Melbourne	82.4%	4.6%	0.38 (MS)	Washington D.C.	82.1%	98.6%	0.20 (RS)
Shanghai	79.6%	14.5%	1.58 (R)	Boston	81.5%	97.2%	0.01 (S)
Kyoto	70.5%	13.8%	0.15 (RS)	Stockholm	78.8%	91.9%	0.00 (S)
Bangkok	70.0%	10.3%	0.85 (R)	Minneapolis	75.8%	99.5%	0.17 (RS)

During the period 2001–2010, with dominant restoration globally observed, the apparent variation is different from the period 1981–2010. From the individual cases (Table 6), most of the variation percentages are higher in the urbanization areas than in the outside buffer. The top 10 restoration cities are not always slow-speed urbanization areas with strong human impact, but the dominant degradation still usually happens in high-speed urbanization cities. On the one hand, the degradation percentage in the urbanization area is consistently higher than the outside buffer in the top 10 degradation cities. On the other hand, in some cities, the restoration percentage can reach more than 95%, but most of the percentage in the outside buffer cannot reach that high level. These phenomena may be explained by the positive human effect on urban vegetation, such as landscape planning for more green space and regular irrigation of the parks. In addition, a large area of cultivated land in some buffers may have a different variation process compared with natural vegetation in other buffers, which results in a complex variation principle.

The result shows clearly that the effect of human impact on vegetation varies greatly. Spatiotemporal difference cannot be ignored when detecting this effect. On temporal characteristics, the degradation trend from the urbanization process is unstable in time series, and the speed of the urbanization process effect on vegetation degradation percentage is difficult to verify. On spatial relations, geographical location may bring consistency to the adjacent metropolis. The adjacent pairs of Berlin and Brussels and New York and Boston both represent similar vegetation variation trends. This finding means that a homogenous climate may strongly influence the vegetation variation direction.

Table 6. The top 10 cities with the highest variation percentages in the urbanization area and the comparison of outside buffers from 2001 to 2010. Dubai is the unique city with mean NDVI < 0.2.

Name	Degradation Inside	Degradation Outside	Slope (Class)	Name	Restoration Inside	Restoration Outside	Slope (Class)
Shanghai	78.5%	38.3%	1.58 (R)	Berlin	100.0%	90.0%	0.21 (RS)
Houston	65.9%	28.8%	0.39 (MS)	Moscow	100.0%	96.1%	0.58 (RF)
Singapore	59.0%	51.5%	0.71 (RF)	Athens	100.0%	88.5%	0.51 (RF)
Buenos Aires	57.7%	53.9%	0.29 (MS)	Johan-nesburg	99.4%	92.8%	0.44 (MS)
Mexico City	52.2%	31.1%	0.34 (MS)	Paris	99.3%	96.0%	0.13 (RS)
Melbourne	39.5%	5.5%	0.38 (MS)	Rio de Janeiro	98.9%	98.7%	0.25 (MS)
Bangkok	38.4%	23.3%	0.85 (R)	Sao Paulo	95.9%	93.7%	0.51 (RF)
St. Louis	36.8%	10.8%	0.20 (RS)	New York	95.3%	91.1%	0.10 (S)
Guangzhou	36.6%	44.0%	0.99 (R)	Detroit	94.7%	97.8%	0.11 (S)
Beijing	35.1%	12.3%	1.14 (R)	Barcelona	93.7%	96.7%	0.29 (MS)

Table 6. The top 10 cities with the highest variation percentages in the urbanization area and the comparison of outside buffers from 2001 to 2010. Dubai is the unique city with mean NDVI < 0.2.

Name	Degradation Inside	Degradation Outside	Slope (Class)	Name	Restoration Inside	Restoration Outside	Slope (Class)
Shanghai	78.5%	38.3%	1.58 (R)	Berlin	100.0%	90.0%	0.21 (RS)
Houston	65.9%	28.8%	0.39 (MS)	Moscow	100.0%	96.1%	0.58 (RF)
Singapore	59.0%	51.5%	0.71 (RF)	Athens	100.0%	88.5%	0.51 (RF)
Buenos Aires	57.7%	53.9%	0.29 (MS)	Johan-nesburg	99.4%	92.8%	0.44 (MS)
Mexico City	52.2%	31.1%	0.34 (MS)	Paris	99.3%	96.0%	0.13 (RS)
Melbourne	39.5%	5.5%	0.38 (MS)	Rio de Janeiro	98.9%	98.7%	0.25 (MS)
Bangkok	38.4%	23.3%	0.85 (R)	Sao Paulo	95.9%	93.7%	0.51 (RF)
St. Louis	36.8%	10.8%	0.20 (RS)	New York	95.3%	91.1%	0.10 (S)
Guangzhou	36.6%	44.0%	0.99 (R)	Detroit	94.7%	97.8%	0.11 (S)
Beijing	35.1%	12.3%	1.14 (R)	Barcelona	93.7%	96.7%	0.29 (MS)

5.2. Confirmation of the Relation by Different NDVI Data Sets

The temporal window and grain size of NDVI may influence the credibility of the result. The conclusion should be verified at different vegetation states. To explore the credibility, MODIS NDVI with 1-km resolution was adopted as a cross-reference data set, and the annual mean values of the pixels have been extracted in Figure 10. Although the values are different between the two data sets of MODIS and VIP, the change tendency in different years is consistent in Beijing’s urbanization area (Figure 10a). With the rapid increase of DN, the ΔMean NDVI calculated by mean NDVI in urbanization areas minus mean NDVI in outside buffers reflects a slight vegetation degradation (Figure 10b). However, the variation difference between Beijing’s urbanization area and the outside buffer is not distinct. The trend of MODIS NDVI or VIP NDVI reflects a slightly increased tendency both in the urbanization area and in the outside buffer during 2000–2010. This phenomenon has two implications. On the one hand, the outside buffer in Beijing was in a dense and highly vegetated state. The difference between the curves of the two data sets is not obvious, which means the vegetation state has been properly considered. On the other hand, on the grain size of 1 km and temporal window of 2000–2010, although the urbanization process does not definitely result in regional vegetation degradation under a vegetation restoration effort, slight vegetation degradation exists when NDVI in the urbanization area is compared with NDVI in the outside buffer (Figure 10b). This result demonstrates that although rapidly urbanized cities have not experienced sharp vegetation degradation, slight vegetation degradation still exists and cannot be ignored.

Figure 10. The vegetation change trends in different statistical objects or data sets for a case in Beijing. (a) Difference between the mean NDVI from two data sets of MODIS and VIP and (b) variation of ΔNDVI (mean NDVI in urbanization area minus mean NDVI in outside buffer) and DN in urbanization area.

Figure 10. The vegetation change trends in different statistical objects or data sets for a case in Beijing. (a) Difference between the mean NDVI from two data sets of MODIS and VIP and (b) variation of ΔNDVI (mean NDVI in urbanization area minus mean NDVI in outside buffer) and DN in urbanization area.

5.3. Urbanization Stage and the Effect on Night Light and Vegetation

The appearance of urbanization can be described as finance aggregation and dwelling expansion, which are both related to the night light variation. However, the relation is not simple. As Elvidge (2014) described, two axes exists, with positive or negative correlation between GDP and population [46]. When the negative correlation appears, the light will shift to uncertainty. Because GDP and population are not in synchronization, and finance aggregation does not strictly rely on the landscape pattern, the direct driving force on light variation in each city is diverse. As Figure 11a shows, if GDP contributes more than population to light growth, the dominant driving force is finance aggregation, or else the dwelling expansion will be the main impact. In developing countries such as China, a large number of rural dwellers’ migrate to cities, resulting in significant urban expansion. On the contrary, in some developed countries, such as the United States, Britain, and Japan, the landscape pattern is relatively stable and finance aggregates without serious land cover changes. If economic depression takes place, the light may synchronously decline without the landscape pattern changing sharply. Thus the driving force is only drawn on the first quadrant. In addition, if GDP increases with a decrease in population, the light variation trend may be more similar to the tendency of GDP. Because electric power consumption at night is an economic indicator, the correlation between GDP and DN has already been confirmed [33,35].

As vegetation activity is more likely to decrease in the circumference of the core urban area, the positive effects generated by the urban environment should be emphasized [12,56]. Liu et al. (2014) separate the mechanism model of vegetation change in built-up areas into an early stage and an advanced stage, labeled, “aggregation” and “expansion,” respectively (Figure 11b) [57]. In the urban expansion stage, the severely changing landscape pattern may result in vegetation biomass decrease, and the impervious surface is increased to support the dwelling space requirement. Nevertheless, although the space requirement also exists in the urban aggregation stage, the negative effect will be countered by the positive demand for urban green space. Consequently, the degradation will convert to restoration with the increased demand for quality in the urban environment. However, the restoration process will not last forever, and the biomass in urban areas will never be as rich as forest. The concept that “urban vegetation variation is the outcome of positive urbanization effect minus negative urbanization effect” is reasonable in describing the relation between urbanization and vegetation degradation [57].

It is important to raise public awareness about how urbanization has an impact on global land-use change [58]. However, the ecological effect triggered by urbanization may be difficult to define, because the results depend on the spatiotemporal scale [59]. The results of our study confirm that, at least for urban biodiversity quantified by vegetation indicator on a macro scale, improving the sustainability and resilience of the urban ecosystem function is possible [60]. Because the urbanization process would not necessarily result in regional vegetation degradation, economic growth may not be adopted as an excuse for urban biodiversity loss. It is worthwhile to pay attention to landscape sustainability and reduce the negative urbanization effect by urban landscape planning.

Figure 11. Urbanization stage and the effect on night light and vegetation: (a) relation between urbanization stage and light [46], and (b) relation between urbanization stage and vegetation [57].

Figure 11. Urbanization stage and the effect on night light and vegetation: (a) relation between urbanization stage and light [46], and (b) relation between urbanization stage and vegetation [57].

6. Conclusions

In this study, we have developed a new approach for delineating the different relations between vegetation degradation and urbanization spatially and temporally through the integrative use of NTL data and NDVI data. The result verifies that the urbanization process would not necessarily result in large-scale vegetation degradation. The correlations between urbanization and vegetation variation are diversified in individual cases. Rapid urbanization cities have a high probability of vegetation degradation. However, although cities in Asia experienced the most rapid urbanization, some cities there did not undergo sharp vegetation degradation. Furthermore, with temporal scale variation, the degradation degree varies, which reflects different urbanization stages and climate history. Consequently, we relate the urbanization effect on night light and vegetation to urbanization stage. In the urban expansion stage, the vegetation biomass decrease may be significant. However, when the urbanization turns to an urban aggregation stage, positive demand of urban green space may increase. Therefore, economic growth may not be adopted as an excuse for urban biodiversity loss. By effective urban landscape planning, urban vegetation restoration would replace the degradation, which may contribute to landscape sustainability.

Our approach has a few limitations. First, the spatial resolution of the VIP data set limits the study to focusing on big metropolises rather than small cities. Second, the limitations from coarse resolution and the spillover effect indicate that the accuracy of urbanization information extracted using NTL data requires further improvement [61,62,63]. The new NPP-VIIRS data, which has overcome the relatively low spatial resolution, the lack of onboard calibration and the pixel saturation issue in DMSP/OLS data, is recommended for obtaining higher-quality urbanization information [64]. Third, the lighting reference area, New York, has been largely stable over time, but that stability does not mean there is no light growth. Fourth, the vegetation variation trend may rely on the NDVI data set, so the trend at the pixel level is, in a sense, uncertain. Nevertheless, the long-term temporal trends of the data set are suitable and unlikely to be replaced by other data sets. The conclusions are credible on the precondition that the result of this article only focuses on the macro scale. Thus, a multi-scale analysis using different data sources will be helpful in reducing uncertainty in future studies.