Responses of Vegetation Cover to Environmental Change in Large Cities of China

Vegetation cover is crucial for the sustainability of urban ecosystems; however, this cover has been undergoing substantial changes in cities. Based on climate data, city statistical data, nighttime light data and the Normalized Difference Vegetation Index (NDVI) dataset, we investigate the spatiotemporal variations of climate factors, urban lands and vegetation cover in 71 large cities of China during 1998–2012, and explore their correlations. A regression model between growing-season NDVI (G-NDVI) and urban land proportion (PU) is built to quantify the impact of urbanization on vegetation cover change. The results indicate that the spatiotemporal variations of temperature, precipitation, PU and G-NDVI are greatly different among the 71 cities which experienced rapid urbanization. The spatial difference of G-NDVI is closely related to diverse climate conditions, while the inter-annual variations of G-NDVI are less sensitive to climate changes. In addition, there is a negative correlation between G-NDVI trend and PU change, indicating vegetation cover in cities have been negatively impacted by urbanization. For most of the inland cities, the urbanization impacts on vegetation cover in urban areas are more severe than in suburban areas. But the opposite occurs in 17 cities mainly located in the coastal areas which have been undergoing the most rapid urbanization. Overall, the impacts of urbanization on G-NDVI change are estimated to be −0.026 per decade in urban areas and −0.015 per decade in suburban areas during 1998–2012. The long-term developments of cities would persist and continue to impact on the environmental change and sustainability. We use a 15-year window here as a case study, which implies the millennia of human effects on the natural biotas and warns us to manage landscapes and preserve ecological environments properly.


Introduction
Over a half of the world's population lives in urban areas which suffer from environmental problems (e.g., air pollution and ecosystem degradation) [1,2].As one of the most important parts in the earth land ecosystem, vegetation provides a wide range of social and environmental services to urban life, and benefits the sustainability of urban ecosystems [2][3][4].For example, urban vegetation can sequestrate carbon [5], regulate microclimates [6], improve air quality [7], preserve biodiversity [8], conserve soil and water and mitigate nature disasters [9].However, palaeoecological and palaeoenvironmental records show that vegetation and ecosystem have been obviously influenced by human disturbances and climate change in long-term landscape evolution [10][11][12].Remarkably, vegetation cover within and around urban areas has experienced obvious transformation during the past decades, which has significantly influenced the sustainability of urban ecosystems [1].

h t t p : / / i r . i s w c . a c . c n
Generally, studies on vegetation variations at macro scales (e.g., city scale and regional scale) are conducted based on the remote sensing technique.Normalized difference vegetation index (NDVI) and vegetation net primary productivity (NPP) derived from remote sensing datasets are main indicators of vegetation activities [13].NPP is the rate of atmospheric carbon uptake by vegetation, which is a useful proxy of vegetation ecosystem services ability [14].Based on MODIS-1 km NPP products, Fu et al. [15] showed a total NPP loss of approximately 167 × 106 g•C from 2001 to 2006 in Guangzhou city, China, while the spatiotemporal pattern of NPP showed obvious variations in local areas.A different approach is to use NDVI to estimate fractional vegetation cover [16], or directly applied to investigate the changes of vegetation cover [17,18].Based on GIMMS NDVI datasets, Eastman et al. [19] indicated that the changing trend of NDVI during 1982-2011 exhibited significant spatial differences over global land surfaces.Differences in spatiotemporal variations of vegetation could be largely attributed to the effects of climate factors [10], anthropogenic activities [20], and their interactions [15,21].
Precipitation and temperature has been demonstrated to be the key climate factors for plant growth and vegetation development in the literatures [22,23].Vegetation cover in humid regions was generally higher than that in arid and semi-arid regions [10,24].But the response of vegetation to precipitation and temperature change were varied in temporal patterns [23,25,26].For example, Zhou et al. [27] demonstrated that temperature change was the leading cause of NDVI decreases, while precipitation played a minor role in high-latitude areas.Sun et al. [28] suggested that increasing precipitation led to an improvement in the vegetation cover, whereas temperature was not a limiting factor in arid and semi-arid regions of China.In cities, anthropogenic activities could impact vegetation cover in a positive or negative way [18,29].On the one hand, large patches of fertile cropland and forest were transformed to urban constructed lands during the process of urbanization [30,31]; on the other hand, people have recognized the importance of vegetation in urban ecosystems, and have given more and more attention to urban greening plans [8,32].However, it is difficult to balance between economic development and ecosystem protection, resulting in obvious constructed land growth and green area reduction [33,34].In addition, the rate of local climate warming would increase if large scale vegetation was removed (e.g., urban heat island) during the process of urbanization [35], which contrarily altered vegetation growing period [36], as well as water and carbon cycles within the urban ecosystem [37].Consequently, local climate change and urban land expansion are the two dominant factors for vegetation variations, even acting as obstructions to reasoned urban planning [8,36,38].
In ancient China, intensive human activities such as deforestation and cultivation have already evidently affected the vegetation and ecosystem [12].In the context of global warming, China has been experiencing a particularly rapid climate change [39][40][41].The impacts of climate change on vegetation cover variation have been frequently explored in previous studies [10,26,28,40], although they mainly focus on the regional scale rather than the city scale.China has experienced remarkably dramatic population growth and rapid urbanization in the past few decades, with the population size of near 1.4 billion in 2012 [29].This implies that the environment (e.g., vegetation cover), particularly in large cities, has been largely impacted by urbanization [42].For example, Li et al. [43] indicated that the rapid urbanization in 1988-2003 led to obvious vegetation cover degradation in the plain areas of Shenzhen city, China.Moreover, urban area in Beijing city nearly doubled during 1997-2002, whereas farmland decreased rapidly as a result of the rapid urbanization [44].Although these studies indicated that urban land expansion significantly influenced the local vegetation variation at individual large cities, the quantitative impacts of urbanization were rarely estimated among different cities.
This study aims to explore the response of vegetation cover to environment change, and assess the impacts of urbanization on vegetation cover change in large cities of China.In this study, we first examine the spatiotemporal changes of climate factors, urban lands and vegetation cover in urban and suburban areas, respectively.Then the correlation analysis among urban land proportion (PU), growing-season NDVI (G-NDVI), temperature (G-T) and precipitation (G-P) is conducted to explain the diversities of G-NDVI change in cities.Finally, the impacts of urbanization on G-NDVI changes are quantified by a regression model.Generally, it would be better to analyse the long-term changes h t t p : / / i r .i s w c . a c .c n of vegetation in cities considering their histories, but we try to use recent observed data (1998-2012) to conduct study in an important "moment" or "window" of human history.It could also reflect long-term changes of the vegetation, landscape, and environment.

