Exploring the Combined Effect of Urbanization and Climate Variability on Urban Vegetation: A Multi-Perspective Study Based on More than 3000 Cities in China

Abstract

More than 3000 cities in China were used to study the effect of urbanization and local climate variability on urban vegetation across different geographical and urbanization conditions. The national scale estimation shows that China’s urban vegetation depicts a trend of degradation from 2000 to 2015, especially in developed areas such as the Yangtze River Delta. According to the panel models, the increase of precipitation (PREC), solar radiation (SRAD), air temperature (TEMP), and specific humidity (SHUM) all enhance urban vegetation, while nighttime light intensity (NLI), population density (POPDEN), and fractal dimension (FRAC) do the opposite. The effects change along the East–West gradient; the influences of PREC and SHUM become greater, while those of TEMP, SRAD, NLI, AREA, and FRAC become smaller. PREC, SHUM, and SRAD play the most important roles in Northeast, Central, and North China, respectively. The role of FRAC and NLI in East China is much greater than in other regions. POPDEN remains influential across all altitudes, while FRAC affects only low-altitude cities. NLI plays a greater role in larger cities, while FRAC and POPDEN are the opposite. In cities outside of the five major urban agglomerations, PREC has a great influence while the key factors are more diversified inside.

1. Introduction

Urban vegetation is a key regulator for social, ecological, and economic processes of cities, modulating their effects on public health, social equity, and the ecological environment [1,2,3]. Urban vegetation in forms of urban parks and forests provides respite and recreation for city residents and is therefore related to multiple quality-of-life indicators [4], socialization, and reduced stress [5,6]. Furthermore, urban vegetation regulates ecosystem services, such as air and water purification, wind and noise reduction, and microclimate regulation as well as social, cultural, and psychological values through its control of biogeochemical cycling [7,8,9,10]. In addition, as vegetation cover mitigates the damage to ecosystem function posed by excessive impervious surfaces, urban vegetation can be interpreted as a trade-off between the amount of vegetation and impervious surfaces, the composition of the urban landscape, and the health status of an urban ecosystem [11]. Therefore, it is of planning significance and ecological significance to study the macro dynamics of urban vegetation and explore its relationship with various influencing factors. In recent decades, the world has experienced unprecedented urban growth and land consumption; the pressures on vital ecosystem functions have escalated rapidly [12,13]. By 2030, approximately five billion people will live in cities [14]. In this context of global urbanization combined with climate variability, understanding the response of urban vegetation to various factors and using it as a basis to protect urban ecosystems will help to achieve the goal of sustainable cities in the United Nations’ sustainability framework [15,16,17].

The interaction between various human factors and climate factors in the city has a complex impact on vegetation, which makes the driving factors of urban vegetation more heterogeneous than those in natural areas [18]. For example, a change of land cover to impervious surfaces may result in eco-environmental threats with net primary production reduction, surface temperature variation, and, thus, vegetation degradation [19,20]. However, as evidenced by recent studies, urbanization factors can also enhance urban vegetation activity [21,22,23]. The plant growth in urban areas might be promoted by warmer temperatures and a greater tropospheric CO₂ concentration [24,25]. Moreover, because the advocation and the practice of promoting the value of urban ecosystem services are rooted in landscape management [26], green infrastructure investments are supported, and some negative effects are mitigated with the urbanization process [27,28,29]. As a result, the interconnection between vegetation degradation and urbanization comes to be a complex and nonlinear system [19,23]. The spatiotemporal relations between urbanization and vegetation degradation may be diversified and related to the stage of urbanization level or geographical location [15,30,31]. However, a systematic evaluation of the key influence factors of urban vegetation covering multiple cities over large areas is still lacking [18,23,32]. In addition, urban form factors are seldomly taken into account in urban vegetation research. It is suggested that urban sprawl leads to greater air pollution, which is related to a larger numbers of commuters owing to functional and spatial separation of places for living and working [33,34,35] and reduced green space coupled with loss of habitats and ecosystem fragmentation [36]. All of these effects have a direct or indirect impact on urban vegetation.

The primary goal of this paper was to study the heterogeneity of the impact of climatic and urbanization factors on urban vegetation based on national scale urban samples in China. We adopted the remote sensed Index (NDVI) as an effective index for describing the dynamics of urban vegetation [37,38], which can indirectly estimate gross and net primary productivity, biomass, and green leaf area in a variety of ecosystems [39,40,41,42,43]. We used more than 3000 refined cities nationwide as research units to carry out this empirical study. The cities included all large-, medium- and small-sized ones across China. We used the data from every five years between 2000 and 2015 to build fixed effect panel models, which effectively took into consideration the impact of individual effects and provided accurate estimations of each impact factor. We modeled the urban NDVI from multiple perspectives, including geographical perspectives, altitude perspectives, urban scale perspectives, and urban agglomeration or non-agglomeration perspectives. Then, the geographic variations and the urbanization gradients of the effects of four climate factors (precipitation, humidity, solar radiation, temperature) as well as four urban characteristic factors (urban area scale, population density, night light intensity, shape spread index) on NDVI were evaluated from each perspective. Finally, the reasons for the high heterogeneity of the impact on urban vegetation were also discussed. The accurate knowledge of urbanization and climate variability effect on urban vegetation not only helps enhance understanding the urban ecosystem responses to local or global environment change but is also essential for formulating mitigation strategies in cities [44].

2. Materials and Methods

2.1. Identify the Spatial Range of Cities