Data
City statistical data included urban built-up land area and population of municipal districts, which were often used to represent the size of cities [39].The city statistical data for 1998-2012 were downloaded from China City Statistical Yearbook issued by the National Bureau of Statistics of China (http://www.stats.gov.cn).
NSL dataset from Defense Meteorological Satellite Program (DMSP) Operational Linescan System (OLS) was obtained from the NOAA's National Geophysical Data Center (NGDC) (https://ngdc.noaa.gov/eog/download.html).The spatial resolution and value range of NSL data were 1 km and 0 to 63, respectively.The background noise of the light images was represented by 0. The NSL dataset has been frequently used to extract urban lands with the empirical thresholding technique [30].For example, Liu et al. [45] demonstrated that the threshold for the year 2000 was 49 in Northern Coastal China where many large cities showed prosperous economy (e.g., Beijing and Tianjin).In our study, The NSL data of 1998 and 2012 was adopted to extract urban lands.
SPOT-VGT NDVI datasets for 1998-2012 were obtained from the Image Processing Centre in Vlaamse Instelling voor Technologisch Onderzock (VITO) (Mol, Belgium) (http://www.vgt.vito.be).These data were 10-day composites of atmospherically corrected maximal values with the actual values ranging from −0.1 to 1.In order to eliminate the impact of large bare lands and water areas, grid cells with NDVI value greater than 0.15 were used.Spatial and temporal resolutions of the NDVI dataset are 1 km and 10 day respectively.Monthly NDVI value was synthesized by the maximum value composite (MVC) method based on the 10-day dataset.Given that the NDVI value is impacted by snow, we only focused on the mean NDVI in growing season (G-NDVI).Annual G-NDVI was generated based on the month mean NDVI for April to October in terms of Zhou et al. [27] and Piao et al. [26].
Surface climate data include monthly mean temperature and monthly total precipitation during 1998-2012, which were obtained from the China Meteorological Data Service Center (http://data.cma.cn/site/index.html).Annual growing-season temperature (G-T) and annual growing-season precipitation (G-P) were calculated based on the monthly climate data for April to October.

City Selection and Quantification of Urban Land Expansion
China has a vast territory with complicated geographical and climate conditions, as well as rich diversities in ecosystems.However, the rapid urbanization has profoundly impacted and changed the environmental, social and economic situations of China.In our study, the large cities of China and nearby meteorological stations were selected based on the following two criteria: The urban population in 1998 was larger than 0.5 million; the meteorological station close to city has a continuous climate data record covering the period of 1998-2012.Seventy one cities were eventually selected (Figure 1).In terms of Liu et al. [45], the different extraction thresholds were employed to extract urban lands of selected large cities in China based on the NSL data of 1998 and 2012 (the specific information of extraction thresholds was shown in Figure S1 of Supplementary Material).The geometric center of urban land in 2012 was determined for each city by the Feature to Point Tool in ArcGIS software, which was considered as the spatial city center.A radius was calculated based on the urban built-up land area of 2012 for each city.By means of the Edit Tool in ArcGIS software and based on the city center as well as the radius, a circle for each city was drawn centered on the city center and including the main city area-this was defined as the inner zone.A buffer area with a width of 10 km was generated outside this circle-this was defined as the outer zone.The inside and outside zones were named as urban area (Z1) and suburban area (Z2), respectively (Figure 2).In terms of Liu et al. [45], the different extraction thresholds were employed to extract urban lands of selected large cities in China based on the NSL data of 1998 and 2012 (the specific information of extraction thresholds was shown in Figure S1 of Supplementary Material).The geometric center of urban land in 2012 was determined for each city by the Feature to Point Tool in ArcGIS software, which was considered as the spatial city center.A radius was calculated based on the urban built-up land area of 2012 for each city.By means of the Edit Tool in ArcGIS software and based on the city center as well as the radius, a circle for each city was drawn centered on the city center and including the main city area-this was defined as the inner zone.A buffer area with a width of 10 km was generated outside this circle-this was defined as the outer zone.The inside and outside zones were named as urban area (Z1) and suburban area (Z2), respectively (Figure 2).In terms of Liu et al. [45], the different extraction thresholds were employed to extract urban lands of selected large cities in China based on the NSL data of 1998 and 2012 (the specific information of extraction thresholds was shown in Figure S1 of Supplementary Material).The geometric center of urban land in 2012 was determined for each city by the Feature to Point Tool in ArcGIS software, which was considered as the spatial city center.A radius was calculated based on the urban built-up land area of 2012 for each city.By means of the Edit Tool in ArcGIS software and based on the city center as well as the radius, a circle for each city was drawn centered on the city center and including the main city area-this was defined as the inner zone.A buffer area with a width of 10 km was generated outside this circle-this was defined as the outer zone.The inside and outside zones were named as urban area (Z1) and suburban area (Z2), respectively (Figure 2).

h t t p : / / i r . i s w c . a c . c n
The areas of extracted urban lands in Z1 (AU Z1 ) and Z2 (AU Z2 ) for the 71 large cities were counted by ArcGIS software.The proportion of urban lands in Z1 (PU Z1 ) and Z2 (PU Z2 ) were calculated by the following formula: where PU Z1 indicated the proportion of urban lands in Z1 (%); AU Z1 indicated the area of extracted urban lands in Z1 (km 2 ); A Z1 indicated the total area of Z1 (km 2 ).PU Z2 was calculated based on the area of extracted urban lands in Z2 (AU Z2 ) and the total area of Z2 (A Z2 ).The change of PU between 1998 and 2012 in Z1 (Z2) of each city, namely ∆PU Z1 (∆PU Z2 ), was calculated to represent the rate of urban land expansion during this period (Figure S1 of Supplementary Materials).

Analysis of Spatiotemporal Variations
The average G-NDVI values in Z1 (G-NDVI Z1 ) and Z2 (G-NDVI Z2 ) indicated the vegetation cover in Z1 and Z2 respectively.Because Z1 and Z2 are in the same city, their climate background was thought to be the same.The G-T and G-P were calculated using the meteorological data observed at a nearby meteorological station to represent the climate in growing-season for a given city (including Z1 and Z2).
This study firstly investigated the spatiotemporal variations of G-T, G-P, PU Z1 and PU Z2 among the selected 71 large cities between 1998 and 2012.Then, the spatiotemporal variations of G-NDVI Z1 and G-NDVI Z2 were also analysed with respect to G-T and G-P.The mean annual G-T, G-P and G-NDVI were calculated based on the datasets during 1998-2012.Linear trends in G-T, G-P and G-NDVI during 1998-2012 were examined using ordinary least-squares regression to analyze the inter-annual variations of G-T, G-P and G-NDVI.The trend rates of G-NDVI were used to reflect the directions of vegetation cover change [40].

Analyzing the Sensitivity of G-NDVI to Climate Change
Relationships between NDVI and other variables are generally examined via a correlation analysis [17,20,22,34,40].In this study, the linear regression analysis was undertaken between mean annual G-NDVI and climate factors (G-P, G-T) to explain the spatial variabilities of G-NDVI.In order to assess the sensitivity of G-NDVI to climate change, the correlations between the inter-annual variations of G-NDVI and climate factors (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) and between the trends of G-NDVI and climate factors were examined by Pearson's correlation analysis.

Quantifying the Impact of Urbanization on G-NDVI Change
The differences of trend rates between G-NDVI Z1 and G-NDVI Z2 (trend rate of G-NDVI Z1-Z2 ) were calculated.Because the climate changes in Z1 and Z2 are similar for a given city, the trend rate of G-NDVI Z1-Z2 could be attributed to the different impacts of urban land expansion between Z1 and Z2.In addition, the differences between ∆PU Z1 and ∆PU Z2 (∆PU Z1-Z2 ) were calculated to show the different rates of urban land expansion between Z1 and Z2.The correlations between the trend rates of G-NDVI Z1 (G-NDVI Z2 ) and ∆PU Z1 (∆PU Z2 ) and between the trend rates of G-NDVI Z1-Z2 and ∆PU Z1-Z2 were detected to analyze the impact of urbanization on vegetation cover change.
All the correlation analyses in this study were the single-variant analysis and performed using SPSS version 19.0.The significance of the correlation coefficients was estimated by t-tests at 0.1, 0.05 and 0.01 significance levels.h t t p : / / i r .i s w c . a c .c n

Climate Factors
The climate condition of 71 cities varied from humid-hot climate to dry-cool climate with the latitude increasing (Figure 3a,b).Moreover, the inter-annual variations of G-T and G-P among the selected 71 cities were greatly different (Figure 3c,d).The trend rate of G-T and G-P ranged from −1.12 to 0.74 • C per decade and from −396 to 467 mm per decade respectively (Table 1).The G-T for 34 cities had ascending trends, which were mainly located in Middle Eastern and Middle Southern China (Figure 3c).More than a half of the cities (42 cities) experienced a descending G-P (Figure 3d).They are mainly located in Middle and Southern China.

Climate Factors
The climate condition of 71 cities varied from humid-hot climate to dry-cool climate with the latitude increasing (Figure 3a,b).Moreover, the inter-annual variations of G-T and G-P among the selected 71 cities were greatly different (Figure 3c,d).The trend rate of G-T and G-P ranged from −1.12 to 0.74 °C per decade and from −396 to 467 mm per decade respectively (Table 1).The G-T for 34 cities had ascending trends, which were mainly located in Middle Eastern and Middle Southern China (Figure 3c).More than a half of the cities (42 cities) experienced a descending G-P (Figure 3d).They are mainly located in Middle and Southern China.Table 1.Descriptive statistics of growing-season normalized difference vegetation index (G-NDVI), G-T, G-P and the change of urban land proportion (ΔPU) for the selected 71 large cities.

Mean Annual Value
Change Trend    4).But the ∆PU Z1 and ∆PU Z2 show great difference among the selected 71 cities.The ∆PU Z1 of nine cities is larger than 60%, while two cities have ∆PU Z2 of larger than 60%.The ∆PU Z2 of cities located in northern region and coast areas of China are larger than in southwestern China.Generally, the ∆PU Z1 of cities located in the inland of China were larger than ∆PU Z2 (i.e., ∆PU Z1-Z2 > 0), while the ∆PU Z1 of some cities located in the coastal areas were smaller than ∆PU Z2 (i.e., ∆PU Z1-Z2 < 0).Overall, the mean values of ∆PU Z1 and ∆PU Z2 are 37% and 22% respectively (Table 1), implying that urban areas of 71 cities have experienced faster urbanization than suburban areas.
decade for G-T and mm per decade for G-P; 5 and 6 the change of urban land proportion between 1998 and 2012 in urban and suburban area respectively, in %.

Urban Lands
PUZ1 and PUZ2 obviously increased from 1998 to 2012, indicating that the selected 71 cities have experienced rapid urbanization (Figure 4).But the ΔPUZ1 and ΔPUZ2 show great difference among the selected 71 cities.The ΔPUZ1 of nine cities is larger than 60%, while two cities have ΔPUZ2 of larger than 60%.The ΔPUZ2 of cities located in northern region and coast areas of China are larger than in southwestern China.Generally, the ΔPUZ1 of cities located in the inland of China were larger than ΔPUZ2 (i.e., ΔPUZ1-Z2 > 0), while the ΔPUZ1 of some cities located in the coastal areas were smaller than ΔPUZ2 (i.e., ΔPUZ1-Z2 < 0).Overall, the mean values of ΔPUZ1 and ΔPUZ2 are 37% and 22% respectively (Table 1), implying that urban areas of 71 cities have experienced faster urbanization than suburban areas.

G-NDVI
Spatial distributions of mean annual growing-season NDVI (including G-NDVIZ1 and G-NDVIZ2) for the 71 cities were remarkably different (Figure 5a,b).G-NDVIZ1 and G-NDVIZ2 ranged from 0.23 to 0.52 and 0.29 to 0.64, with standard deviations of them being 0.06 and 0.08 respectively (Table 1).The G-NDVIZ1 of 22 cities was lower than 0.3, while one cities had G-NDVIZ1 of greater than 0.5.Two cities had G-NDVIZ2 of lower than 0.3, while 34 cities had G-NDVIZ2 of greater than 0.5.G-NDVIZ2 was lower in northern region and coast areas of China than in southern China.For a given city, the mean annual G-NDVIZ1 was generally lower than G-NDVIZ2.The mean values of G-NDVIZ1 and G-NDVIZ2 were 0.34 and 0.48 respectively.Moreover, some adjacent cities with similar climate conditions had different G-NDVI.For example, G-NDVIZ2 of Shanghai and Nantong were 0.36 and 0.51 respectively.for the 71 cities were remarkably different (Figure 5a,b).G-NDVI Z1 and G-NDVI Z2 ranged from 0.23 to 0.52 and 0.29 to 0.64, with standard deviations of them being 0.06 and 0.08 respectively (Table 1).The G-NDVI Z1 of 22 cities was lower than 0.3, while one cities had G-NDVI Z1 of greater than 0.5.Two cities had G-NDVI Z2 of lower than 0.3, while 34 cities had G-NDVI Z2 of greater than 0.5.G-NDVI Z2 was lower in northern region and coast areas of China than in southern China.For a given city, the mean annual G-NDVI Z1 was generally lower than G-NDVI Z2 .The mean values of G-NDVI Z1 and G-NDVI Z2 were 0.34 and 0.48 respectively.Moreover, some adjacent cities with similar climate conditions had different G-NDVI.For example, G-NDVI Z2 of Shanghai and Nantong were 0.36 and 0.51 respectively.
The inter-annual variations of G-NDVI during 1998-2012 also showed great difference among the 71 selected large cities (Figure 5c,d).The G-NDVI Z1 for 40 cities experienced a decreasing trend (Figure 5c).The trend rates of G-NDVI Z1 for 3 cities were greater than 0.05 per decade.Contrarily, the G-NDVI Z2 for most of the cities (51 cities) experienced an increasing trend (Figure 5d).The trend rates of G-NDVI Z2 for 13 cities were greater than 0.05 per decade.Cities with decreasing G-NDVI (particularly decreasing G-NDVI Z2 ) were mainly located in the northern region and eastern coastal areas, where rapid urbanization and economic development occurred during 1998-2012.The trend rate of G-NDVI Z1 and G-NDVI Z2 ranged from −0.099 to 0.067 per decade and from −0.082 to 0.092 per decade, respectively (Table 1).Overall, the average trend rates of G-NDVI Z1 and G-NDVI Z2 were −0.007 and 0.016 per decade, respectively.The inter-annual variations of G-NDVI during 1998-2012 also showed great difference among the 71 selected large cities (Figure 5c,d).The G-NDVIZ1 for 40 cities experienced a decreasing trend (Figure 5c).The trend rates of G-NDVIZ1 for 3 cities were greater than 0.05 per decade.Contrarily, the G-NDVIZ2 for most of the cities (51 cities) experienced an increasing trend (Figure 5d).The trend rates of G-NDVIZ2 for 13 cities were greater than 0.05 per decade.Cities with decreasing G-NDVI (particularly decreasing G-NDVIZ2) were mainly located in the northern region and eastern coastal areas, where rapid urbanization and economic development occurred during 1998-2012.The trend rate of G-NDVIZ1 and G-NDVIZ2 ranged from −0.099 to 0.067 per decade and from −0.082 to 0.092 per decade, respectively (Table 1).Overall, the average trend rates of G-NDVIZ1 and G-NDVIZ2 were −0.007 and 0.016 per decade, respectively.
A great difference in the changing trend of G-NDVIZ1 and G-NDVIZ2 was found for a given city (Figure 5c,d).For example, the trend rates of G-NDVIZ1 and G-NDVIZ2 for Beijing city were 0.026 and −0.027 per decade, respectively.Moreover, some adjacent cities with similar climate change conditions also showed great differences in the G-NDVI change.For instance, the trend rates of G-NDVIZ1 for Shanghai and its adjacent city (Nantong) were 0.001 and −0.098 per decade respectively.

Relationship between Mean Annual G-NDVI and Climate
The G-NDVI was lower when the G-T was below 19 °C and the G-P was below 500 mm (Figure 6).G-NDVI generally increased with G-T and G-P.However, the G-NDVI value was relatively low when G-T exceeded 26 °C and G-P exceeded 1200 mm.High values of G-NDVIZ1 occurred when G-T ranged between 23 and 24.5 °C and G-P ranged between 750 and 1000 mm (Figure 6a).High values of G-NDVIZ2 occurred when G-T ranged between 23 and 25 °C and G-P ranged between 750 and 1200 mm (Figure 6b).A great difference in the changing trend of G-NDVI Z1 and G-NDVI Z2 was found for a given city (Figure 5c,d).For example, the trend rates of G-NDVI Z1 and G-NDVI Z2 for Beijing city were 0.026 and −0.027 per decade, respectively.Moreover, some adjacent cities with similar climate change conditions also showed great differences in the G-NDVI change.For instance, the trend rates of G-NDVI Z1 for Shanghai and its adjacent city (Nantong) were 0.001 and −0.098 per decade respectively.

Relationship between Mean Annual G-NDVI and Climate
The G-NDVI was lower when the G-T was below 19 • C and the G-P was below 500 mm (Figure 6).G-NDVI generally increased with G-T and G-P.However, the G-NDVI value was relatively low when G-T exceeded 26 • C and G-P exceeded 1200 mm.High values of G-NDVI Z1 occurred when G-T ranged between 23 and 24.5 • C and G-P ranged between 750 and 1000 mm (Figure 6a).High values of G-NDVI Z2 occurred when G-T ranged between 23 and 25 • C and G-P ranged between 750 and 1200 mm (Figure 6b).The result of linear regression analysis showed that G-NDVI was positively correlated with climate variables (i.e., G-T and G-P) (Figure 7).The correlation coefficients between G-NDVI and climate variables were significant at the 0.01 significance level.Overall, the values of mean annual G-NDVI are closely related to climate factors.The mean annual G-NDVI of the selected 71 cities increased from areas with dry-cool climate to those with humid-hot climate.

Relationship between G-NDVI Variations and Climate Change
For the G-T, the correlation for the inter-annual variation between G-T and G-NDVI was insignificant for most of the selected cities (Figure 8a,c).For the G-P, cities with significantly positive relationships between G-NDVI and G-P were mainly located in Northern China where the mean The result of linear regression analysis showed that G-NDVI was positively correlated with climate variables (i.e., G-T and G-P) (Figure 7).The correlation coefficients between G-NDVI and climate variables were significant at the 0.01 significance level.Overall, the values of mean annual G-NDVI are closely related to climate factors.The mean annual G-NDVI of the selected 71 cities increased from areas with dry-cool climate to those with humid-hot climate.The result of linear regression analysis showed that G-NDVI was positively correlated with climate variables (i.e., G-T and G-P) (Figure 7).The correlation coefficients between G-NDVI and climate variables were significant at the 0.01 significance level.Overall, the values of mean annual G-NDVI are closely related to climate factors.The mean annual G-NDVI of the selected 71 cities increased from areas with dry-cool climate to those with humid-hot climate.

Relationship between G-NDVI Variations and Climate Change
For the G-T, the correlation for the inter-annual variation between G-T and G-NDVI was insignificant for most of the selected cities (Figure 8a,c).For the G-P, cities with significantly positive relationships between G-NDVI and G-P were mainly located in Northern China where the mean

Relationship between G-NDVI Variations and Climate Change
For the G-T, the correlation for the inter-annual variation between G-T and G-NDVI was insignificant for most of the selected cities (Figure 8a,c).For the G-P, cities with significantly positive relationships between G-NDVI and G-P were mainly located in Northern China where the mean annual precipitation was lower than 500 mm (Figure 8b,d).However, the correlation between the inter-annual variations of G-P and G-NDVI was insignificant for most of the cities (Figure 8).Moreover, the results of Pearson correlation test indicate that the correlations between the trend rates of G-NDVI and G-T and between the trend rates of G-NDVI and G-P are insignificant (Table 2).Overall, the inter-annual variations of G-NDVI were less sensitive to the changes of G-P and G-T.
Sustainability 2018, 10, 270 10 of 16 annual precipitation was lower than 500 mm (Figure 8b,d).However, the correlation between the inter-annual variations of G-P and G-NDVI was insignificant for most of the cities (Figure 8).Moreover, the results of Pearson correlation test indicate that the correlations between the trend rates of G-NDVI and G-T and between the trend rates of G-NDVI and G-P are insignificant (Table 2).
Overall, the inter-annual variations of G-NDVI were less sensitive to the changes of G-P and G-T.

Impact of Urbanization on G-NDVI Change
Table 2 shows that the negative correlations between ΔPUZ1 and the trend of G-NDVIZ1 and between ΔPUZ2 and the trend of G-NDVIZ2 are significant (P < 0.01), indicating that urban land expansion is closely associated with vegetation cover change in 71 large cities.In addition, the trends of G-NDVIZ1-Z2 were negative for 54 cities mainly located in the inland of China, implying that the trend of G-NDVI in urban area was generally lower than that in suburban area for a given inland city (Figure 9).But the opposites were found in 17 cities mainly located in the coastal areas of China.The smallest trend of G-NDVIZ1-Z2 was found in Xinyang city with the rate being −0.093 per decade.The largest trend of G-NDVIZ1-Z2 was found in Shanghai city with the rate being 0.083 per decade.Overall,

Impact of Urbanization on G-NDVI Change
Table 2 shows that the negative correlations between ∆PU Z1 and the trend of G-NDVI Z1 and between ∆PU Z2 and the trend of G-NDVI Z2 are significant (P < 0.01), indicating that urban land expansion is closely associated with vegetation cover change in 71 large cities.In addition, the trends of G-NDVI Z1-Z2 were negative for 54 cities mainly located in the inland of China, implying that the trend of G-NDVI in urban area was generally lower than that in suburban area for a given inland city (Figure 9).But the opposites were found in 17 cities mainly located in the coastal areas of China.

h t t p : / / i r . i s w c . a c . c n
The smallest trend of G-NDVI Z1-Z2 was found in Xinyang city with the rate being −0.093 per decade.The largest trend of G-NDVI Z1-Z2 was found in Shanghai city with the rate being 0.083 per decade.Overall, vegetation cover in suburban areas has been greatly influenced by urbanization for cities located in the coastal area of China.But vegetation cover in urban areas has experienced more severe urbanization impact than in suburban areas for most of the inland cities.  Figure 10 shows that the trend of G-NDVIZ1-Z2 is negatively correlated with ΔPUZ1-Z2 (R 2 = 0.3626, P < 0.01).The regression equation between the trend of G-NDVIZ1-Z2 and ΔPUZ1-Z2 is y = −0.0007x− 0.0012.Through substituting the ΔPUZ1 and ΔPUZ2 of the 71 cities into this regression equation respectively (the constant being zero), the trends of G-NDVI induced by urbanization were calculated.We found that G-NDVI changes have been negatively impacted by urbanization in urban and suburban areas for most of the selected cities (Figure S1 of Supplementary Materials).On average, the impact of urbanization on G-NDVI change was estimated to be −0.026 per decade in Z1 and −0.015 per decade in Z2 during 1998-2012.Figure 10 shows that the trend of G-NDVIZ1-Z2 is negatively correlated with ΔPUZ1-Z2 (R 2 = 0.3626, P < 0.01).The regression equation between the trend of G-NDVIZ1-Z2 and ΔPUZ1-Z2 is y = −0.0007x− 0.0012.Through substituting the ΔPUZ1 and ΔPUZ2 of the 71 cities into this regression equation respectively (the constant being zero), the trends of G-NDVI induced by urbanization were calculated.We found that G-NDVI changes have been negatively impacted by urbanization in urban and suburban areas for most of the selected cities (Figure S1 of Supplementary Materials).On average, the impact of urbanization on G-NDVI change was estimated to be −0.026 per decade in Z1 and −0.015 per decade in Z2 during 1998-2012.The urbanization is an increasing active driving force of landscape change worldwide, and their long-term impacts on environmental change and sustainability could be described through concrete evidences within a relative short period.This study indicated that the mean annual G-NDVI were remarkably different among the selected 71 cities of China.This may be because the difference of NDVI values in China's cities is related to not only the climate conditions but also the development and management of cities (e.g., urbanization level) [34].For larger cities, urban areas include a large proportion of constructed lands and a high density of building [46], resulting in relative lower NDVI in urban areas [34].In addition, larger cities generally have more satellite towns and vaster intercity highway network around cities, implying that vegetation cover in suburban areas of larger city would be lower than that in smaller cities [47,48].This is confirmed by the results presented in Figure 5, where the values of G-NDVI Z2 in developed coast areas are generally lower than that in less developed inland areas.Despite this, the result of linear regression analysis showed a positive correlation between mean annual G-NDVI and climate variables (Figure 7).Consequently, the spatial difference of mean annual G-NDVI among the selected 71 cities is closely related to climate variabilities, but city development may play a certain role in local scales.

Driving Forces of the Temporal Variability of G-NDVI
The change of G-NDVI in the selected 71 cities was found to be less sensitive to climate change.Only few cities located in Northern China with the mean annual G-P being lower than 500 mm showed significant positive relationships between G-NDVI and G-P.Similar phenomenon was found by previous studies [25,26].For example, there was a strong sensitivity of vegetation to precipitation in the regions where the annual rainfall ranged from 100 to 400 mm in Central Asia [23].In general, the weak responses of NDVI to temperature change were attributed to increased water evaporation [26] and different vegetation types [10,24].While the lower sensitivity of vegetation to precipitation in humid regions were ascribed to the limited water-use capacity of vegetation [49], the lower temperature and sunshine levels as precipitation increased [26], as well as the uncertainties of data adopted (e.g., cloud effects in the NDVI time-series) [50].
In cities, urbanization may be the crucial driver of effects on vegetation cover change and can be greater than climate factors [20,34] and hence mask the climate effects.For instance, urban land expansion, driven by urbanization in cities, could directly transform a large amount of agricultural lands to urban lands [43,44,47].As a result, the effects of urban land expansion are able to disrupt the coupling between vegetation and precipitation [38].The difference of inter-annual variation of G-NDVI between Z1 and Z2 for a given city also confirmed the impacts of urban land expansion in this study (Figure 9).The trend rate of G-NDVI change in some adjacent cities also showed a great difference.Climate conditions in a given city and among adjacent cities are similar, the difference was deemed to result from urban land expansion rather than climate variations particularly given that the selected 71 cities underwent rapid urbanization during the study period.
Our results indicated that for most of the selected cities changes in the G-NDVI were negatively impacted by urban land expansion in urban and suburban areas.Similar results were reported by previous studies based on different methods.For example, an analysis of NPP dynamics conducted by Peng et al. [18] showed that urbanization resulted in a lasting and observable loss of NPP over time and space.From the perspective of landscape patterns, Xia et al. [44] suggested that the urban area in Beijing city nearly doubled during 1997-2002, and farmland decreased rapidly as a result of urbanization.Using the correlation analysis, Sun et al. [34] found that NDVI change trends of China metropolises were negatively related to the change trends of urban area, population and GDP.While other studies qualitatively analyzed the impacts of urbanization on vegetation cover change, we quantified it based on a regression model between G-NDVI and PU.This method has excluded the h t t p : / / i r .i s w c . a c .c n impact of climate change on G-NDVI through calculating the difference of G-NDVI trend rate between Z1 and Z2 in a given city.
Certainly, the method presented here also has its limitations.One such example is that we have to rely on the assumption that the trend rate of G-NDVI Z1-Z2 is mainly caused by urban land expansion rather than other anthropogenic activities (i.e., agricultural activities and greening measures) [18].Because our sample size is big enough relative to the influences of these minor factors, we feel that this assumption is pretty reasonable.Another improvement based on the more specific land cover change will be conducted in further studies.

Conclusions
In this study, the spatiotemporal variations of vegetation cover and its sensitivity to climate change was investigated based on the growing-season NDVI, temperature and precipitation in 71 large cities of China during 1998-2012.The impact of urban expansion on vegetation cover change in urban and suburban areas was also quantitatively assessed.The main findings are as follows: (1) The mean annual G-T, G-P and G-NDVI of the selected 71 cities are found to be greatly different.
The spatial difference of G-NDVI is closely related to diverse climate conditions.Overall, the mean annual G-NDVI of 71 cities increases from dry-cool climate to humid-hot climate.Cities generally have a long history and would keep on developing.However, many modern cities have encountered obstacles to maintain urban ecosystem sustainability which caused by intensive human activities (e.g., destroying vegetation).This study not only implies a significant long-term impact of human activities on the landscape, but also provides us information about human effects on the natural biotas.Therefore, we need to properly use these knowledges to plan the landscape management and ecological environment conservation for the future.

Figure 1 .
Figure 1.Locations of the selected 71 large cities in China and the corresponding urban populations in 1998.

Figure 2 .
Figure 2.An example of the extracted urban lands of 1998 and 2012 in the urban area (Z1) and suburban area (Z2) of Nanchang city.

Figure 1 .
Figure 1.Locations of the selected 71 large cities in China and the corresponding urban populations in 1998.

Sustainability 2018, 10 , 270 4 of 16 Figure 1 .
Figure 1.Locations of the selected 71 large cities in China and the corresponding urban populations in 1998.

Figure 2 .
Figure 2.An example of the extracted urban lands of 1998 and 2012 in the urban area (Z1) and suburban area (Z2) of Nanchang city.

Figure 2 .
Figure 2.An example of the extracted urban lands of 1998 and 2012 in the urban area (Z1) and suburban area (Z2) of Nanchang city.

Figure 3 .
Figure 3. Spatial distribution of mean annual growing-season temperature (G-T) (a) and growingseason precipitation (G-P) (b) and inter-annual variations of G-T (c) and G-P (d) during 1998-2012 for the selected 71 large cities.

Figure 3 .
Figure 3. Spatial distribution of mean annual growing-season temperature (G-T) (a) and growing-season precipitation (G-P) (b) and inter-annual variations of G-T (c) and G-P (d) during 1998-2012 for the selected 71 large cities.

Figure 4 .
Figure 4. Spatial distribution of the ΔPUZ1 (a) and ΔPUZ2 (b) during 1998-2012 for the selected 71 large cities.

Figure 5 .
Figure 5. Spatial distribution of mean annual G-NDVIZ1 (a) and G-NDVIZ2 (b) and inter-annual variations of G-NDVIZ1 (c) and G-NDVIZ2 (d) during 1998-2012 in the selected 71 large cities.

Figure 5 .
Figure 5. Spatial distribution of mean annual G-NDVI Z1 (a) and G-NDVI Z2 (b) and inter-annual variations of G-NDVI Z1 (c) and G-NDVI Z2 (d) during 1998-2012 in the selected 71 large cities.

Figure 6 .
Figure 6.The relationship of mean annual G-T, G-P and G-NDVIZ1 (a) G-NDVIZ2 (b) of the selected 71 large cities.Each dot in the panel corresponds to three values (i.e., the mean annual G-NDVI, G-T and G-P).

Figure 7 .
Figure 7. Correlations between the mean annual G-T and G-NDVI (a) and between the mean annual G-P and G-NDVI (b) of the selected 71 large cities.

Figure 6 .
Figure 6.The relationship of mean annual G-T, G-P and G-NDVI Z1 (a) G-NDVI Z2 (b) of the selected 71 large cities.Each dot in the panel corresponds to three values (i.e., the mean annual G-NDVI, G-T and G-P).

Sustainability 2018, 10 , 270 9 of 16 Figure 6 .
Figure 6.The relationship of mean annual G-T, G-P and G-NDVIZ1 (a) G-NDVIZ2 (b) of the selected 71 large cities.Each dot in the panel corresponds to three values (i.e., the mean annual G-NDVI, G-T and G-P).

Figure 7 .
Figure 7. Correlations between the mean annual G-T and G-NDVI (a) and between the mean annual G-P and G-NDVI (b) of the selected 71 large cities.

Figure 7 .
Figure 7. Correlations between the mean annual G-T and G-NDVI (a) and between the mean annual G-P and G-NDVI (b) of the selected 71 large cities.

Figure 8 .
Figure 8. Distributions of the correction coefficients for the inter-annual variation between G-NDVIZ1 and G-T (a), G-NDVIZ1 and G-P (b), G-NDVIZ2 and G-T (c), and G-NDVIZ2 and G-P (d) in the selected 71 large cities during 1998-2012.

ΔPUZ1 1 ΔPUZ2 2
Trend rate of G-T3 Trend rate of G-P4  Trend rate of G-NDVIZ1 5 −0.37 *the change of urban land proportion between 1998 and 2012 in urban and suburban area respectively; 3 and 4 growing-season temperature and precipitation respectively; 5 and 6 change trends of growing-season NDVI in urban and suburban area respectively; ** means the correlation is significant at the 0.01 significance level.

Figure 8 .Table 2 .
Figure 8. Distributions of the correction coefficients for the inter-annual variation between G-NDVI Z1 and G-T (a), G-NDVI Z1 and G-P (b), G-NDVI Z2 and G-T (c), and G-NDVI Z2 and G-P (d) in the selected 71 large cities during 1998-2012.

Sustainability 2018 ,
10, 270 11 of 16vegetation cover in suburban areas has been greatly influenced by urbanization for cities located in the coastal area of China.But vegetation cover in urban areas has experienced more severe urbanization impact than in suburban areas for most of the inland cities.

Figure 10 .
Figure 10.The relationship between trends of G-NDVIZ1-Z2 and ΔPUZ1-Z2 for the selected 71 large cities.

Figure 10
Figure10shows that the trend of G-NDVI Z1-Z2 is negatively correlated with ∆PU Z1-Z2 (R 2 = 0.3626, P < 0.01).The regression equation between the trend of G-NDVI Z1-Z2 and ∆PU Z1-Z2 is y = −0.0007x− 0.0012.Through substituting the ∆PU Z1 and ∆PU Z2 of the 71 cities into this regression equation respectively (the constant being zero), the trends of G-NDVI induced by urbanization were calculated.We found that G-NDVI changes have been negatively impacted by urbanization in urban and suburban areas for most of the selected cities (FigureS1of Supplementary Materials).On average, the impact of urbanization on G-NDVI change was estimated to be −0.026 per decade in Z1 and −0.015 per decade in Z2 during 1998-2012.

Figure 10 .
Figure 10.The relationship between trends of G-NDVIZ1-Z2 and ΔPUZ1-Z2 for the selected 71 large cities.

( 2 )
The changes of G-T, G-P, PU and G-NDVI during 1998-2012 are different among the selected 71 cities.The mean values of ∆PU Z1 and ∆PU Z2 were 37% and 22% respectively, indicating that the selected 71 cities have experienced rapid urbanization during 1998-2012.The trend rates of G-NDVI Z1 and G-NDVI Z2 range from −0.099 to 0.067 per decade and −0.082 to 0.092 per decade respectively.G-NDVI changes are less sensitive to climate change, while closely related to urban land expansion.There is a negative correlation between G-NDVI trend and PU change, indicating vegetation cover in cities has been negatively impacted by urbanization.(3)For most of the inland cities, vegetation cover in urban areas has experienced more severe urbanization impact than in suburban areas.But opposites occur in the 17 cities mainly located in the coastal areas of China.The average impacts of urbanization on G-NDVI change were estimated to be −0.026 per decade in Z1 and −0.015 per decade in Z2 during 1998-2012.

Table 1 .
Descriptive .1.2.Urban Lands PU Z1 and PU Z2 obviously increased from 1998 to 2012, indicating that the selected 71 cities have experienced rapid urbanization (Figure statistics of growing-season normalized difference vegetation index (G-NDVI), G-T, G-P and the change of urban land proportion (∆PU) for the selected 71 large cities. 1 and 2 mean the growing-season NDVI of urban area and suburban area, respectively; 3 growing-season temperature, in • C; 4 growing-season precipitation, in mm; The unit of change trend is • C per decade for G-T and mm per decade for G-P; 5 and 6 the change of urban land proportion between 1998 and 2012 in urban and suburban area respectively, in %. h t t p : / / i r .i s w c . a c .c n 3

Table 2 .
The results of Pearson correlation test for the selected 71 large cities.