Our study focused on the city areas in China, which refer to areas of the urban land according to the land use classification system smaller than a city’s administrative boundary but larger than its impervious surface (Figure 1a). Cities are traditionally defined based on administrative boundaries, such as municipalities, prefectural cities, and country-level cities in China. However, modern city systems are complex and consist of strongly interrelated components [45,46,47,48]. In addition, urban form has been dramatically changed during the past decades due to rapid urbanization globally [46,49,50], especially in China. There is an apparent inconsistency between traditional administrative boundaries of cities and the real densely populated or functional central urban extent. In some areas of China, multiple urban units are connected into a mega city across administrative regions, while other cities are leaping to form new independent cities far away from the core urban areas. Because we need to calculate the geometric shape index of each city, such as the fractal dimension of the boundary, it is key to effectively identify the mega city across the administrative region and the multi-cities within the same administrative region.

In this paper, we define a city as an independent urban patch closely connected or adjacent in space and closely related in socioeconomic function. Based on the high-resolution land use data and the recognized spatial ranges of natural cities according to the previous study [51], we extracted the spatial range of cities. According to Song et al., natural cities are redefined as real spatially, closely connected urban central regions where population is dense, human activities are abundant and active, and built environment and infrastructure conditions are relatively complete. The process of natural city identification includes the following steps. The first step is to calculate point of interest (POI) density and generate a POI density map with Kernel density functions with the spatial resolution of 500 m. The second step is to determine the threshold to classify the urban regions according to Kernel and POI density maps. The final step is to convert the raster map to polygons to derive the initial city boundaries.

However, the spatial resolution of the identified city, which influences the accurate estimation of the fractal dimensions of the city boundary, is not enough. In addition, due to the discarding process, the spatial range of urban land covered by the identified natural cities is not complete. Therefore, we combined the spatial range of identified natural cities with the high-resolution urban land use map to get the detailed spatial boundary of our study units. Firstly, we used the independent urban patches obtained from the 30 m-resolution land-use map as the candidate objects of urban units. Then, we defined the independent land use patches of multiple cities with overlapping parts with one natural city as the parts of one city, which could effectively avoid the problem that some cities are separated by rivers or that cities with a certain distance in space but closely related in function or population distribution are identified as multiple city units. Finally, we regarded urban patches derived from the urban land map that do not overlap with the spatial range of natural cities as additional urban units. Through the above processes, 3189 city units were finally obtained (Figure 1b).

The urban categories derived from 30 m land use data of 2000, 2005, 2010, and 2015 were obtained from the Resource and Environment Data Cloud Platform, CAS (http://www.resdc.cn/DataList.aspx), and the corresponding spatial ranges of identified cities for the four periods were obtained accordingly (Figure 1c–e).

2.2. Metric of Urban Vegetation Cover

For the urban vegetation coverage, the normalized difference vegetation index (NDVI) was derived from MODIS dataset MOD13Q1 with 250 m spatial resolution. In order to effectively match with other annual scale data and to take into account the vegetation growth duration, we used the averaged 16 days NDVI throughout the whole year to represent the overall greening level of each natural city. The MOD13Q1 dataset was used because the MODIS satellite has a larger swath width, which can effectively meet the needs of region-scale cross-section data analysis and modeling [52]. The cloud-contaminated area was removed in data preprocessing. The NDVI is a normalized transform of the near-infrared radiation (NIR) to red reflectance (Red) ratio, ρNIR/ρRed. The NDVI has the advantage of minimizing types of band-correlated noise and influences attributed to variations in direct/diffuse irradiance, clouds and cloud shadows, sun and view angles, topography, and atmospheric attenuation. The ratio process can also help to reduce calibration and instrument-related errors. It is commonly expressed as Equation (1), which is as follows:

$$\mathrm{NDVI}=\frac{\mathrm{NIR}-\mathrm{Red}}{\mathrm{NIR}+\mathrm{Red}}$$

(1)

2.3. Measuring the Influencing Factors

We used precipitation rate (PREC), solar radiation (SRAD), air temperature (TEMP), and specific humidity (SHUM) as independent variables to measure the impact of natural environmental changes on urban vegetation [53]. The annual average wind speed and precipitation data are based on the existing Princeton reanalysis data, GLDAS (Global Land Data Assimilation System) data, GEWEX-SRB (Global Energy and Water Cycle Experiment—Surface Radiation Budget) radiation data, and TRMM (Tropical Rainfall Measuring Mission) precipitation data, which are integrated with the records of China’s meteorological stations [54]. The meteorological data with the spatial resolution of 0.1 degree was resampled to 1 km for the subsequent zonal statistics.

We used four indicators, including night light intensity (NLI), population density (POPDEN), area scale (AREA), and morphological fractal dimension (FRAC), as independent variables to measure the impact of urbanization level or urban expansion pattern on urban vegetation. The AREA index was directly calculated by using the aforementioned 30 m land use data. POPDEN was measured based on the population density dataset from the Center for International Earth Science Information Network at Columbia University (available at http://sedac.ciesin.columbia.edu/data/set/gpw-v4-population-density-rev10). We used the NLI obtained by satellite remote sensing to measure the impact of anthropogenic factors, because the traditional social survey data are difficult to match with the identified urban spatial range. The NLI mean value of each city unit extracted by remote sensing zonal statistics proved to be able to effectively approximate the city’s power consumption [55,56], GDP, intensity of social and economic activities [57,58], urbanization level [59,60,61], etc. We used NLI as a comprehensive urbanization indicator to measure the impact of comprehensive anthropogenic factors, which also helped us control climate and other factors in the regression model. The Defense Meteorological Satellite Program (DMSP)-Operational Linescan System (OLS) nighttime light (NTL) version 4 stable average visible data were obtained from the NOAA-National Geophysical Data Center (http://ngdc.noaa.gov/eog/dmsp/downloadV4composites.html), and then the NTL data were calibrated via the ridgeline sampling regression method to obtain a consistent NLI time series [62]. We used NLI of 2013 as an alternative to that of 2015 due to the lack of data. Moreover, we calculated the FRAC of each city based on the city boundary layer obtained from land use data. Compared with the basic measures of urban forms, such as length, area, and density, fractal dimension is used to differentiate shapes of cities, which tend to be circular or striped, and it is more efficient in describing space filling or sprawl of urban evolution [63,64]. In our context, a larger FRAC corresponded to a higher level of urban sprawl. According to the spatial range of the identified urban units, we made the zonal statistics for four meteorological indicators, NLI and POPDEN, and we obtained the mean value of the corresponding urban and climatic indicators of each city unit as regression variables. AREA and FRAC were calculated directly based on the vector layer of each urban patch.

3. Econometric Model

We used panel data models, which took the period of 2000–2015 into consideration, to estimate the impact of eight influencing factors on urban vegetation. In the panel model, more data points could be considered for analysis so that the degrees of freedom were increased while the collinearity was reduced [65,66]. In addition, the panel model usually had the power to control individual heterogeneity [67]. When implementing a series of panel data models, an important premise is that the natural logarithm transformation should be used for all dependent and independent variables to avoid the effects of nonstationarity and heteroscedasticity [68,69,70,71,72]. We used the following equation to create a fixed-effect model based on the panel data of the four study years:

ln(UVI_it) = β_itln(IF_it) + λ_i + u_it

(2)

where UVI_it is the urban vegetation index of urban extent i at time t; IF_it is the value of the influencing factors of urban extent i at time t; β_it is the slope coefficient; λ_i is the fixed effects in urban extent i; and u_it is the residual error of urban extent i at time t. A Hausman test was applied to examine the validity of the fixed-effect model (the result is shown in

Table A1. Result of fixed effect test.

	Coefficient
	Fixed Effect	Random Effect	Difference	S.E.
ln(PREC)	0.043	0.064	−0.021	0.003
ln(SHUM)	0.222	0.289	−0.067	0.014
ln(SRAD)	0.118	0.146	−0.027	0.016
ln(TEMP)	8.178	2.985	5.193	0.567
ln(FRAC)	−0.448	−0.092	−0.356	0.033
ln(AREA)	0.013	0.005	0.008	0.002
ln(NLI)	−0.080	−0.090	0.010	0.001
ln(POPDEN)	−0.026	−0.003	−0.024	0.015
Cons.	−37.001	−7.489	−29.512	3.152
Chi2	335.000
Prob.	0.000

in Appendix A).

4. Results

4.1. Spatiotemporal Dynamics of Urban Vegetation and Influencing Factors for All Natural Cities in China

According to the descriptive statistics (

Table 1. Statistics of status and temporal change for all cities.

		NDVI	PREC	SHUM	TEMP	SRAD	NLI	POPDEN	AREA	FRAC
	Units	1	mm/hr	10−2kg/kg	K	W/m2		pc/km2	km2
Status (2015)	Mean	0.29	0.14	0.90	287.94	158.83	29.20	1618.18	20.25	1.14
	Std	0.07	0.08	0.30	5.02	19.49	13.08	1972.62	82.59	0.04
	25%	0.25	0.07	0.67	285.62	144.07	18.19	428.78	2.78	1.12
	Median	0.29	0.12	0.94	288.93	156.17	29.09	946.93	6.87	1.14
	75%	0.33	0.21	1.10	291.04	169.63	39.72	2002.11	16.44	1.17
Difference (2000–2015)	Mean	0.00	0.03	0.03	0.62	0.50	11.20	144.36	12.03	0.01
	Std	0.05	0.05	0.07	0.69	8.97	8.47	436.31	53.85	0.03
	25%	−0.03	0.00	−0.01	0.19	−5.58	4.88	−7.47	0.91	−0.01
	Median	0.00	0.02	0.02	0.56	−0.04	10.02	45.54	3.00	0.01
	75%	0.03	0.05	0.07	1.01	5.62	16.29	186.56	9.16	0.03

NDVI: Normalized Difference Vegetation Index; PREC: precipitation; SRAD: solar radiation; TEMP: air temperature; SHUM specific humidity (SHUM); NLI: nighttime light intensity; POPDEN: population density; FRAC: fractal dimension.

), NDVI was generally stable from 2000 to 2015 (the median and the mean values of NDVI differences were close to zero), but there were large differences in the changes of many cities, with a maximum and minimum difference of 0.43. The urbanization index and the temperature showed obvious increments. The urban density and the night light intensity increased drastically, with some cities shrinking. According to the descriptive statistics of meteorological indicators over the years, we found that the overall temperature of the region had significant changes, with the 15 years annual difference up to 0.62 K.

The map of NDVI of all identified cities in 2015 is shown in Figure 2, which is basically consistent with the distribution pattern of NDVI in China’s natural regions. However, the spatial distribution of NDVI changes from 2000 to 2015 shows the urban vegetation cover in the Yangtze River Delta and the southeast coastal areas declined sharply, while the urban vegetation in the northeast and the south-central areas had a certain trend of enhancement. The spatial-temporal changes of the eight influencing factors were quite different. Among them, NLI and AREA were generally on the rise, with the most dramatic changes in the Beijing-Tianjin-Hebei agglomeration. FRAC grew in most parts of China, with the exception of some central and western regions. The increase and the decrease of population density coexisted. The population density of several developed urban agglomerations increased significantly, while the population density of some parts of central and western regions decreased. Several climate factors showed inconsistent trends in which PREC generally increased in Southwest China, while TEMP increased in most parts of the country. Cities in the north tended to be drier, while those that lie south of the Yangtze River tended to be wetter. In addition, the increase and the decrease of solar radiation intensity received from the ground were not stable, but the absolute amount of change was not substantial.

4.2. Panel Regression Results at the National Scale

The results of the Hausman test based on national scale samples showed that the fixed effect model was reasonable (Table S1). The regression results (

Table 2. Panel model estimations at the national scale.

	Coef.	Std. Err.	t	P > |t|	[95% Conf. Interval]
ln(PREC)	0.043	0.005	9.040	0.000	0.034	0.052
ln(SHUM)	0.222	0.022	9.960	0.000	0.178	0.265
ln(SRAD)	0.118	0.025	4.730	0.000	0.069	0.167
ln(TEMP)	8.178	0.642	12.740	0.000	6.920	9.437
ln(FRAC)	−0.448	0.069	−6.490	0.000	−0.583	−0.313
ln(AREA)	0.013	0.003	4.760	0.000	0.007	0.018
ln(NLI)	−0.080	0.003	−24.330	0.000	−0.087	−0.074
ln(POPDEN)	−0.026	0.015	−1.790	0.074	−0.056	0.003
cons	−37.001	3.608	−10.260	0.000	−44.073	−29.929
R2-within	0.12
R2-between	0.52
R2-overall	0.49

) showed that all eight factors had significant estimation coefficients. PREC, SHUM, SRAD, and TEMP all promoted NDVI (coefficient less than 0.3), and the coefficient of TEMP (coefficient larger than 8.0) was the largest. FRAC, NLI, and POPDEN all had a negative effect on NDVI, while AREA had a slightly positive effect. It is worth noting that the absolute value of the elasticity coefficient corresponding to FRAC in the four urbanization indicators was more than five times greater than that of the others.

4.3. The Relationships between NDVI and Influencing Factors with Different Geographical Locations

We compared the influence of various factors according to three geographical regions and six geographical regions in China (Figure 3). When the coefficient was significant, the greater the absolute value was and the greater the influence was. The results of three geographical regions showed that the impact of climate factors and urbanization factors on NDVI presented an apparent east–west gradient. From the east to the west, PREC and SHUM became more and more influential, while TEMP and SRAD had less and less of an influence. Among the urbanization related factors, NLI and FRAC became more influential in the east than in the west. POPDEN had a negative effect on the east and a positive effect centrally but no significant effect on the west. AREA mainly played a role in NDVI in the central region.

From the perspective of six more detailed geographical regions (

Table 3. The estimated coefficient of urbanization and climate factors with different geographical location: (a) with different geographical regions; (b) with different altitudes.

(a)
			Cities in Three Geographical Regions			Cities in Six Geographical Regions
			Eastern	Central	Western	Northeast	North	East	Central-South	Northwest	Southwest
ln(PREC)			0.0344 ***	0.0441 ***	0.0579 ***	0.0660 ***	0.0211	0.0397 ***	0.0275 *	0.0340 **	0.0119
			(4.63)	(4.84)	(6.73)	(4.75)	(1.90)	(4.51)	(2.49)	(2.63)	(0.67)
ln(SHUM)			0.177 ***	0.215 ***	0.225 ***	−0.336 ***	0.311 ***	0.277 ***	0.415 ***	0.224 ***	0.123
			(4.79)	(4.94)	(5.82)	(−3.95)	(5.14)	(6.40)	(8.26)	(3.40)	(1.95)
ln(SRAD)			0.140 ***	0.122 **	0.105 *	−0.196 **	0.379 ***	0.150 **	0.258 ***	0.0834	0.0560
			(3.68)	(2.61)	(2.14)	(−2.69)	(6.11)	(2.95)	(5.03)	(0.87)	(0.83)
ln(TEMP)			15.25 ***	7.113 ***	3.048 **	14.95 ***	10.93 ***	11.96 ***	−0.471	0.598	4.213 *
			(12.99)	(6.21)	(2.74)	(10.29)	(7.59)	(7.08)	(−0.30)	(0.33)	(2.19)
ln(FRAC)			−0.564 ***	−0.385 **	−0.290 *	−0.398	−0.0168	−0.842 ***	−0.415 **	−0.348	−0.533 *
			(−5.64)	(−3.25)	(−2.01)	(−1.66)	(−0.11)	(−7.50)	(−2.77)	(−1.40)	(−2.49)
ln(AREA)			0.00650	0.0355 ***	−0.00262	0.00963	0.00592	0.0155 ***	0.0253 ***	−0.00700	0.00763
			(1.76)	(7.40)	(−0.44)	(0.84)	(0.90)	(3.82)	(4.67)	(−0.63)	(0.90)
ln(NLI)			−0.0966 ***	−0.0852 ***	−0.0578 ***	−0.0243	−0.0471 ***	−0.110 ***	−0.0654 ***	−0.0422 ***	−0.0887 ***
			(−17.59)	(−14.51)	(−9.56)	(−1.57)	(−6.08)	(−19.41)	(−8.13)	(−4.12)	(−9.55)
ln(POPDEN)			−0.0841 ***	0.130 ***	0.0427	−0.0379	0.208 ***	−0.142 ***	0.0323	0.313 ***	−0.373 ***
			(−4.23)	(3.74)	(1.42)	(−0.53)	(6.24)	(−5.43)	(1.10)	(6.33)	(−7.83)
Cons.			−76.86 ***	−32.12 ***	−8.476	−76.64 ***	−55.31 ***	−57.26 ***	11.60	3.567	−12.31
			(−11.65)	(−5.00)	(−1.36)	(−9.66)	(−6.84)	(−5.98)	(1.30)	(0.36)	(−1.14)
N			4278	3279	3204	928	1473	3551	2323	1113	1394
R2 (within)			0.213	0.155	0.076	0.284	0.218	0.248	0.161	0.084	0.151
R2 (between)			0.324	0.306	0.617	0.171	0.321	0.100	0.001	0.157	0.054
R2 (overall)			0.311	0.295	0.574	0.163	0.289	0.107	0.003	0.140	0.083
rho			0.890	0.907	0.807	0.791	0.969	0.862	0.823	0.969	0.962
(b)
		Altitudes
		Low		Medium-Low		Medium		High		Ultra-High
ln(PREC)		0.0390 ***		0.0183		0.0435 ***		0.0856 ***		0.0668
		(7.04)		(1.06)		(3.52)		(3.82)		(1.03)
ln(SHUM)		0.244 ***		−0.105		0.135 *		0.148		−0.0214
		(9.31)		(−1.38)		(2.09)		(1.66)		(−0.09)
ln(SRAD)		0.128 ***		0.0653		0.250 *		−0.498 ***		−1.009
		(4.58)		(0.91)		(2.36)		(−3.96)		(−1.12)
ln(TEMP)		10.79 ***		5.547 **		−2.585		−3.369		10.12
		(14.45)		(2.96)		(−1.24)		(−1.21)		(0.90)
ln(FRAC)		−0.679 ***		0.110		0.146		0.258		−0.194
		(−8.66)		(0.49)		(0.66)		(0.82)		(−0.21)
ln(AREA)		0.0161 ***		0.0242 *		−0.00476		−0.0323 *		−0.0470
		(5.55)		(2.48)		(−0.52)		(−2.52)		(−1.34)
ln(NLI)		−0.0937 ***		−0.0285 **		−0.0531 ***		−0.0618 ***		−0.000338
		(−24.28)		(−2.76)		(−5.23)		(−4.85)		(−0.01)
ln(POPDEN)		−0.0529 **		0.247 ***		0.158 **		0.156 *		−0.644 *
		(−3.26)		(4.37)		(3.03)		(2.25)		(−2.37)
Cons.		−51.42 ***		−25.77 *		21.29		29.99		−41.18
		(−12.23)		(−2.44)		(1.83)		(1.94)		(−0.66)
N		8177		1069		967		641		61
R2 (within)		0.182		0.086		0.068		0.163		0.653
R2 (between)		0.409		0.030		0.114		0.189		0.105
R2 (overall)		0.382		0.032		0.113		0.180		0.155
rho		0.854		0.956		0.937		0.936		0.989

Note: t statistics in parentheses = “* p < 0.05, ** p < 0.01, *** p < 0.001”; rho refers to the fraction of variance.

a), PREC had the greatest influence in Northeast China, while it had the smallest influence in Southwest China. SHUM had the strongest promotion effect in Central China (coefficient 0.415) but the smallest effect in Southwest China (coefficient insignificant). SRAD had the largest promotion effect in North China but the smallest effect in Northwest and Southwest China (insignificant). In Northeast, East, and North China, the effect of TEMP on NDVI was significant, with the coefficients as 14.95, 11.96, and 10.93, respectively, while in South Central, Northwest, and Southwest China, the effect was small. The negative effect of FRAC on NDVI in East China (coefficient is −0.842) was significantly higher than that in other regions (coefficient absolute value was less than 0.53). The influence of AREA was significant in East China and Central South China, but the absolute value of coefficient was very small at 0.0155 and 0.0253, respectively. Except in Northeast China, NLI had a significant negative effect on other geographical regions, especially in East China, where the corresponding coefficient was −0.110 higher than that in other regions where most coefficients were lower than −0.06. The impact of POPDEN was the most complex, with significant positive impact (in North and Northwest China), significant negative impact (in East and Southwest China), and non-significant impact (Northeast, South Central China).

From the perspective of altitude, we divided the cities into five categories (Table 3b): low altitude cities (altitude below 500 m), medium-low altitude cities (500 m ≤ altitude < 1000 m), medium altitude cities (1000 m ≤ altitude < 1500 m), high altitude cities (1500 m ≤ altitude < 3500 m), and ultra-high altitude cities (altitude above 3500 m). The results showed that the regression coefficients of most factors were significant with low and medium altitude cities, but the regression coefficients of the factors corresponding to cities with ultra-high altitude were not significant. Only the coefficients of POPDEN were significantly positive. It is worth noting that climate factors such as SHUM and TEMP mainly played a role in cities at lower altitudes, while POPDEN and NLI had definite influence in almost all altitudes. In particular, the impact of POPDEN gradually increased with the elevation, and the regression coefficient gradually transitioned from −0.05 to −0.64, whose value increased by more than ten times. Although NLI had insignificant influence in the area of ultra-high altitude, its influence tended to be intensified from low and medium altitude to high altitude, whose regression coefficient gradually changed from −0.03 to −0.06. In contrast, FRAC only had a significant impact on the low altitude plain area with altitude below 500 m.

4.4. The Relationships between NDVI and Influencing Factors with Different Urban Characteristics

According to existing research classification standards, cities are divided into four categories according to their area sizes [69]: small-sized cities (urban area size ≤ 50 km²), medium-sized cities (50 km² < urban area size ≤ 150 km²), large-sized cities (150 km² < urban area size ≤ 250 km²), and mega cities (urban area size > 250 km²) (Figure 4a). In addition, according to the administrative region scope of urban agglomerations planned by the Chinese government, we extracted the cities of five mature urban agglomerations, including the Beijing-Tianjin-Hebei urban agglomeration (BTH), the Yangtze River Delta urban agglomeration (YRD), the Pearl River Delta urban agglomeration (PRD), the Chengdu-Chongqing urban agglomeration (CC), and the Triangle of Central China urban agglomeration (TCC) (Figure 4b).

The regression results are shown in

Table 4. The relationships between NDVI and urbanization/climate factors with different urban size and urban agglomerations.

	Area Size			Urban Agglomeration
	Small	Medium	Large	BTH	YRD	PRD	MYR	CC	Non
ln(PREC)	0.0458 ***	0.0540 *	0.0645	0.0186	0.0260	−0.0438	−0.00835	0.0161	0.126 ***
	(9.27)	(2.20)	(0.71)	(1.21)	(1.79)	(−1.07)	(−0.52)	(0.44)	(8.23)
ln(SHUM)	0.214 ***	0.221	1.003	−0.00280	0.601 ***	0.427	0.526 ***	0.186	−0.0465
	(9.30)	(1.71)	(1.01)	(−0.03)	(8.79)	(1.92)	(7.26)	(1.29)	(−0.69)
ln(SRAD)	0.113 ***	0.0774	0.578	0.392 ***	0.355 ***	0.757 ***	0.376 ***	−0.0660	0.192 *
	(4.40)	(0.63)	(0.93)	(4.73)	(4.57)	(4.64)	(3.35)	(−0.56)	(1.99)
ln(TEMP)	7.924 ***	7.665 *	−4.458	11.05 ***	27.39 ***	1.313	−5.260	11.74 **	−3.139
	(11.97)	(2.07)	(−0.26)	(5.11)	(7.90)	(0.18)	(−1.77)	(2.82)	(−1.41)
ln(FRAC)	−0.489 ***	−0.0495	0.258	−0.384	−0.738 ***	−0.836	−1.164 ***	−0.310	−0.185
	(−6.85)	(−0.10)	(0.12)	(−1.84)	(−4.43)	(−1.72)	(−5.32)	(−0.71)	(−0.92)
ln(AREA)	0.0122 ***	0.0911 ***	0.0168	0.00233	−0.00178	−0.00735	0.0359 ***	0.0263	0.00999
	(4.43)	(4.27)	(0.16)	(0.23)	(−0.30)	(−0.52)	(4.89)	(1.45)	(1.28)
ln(NLI)	−0.0780 ***	−0.344 ***	−0.471	−0.0593 ***	−0.113 ***	−0.0382	−0.127 ***	−0.119 ***	−0.0711 ***
	(−23.36)	(−7.82)	(−1.72)	(−4.33)	(−14.44)	(−1.27)	(−12.16)	(−6.54)	(−6.36)
ln(POPDEN)	−0.0338 *	−0.0325	0.149	0.124 **	−0.0417	0.104	0.123 *	−0.256 **	−0.0810
	(−2.18)	(−0.43)	(0.89)	(3.02)	(−1.24)	(1.83)	(2.26)	(−3.03)	(−1.76)
Cons.	−35.53 ***	−32.03	37.61	−56.99 ***	−144.9 ***	−1.848	38.50 *	−54.54 *	25.75 *
	(−9.54)	(−1.55)	(0.39)	(−4.70)	(−7.39)	(−0.05)	(2.31)	(−2.27)	(2.08)
N	10441	420	54	696	1491	286	1032	458	1162
R2 (within)	0.12	0.38	0.49	0.23	0.35	0.36	0.26	0.24	0.12
R2 (between)	0.51	0.37	0.30	0.18	0.04	0.11	0.13	0.11	0.00
R2 (overall)	0.49	0.39	0.26	0.17	0.09	0.03	0.06	0.10	0.00
rho	0.81	0.87	0.82	0.85	0.79	0.88	0.89	0.85	0.93

Note: * is significant at the 0.05 level (2-tailed); ** is significant at the 0.01 level (2-tailed); and *** is significant at the 0.001 level (2-tailed).

. Most of the estimated coefficients of the factors of mega cities were not significant, which may have been due to the small number of city clusters. In general, all climate factors, including PREC, SHUM, SRAD, and TEMP, showed greater impact on smaller-sized cities, while NLI had significant constant impacts on cities of all sizes but had a greater impact on larger cities (coefficient of small city was −0.08, while that of a big city was about −0.6). FRAC, POPDEN, and AREA, on the other hand, had a greater impact on smaller-sized cities.

From the perspective of urban agglomerations, urban NDVI of five mature urban agglomerations was not sensitive to PREC, while PREC played an important role in non-urban agglomeration areas with the elastic coefficient up to 0.126. While POPDEN and FRAC played an important role in urban agglomeration areas, their coefficients of non-urban agglomerations were not significant. The coefficient estimation of the factors in each urban agglomeration were distinct. Compared with other agglomerations, the coefficient value of POPDEN in Beijing-Tianjin-Hebei urban agglomeration was the largest among all of the agglomerations, while the coefficients of SHUM and TEMP in the Yangtze River Delta urban agglomeration were the largest. The Triangle of Central China urban agglomeration consisted of many scattered urban patches. Compared with other urban agglomerations, its coefficients of FRAC, NLI, and AREA were the largest. A noteworthy characteristic of the Chengdu-Chongqing urban agglomeration is that SRAD had no significant impact on its NDVI, but POPDEN had a distinct negative impact over other agglomerations whose absolute value of coefficient was more than twice that of others.

5. Discussions

5.1. Impact of Urbanization and Climate Factors from a National Perspective

The spatial pattern of NDVI changes in the past 15 years is almost the opposite of the current pattern in 2015 (Figure 2). For example, the Yangtze River Delta was the area with the best vegetation growth in 2015, but it was also the area with the most serious urban vegetation degradation. On the other hand, although the total amount of vegetation in a large number of northern cities was low, it showed a high increase in the past. This pattern difference was caused by the complex combined effect of climate factors and urbanization factors, such as the severe urbanization and urban sprawl in the Yangtze River Delta and the strong promotion of climate warming on vegetation growth in the northern region.

According to the estimation of the national scale regression model, precipitation, humidity, solar radiation, and temperature all play an important role in promoting urban vegetation because of their enhanced transpiration and photosynthesis [73]. Among the climate factors, the elasticity coefficient of TEMP was the largest, which was far higher than the other three factors, indicating that climate warming had a strong impact on urban vegetation. Among the urbanization factors, the elasticity coefficient of FRAC was negative with the absolute value far greater than that of other factors. This is because the process of sprawling urbanization often results in inefficient development and utilization of land, and, thus, more natural surface is occupied [74,75,76,77,78]. In addition, a wealth of research shows that, on the urban scale, sprawl strongly increases traffic flow, reduces traffic efficiency, increases traffic pollutant emissions, etc. [34,79,80,81,82]. Air pollution directly harms the growth of plants and indirectly affects the growth of vegetation by reducing air visibility and the intensity of solar radiation received by vegetation [83,84].

In general, the level of urban economic development represented by nighttime lighting restrains the growth of vegetation, and its mechanism is similar to that of population density, which suggests increased buildings, high-intensity development or land use, etc., may have a negative effect on the urban ecosystem, compressing the growth space of urban vegetation [85,86]. However, the expansion of urban area does not entirely mean the degradation of urban vegetation. On the whole, the expansion of urban area promotes the increase of overall urban vegetation to a certain extent, which may be related to the enhancement of green infrastructure construction in the process of urban growth [87]. In addition, urban expansion may bring or build more ecological spaces in their edging areas near suburbs. It is worth noting that our AREA is based on the area of urban categories derived from land use data rather than the impervious surface area, thus it is inconsistent with some existing conclusions based on the impervious surface area [88]. Our results provide new insights to the relationship between urban size and vegetation on an urban scale other than a pixel scale.

5.2. The Reasons and Implications of the Influence by Geographical Location

According to the subpanel models of the three geographical regions, the elastic coefficient values corresponding to climate factors and urbanization factors all showed clear gradient laws of gradual change from east to west. This gradient feature was closely related to natural and urban fundamental features underlying the geographical location [89]. For example, water resources in the eastern region were generally more abundant than those in the western region, and further westward, the more likely water was to become a limiting factor for vegetation growth. Therefore, from the east to the west, PREC and SHUM became more and more influential. However, the solar radiation in the west may have become too strong, which may have had a nonlinear effect on vegetation growth, such as the negative effect of drying caused by high temperature on vegetation growth. Therefore, from the east to the west, the elastic coefficients of TEMP and SRAD were growing smaller and smaller.

The geographical distribution of urbanization factors may be related to background social and economic characteristics of eastern and western cities. The economic development level of eastern cities is higher, and the corresponding urban building density and economic intensity in the same area are larger [52]. Therefore, in the eastern cities, the changes of biophysical and socio-economic activities brought by NLI, AREA, and FRAC with the same proportion increment were more significant, which led to the increasing trend of the elasticity coefficient corresponding to NLI, AREA, and FRAC from the west to the east.

China’s six geographical divisions are a natural environmental division as well as an economic division. Each geographical region corresponds to a unique climate, vegetation type, comprehensive development level, and local culture, which directly or indirectly affect the relationship between natural or urbanization factors and the urban ecosystem. For example, Southwest China has abundant water resources, and a small amount of additional precipitation change is difficult to have a further significant impact on urban vegetation, thus the coefficients of PREC and SHUM were small or not significant in the southwest. In contrast, the temperature in the north is generally lower than that in the south, which may make TEMP the limiting factor for vegetation growth in the north. Therefore, the increase of SRAD and TEMP in Northeast, East, and North China showed a greater role in urban vegetation enhancement. The impact of NLI on East China was much higher than other regions. This was because the NLI changes in East China were a very intense urban development process, and the occupation of natural land and the transformation of natural surfaces were also the strongest. It suggests that the investment of vegetation renovation and ecosystem improvement should be increased in developed areas such as East China. The relationship between several indices and NDVI showed a dependence on altitude (TEMP, AREA, NLI, POPDEN, FRAC), while others (PREC, SHUM, SRAD) did not. POPDEN had a certain degree of influence at all altitudes, and it had a greater influence at high altitudes (corresponding to −0.64 at ultra-high altitudes, while the absolute value of low-altitude coefficient was less than 0.3), which is related to the more vulnerable urban ecosystem in high-altitude cities. While FRAC only played a significant role in low-altitude areas, the urban sprawl in low-altitude plain areas meant the emergence of a large number of inefficient spaces on the edge of the city. The results indicate that controlling urban density at a high altitude and urban sprawl at a low altitude is of significance to urban ecosystem enhancement.

5.3. The Reasons and Implications of the Influence by Urbanization Characteristics

With the development and the expansion of the city, the fundamental characteristics of the city, including building density, industrial structure, etc., were altered, further changing the intensity and the pattern of land use, green infrastructure construction, and microclimate change and then indirectly affecting the relationship between various indicators and urban vegetation. In general, due to the larger impervious surface area and the drier microclimate [90,91], the impact of precipitation on vegetation may be greater in larger cities. However, the change of land use intensity caused by urbanization may be nonlinear, which leads to NLI playing a more important role in vegetation in larger cities. The urban sprawl of small cities and the increase of population density are also noteworthy. Due to the small size of these cities, the small proportion increase of these two indicators may have a greater negative effect on the overall urban ecosystem.

The five mega urban agglomerations we selected are the most mature, urbanized areas in China, whose development has regional effects on economic society and natural ecology. The results show that urban vegetation is very dependent on precipitation in non-urban agglomeration areas, while in urban agglomeration regions, due to more artificial irrigation and the intense construction of green infrastructure [29,92], the influencing factors of vegetation become more diverse. The key influencing factors of each city group, though, are distinct. For example, the wealth of small and medium-sized cities as well as the ultra-high population capacity in the Beijing-Tianjin-Hebei urban agglomeration is densely distributed, and, thus, the additional increase of high-density population there contributes to a far greater impact on vegetation than other urban agglomerations. Therefore, the population density control in this area is significant to the improvement of urban ecology, but it should be considered in accordance with local conditions for other urban agglomerations.

5.4. Limitations and Future Studies

There are several limitations in the research. First, the panel models were used in this paper to evaluate the impact differences of specific factors across regions. The limitation is that the statistics depend on the panel models, which lack a consideration of the nonlinear relationship between urban morphology, climatic factors, and urban vegetation. In the future, it would be interesting to study and simulate the overall impact of urban development and climate variability on an urban ecosystem with co-variates controlled by constructing refined mechanism models or machine learning models at the urban scale with full consideration of the interactions among closely related factors.

Second, we used the spatial range of the cities changing with the year to extract the average NDVI. The average NDVI of each spatiotemporal sample ensures that the samples can be compared vertically and horizontally at the same time so as to effectively carry out panel analysis. However, due to the rapid change of land use within cities, it is difficult to confirm whether the change of NDVI is caused by vegetation growth or green area change without considering more detailed land use within cities. In this paper, we could not get the data of land use within cities all over the country. The addition of these data could help us better distinguish the specific impact of urbanization on urban ecosystem for future studies.

Third, NLI is used to study the impact of the comprehensive urbanization level on an urban scale, but NLI indicators cannot more finely separate the impact of specific anthropogenic factors, such as GDP, green infrastructure investment, and so on. In the future, more specific anthropogenic indicators could be incorporated into the statistical modeling based on remote sensing inversion as well as downscaling grids of socio-economic survey data, which would provide us with deeper insight into the impact of urbanization on vegetation from multiple perspectives.

Furthermore, the study only used NDVI to measure the state of an urban ecosystem, but the net primary productivity of the vegetation and other spatial fine inversion products (such as Landsat or worldview satellite data) could be used as improved indicators. It is also necessary to further expand the observation range because, from a global scale, the impact of climate and urban development may have more feedback, which will lead to more complex spatial differentiation of impact.

6. Conclusions

In this paper, more than 3000 city samples were used to study the impact of urbanization and climate factors on urban vegetation. For the first time, the research samples covered all sizes of cities in China and supplemented the understanding of the heterogeneity of various factors from the perspective of geographical location and urbanization. The results show that the changes of urban vegetation in China from 2000 to 2015 exhibited complex spatial differences, and the degradation in the Yangtze River Delta was most significant. Panel regression estimation shows how the increase of precipitation, light, temperature, and humidity all improved urban vegetation. Nighttime light intensity, population density, and morphological sprawl of the city had serious, negative effects on the urban vegetation.

The magnitude of the effect followed obvious gradient change patterns from east to west. Among the climate factors, PREC and SHUM had more influence, while TEMP and SRAD had less influence. Among the urbanization factors, NLI, AREA, and FRAC became less influential. PREC had the greatest effect in Northeast China, SHUM had the strongest promoting effect in Central China (coefficient 0.415), SRAD had the largest promoting effect in North China, and TEMP had a great promoting effect in Northeast China, East China, and North China (14.95, 11.96, and 10.93, respectively), while FRAC and NLI played a far greater role in East China than in other regions. POPDEN had influence at all altitudes, especially at high altitudes. On the contrary, FRAC had significant influence only in low-altitude cities. NLI played a greater role in larger cities, while the urban sprawl and the population density of small cities played greater roles. Urban vegetation was very dependent on precipitation in non-urban agglomeration areas, but the influence factors of vegetation in urban areas were more diverse. Among them, the population density of Beijing-Tianjin-Hebei urban agglomeration had a far greater impact on vegetation than other areas. The heterogeneity was driven by background climate and urbanization condition.