Exploring the Dynamics of Urban Greenness Space and Their Driving Factors Using Geographically Weighted Regression: A Case Study in Wuhan Metropolis, China

Urban greenness plays a vital role in supporting the ecosystem services of a city. Exploring the dynamics of urban greenness space and their driving forces can provide valuable information for making solid urban planning policies. This study aims to investigate the dynamics of urban greenness space patterns through landscape indices and to apply geographically weighted regression (GWR) to map the spatially varied impact on the indices from economic and environmental factors. Two typical landscape indices, i.e., percentage of landscape (PLAND) and aggregation index (AI), which measure the abundance and fragmentation of urban greenness coverage, respectively, were taken to map the changes in urban greenness. As a case study, the metropolis of Wuhan, China was selected, where time-series of urban greenness space were extracted at an annual step from the Landsat collections from Google Earth Engine during 2000–2018. The study shows that the urban greenness space not only decreased significantly, but also tended to be more fragmented over the years. Road network density, normalized difference built-up index (NDBI), terrain elevation and slope, and precipitation were found to significantly correlate to the landscape indices. GWR modeling successfully captures the spatially varied impact from the considered factors and the results from GWR modeling provide a critical reference for making location-specific urban planning.


Introduction
More than half of the world's population lives in cities [1]. The World Urbanization Prospects 2018 estimated that 68% of the world's population would dwell in cities by 2050 [2]. In the past few decades, urban land across the globe has expanded dramatically from urbanization [3], which is also often accompanied by spatially fragmented and socially uneven urbanized areas [4]. The process of urbanization deeply affected vegetation patterns and ecosystem services, as it could significantly update the material circulation and energy flow of an area [5,6]. As a typical outcome in terms of land fragmentation from urbanization, a large vegetation patch is often cut into smaller pieces, which could easily result in diminished urban greenness. It is well recognized that urban greenness space provides various ecosystem services to inhabitants in cities [7,8]. Compared to other land cover types, urban greenness coverage has a more profound impact on the sustainable development within and around a city [9]. For example, Dennis et al. confirmed that urban green infrastructure had a key Land 2020, 9, 500 3 of 21 namely, percentage of landscape (PLAND) and aggregation index (AI), were computed to measure the urban greenness patterns at regional and local scales. Because PLAND reflects the magnitude or abundance of urban greenness space [35], while AI indicates the structure or fragmentation in vegetation [37], the adoption of the two indices could provide a better understanding of the dynamics of urban greenness space. There have been a variety of approaches to explore the factors for urban greenness space. Traditional multivariate linear regression (MLR) provides a simple method to study the relationships between multiple variables. The advantage of MLR is that the result is clear and easy to interpret. However, in most cases, spatial phenomenon demonstrates spatial heterogeneity or non-stationarity, which could not be revealed by an MLR approach. More advanced approaches, such as machining learning or an artificial neural network (ANN), are useful in understanding complicated patterns in geospatial fields but they are essentially a black-box, which makes the result difficult to interpret in terms of which predictors are more important than others or how they are related to the target variable [38]. We hypothesize that both natural and human factors could impose complicated and spatially varied impact on urban greenness space. Therefore, the factors that affect the dynamics of urban greenness space were analyzed using geographically weighted regression (GWR) with input from economic and environmental variables. The benefit of the GWR model is that the outcome not only presents clear quantitative relationships between the predictors and the target property, but also highlights the spatial variation of the relationships. We expect that the spatially varied impact of the factors from the GWR modeling reveals the interactions between the urban greenness landscapes and the economic and environmental factors and that the output from the study may provide a basis for making decisions on sustainable urban development.

Study Area
The study was performed for Wuhan metropolis, Hubei province, covering nearly 10,000 km 2 in central China, as shown in Figure 1. Wuhan, as the capital city of Hubei province, is a megacity in central China. It is one of the 17 cities included in the global sustainable development plan issued by the United Nations Development Programme and the United Nations Environment Programme [39]. Located at the intersection of the Yangtze River and Han River, the metropolis is spatially divided into three districts by the rivers, including Hankou, Wuchang, and Hanyang. Famous for its hundreds of water bodies (lakes) and rivers, it is well-known as a "river city". The city is also characterized by uneven topography (ranging between 13 and 262 m above sea level), complicated water network, dense population (over 10 million), and intensified human-related activities. The area is covered by various land cover and land use types. Urbanized or built-up areas dominate mainly around the central part of the region. Croplands, forests, and grasslands spread mainly over the peri-urban or rural areas. Continuous urban expansion has transformed vegetation into built-up areas, causing significant shrinkage in the urban greenness space and imposing negative impacts on the urban ecosystem [40]. Therefore, understanding the changes of urban greenness space and the driving factors of the changes have been urgently needed to promote a sustainable urban ecosystem for Wuhan metropolis.

Data Sources
The data sources used in the study can be grouped into 4 categories-human-related socio-economic activities, climate, topography, and urban greenness indicators-with the first three serving as explanatory variables and the last as target or dependent variables, as shown in Table 1. The details of the datasets are described below.
Google Earth Engine (GEE) is a planetary-scale platform for earth science data and analysis, useful for researchers to detect changes of the earth's surface [41]. The latest updated transportation dataset of road network and railways in 2018 was obtained from the transportation bureau of the local government. Furthermore, rivers and waterways were classified from Landsat images on the  Figure 1. The study area, Wuhan metropolis; (a) the location of the province, (b) boundary of the province, (c) Wuhan metropolis area.

Data Sources
The data sources used in the study can be grouped into 4 categories-human-related socioeconomic activities, climate, topography, and urban greenness indicators-with the first three serving as explanatory variables and the last as target or dependent variables, as shown in Table 1. The details of the datasets are described below. Google Earth Engine (GEE) is a planetary-scale platform for earth science data and analysis, useful for researchers to detect changes of the earth's surface [41]. The latest updated transportation dataset of road network and railways in 2018 was obtained from the transportation bureau of the local government. Furthermore, rivers and waterways were classified from Landsat images on the GEE platform. Those variables are closely related to the local economy and thus are related to urban land changes.
Climate conditions can alter urban greenness landscapes. Precipitation and temperature are two key variables, as they may affect the urban greenness based on a previous study [42]. Precipitation and temperature at a monthly scale composed of more than 650 climatic stations were obtained from the China Meteorological Data Sharing Platform during the years 2000-2018 (http://cdc.cma.gov.cn) and spatially interpolated using thin-plate smoothing splines [43]. Annually accumulated precipitation and mean temperature were then clipped by the study boundary. A digital elevation model (DEM) was used to represent elevation. DEM data, which are Global Multi-resolution Terrain Elevation Data 2010 with a spatial resolution of about 30 m, were acquired  Climate conditions can alter urban greenness landscapes. Precipitation and temperature are two key variables, as they may affect the urban greenness based on a previous study [42]. Precipitation and temperature at a monthly scale composed of more than 650 climatic stations were obtained from the China Meteorological Data Sharing Platform during the years 2000-2018 (http://cdc.cma.gov.cn) and spatially interpolated using thin-plate smoothing splines [43]. Annually accumulated precipitation and mean temperature were then clipped by the study boundary. A digital elevation model (DEM) was used to represent elevation. DEM data, which are Global Multi-resolution Terrain Elevation Data 2010 with a spatial resolution of about 30 m, were acquired from GEE. The topographic slope was derived from DEM on GEE [44]. The terrain characteristics may affect vegetation cover and thus should be examined in terms of their impact on urban vegetation cover.
The normalized difference vegetation index (NDVI) is often used as a proxy for vegetation coverage [45]. It is computed through the two reflectance bands, near infrared (NIR) and red (R), from remote sensing imagery in the form of GEE provides online access to archived Landsat data, which includes Landsat 5 Thematic Mapper (TM) from 1984 to 2012, Landsat 7 Enhanced Thematic Mapper Plus (ETM+) from 1999 to present, Land 2020, 9, 500 5 of 21 and Landsat 8 Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS) from 2013 to present [41,46,47]. ETM+ suffered from the failure of the scan line corrector (SLC) starting from May 2003, which resulted in data gaps in the acquired images. Image preprocessing included gap filling for ETM+ SLC-off images and cloud filtering. A simple and effective method for filling gaps in ETM+ SLC-off images using the neighborhood similar pixel interpolator (NSPI), which works especially well in heterogeneous regions as detailed by Chen et al. [48] and Asare et al. [49], was applied to improve ETM+ image quality. To minimize the effects of clouds and cloud shadows from the images, a cloud score algorithm based on a combination of brightness, temperature, and the normalized difference snow index was applied to compute a cloud-likelihood score based on GEE's built-in algorithm (i.e., ee.Algorithms.Landsat.simpleCloudScore) and pixels with the score greater than 10 were excluded [50,51]. The preprocessed Landsat imagery series enabled the computation of long-term NDVI series for the study region. An annual maximum composite NDVI for each year during 2000-2018 was derived from the computed NDVI time-series.
The normalized difference built-up index (NDBI) is an index representing the cover of the impervious surface of a city and is an indicator of human activities in urban areas [52]. NDBI is strongly correlated to the impervious surface of the city [53]. To compute NDBI, the reflectance of the two bands of the middle infrared (MIR) and near infrared band are taken, in the form of [53,54] A composite of maximum NDBI was derived for each year during 2000-2018 based on the multiple time-series maps of NDBI with that year.

Methods
The overall framework for the analysis of the dynamics and the driving factors of the urban greenness space is shown in Figure 2.
First, urban greenness patches were extracted based on image classification using GEE [40]. An urban greenness patch was defined as a continuous vegetated area covered by either grassland, forestry, croplands, or any combination of the above. Two urban greenness landscape indices, namely percentage of landscape (PLAND hereafter) and aggregation index (AI hereafter) were computed to build a time-series dataset. PLAND measures the area of the land cover type standardized by the total area of the landscape [34], while AI reflects the degree of aggregation or extension of patches for that type [55]. PLAND is an indicator used to measure the proportion of urban greenness landscapes compared to a region of interest. The abundance of greenness space is indicated by PLAND. The temporal changes, either an increase or decrease, in the area of greenness space can be illustrated by inter-annual dynamics of PLAND. For a given region R, PLAND is computed as where A R is the total area of R, a i is the area of the greenness patch i (i = 1, 2, . . . , n, where n is the total number of urban greenness patches). The aggregation index (AI) indicates the clustering pattern of urban greenness cover, which is a useful measurement for vegetation fragmentation. For the urban greenness landscape, AI is defined as the ratio of the edges shared by greenness pixels to the total number of edges, in the form of AI green = g green−green /max green (4) where g green-green is the number of all shared edges between pixels of greenness patches, max green is the maximum number of edges between pixels of greenness patches. The highest level of clustering indicates greenness pixels sharing the most possible edges (AI will be close to 1), while the lowest Land 2020, 9, 500 6 of 21 level of aggregation means greenness pixels sharing no edges (AI = 0, greenness pixels completely disaggregated) [37]. To examine the temporal dynamics of the urban greenness space, Sen's slope and the Mann-Kendall trend test were performed on the time-series of the NDVI and the landscape indices. To account for the possible driving forces, Pearson correlation analysis was applied to quantify the correlation between the landscape indices and the explanatory variables. Then two regression models, namely, ordinary least square (OLS) and the geographically weighted regression (GWR) model, were performed to reveal the global and locally varied relationships between the explanatory variables and the target landscape indices. The coefficients between the explanatory variables and the indices of urban greenness landscape from GWR were mapped so that the location-dependent impact from the driving forces could visualized.  First, urban greenness patches were extracted based on image classification using GEE [40]. An urban greenness patch was defined as a continuous vegetated area covered by either grassland, forestry, croplands, or any combination of the above. Two urban greenness landscape indices, namely percentage of landscape (PLAND hereafter) and aggregation index (AI hereafter) were computed to build a time-series dataset. PLAND measures the area of the land cover type standardized by the total area of the landscape [34], while AI reflects the degree of aggregation or extension of patches for that type [55]. PLAND is an indicator used to measure the proportion of urban greenness landscapes compared to a region of interest. The abundance of greenness space is indicated by PLAND. The temporal changes, either an increase or decrease, in the area of greenness space can be illustrated by inter-annual dynamics of PLAND. For a given region R, PLAND is computed as where AR is the total area of R, ai is the area of the greenness patch i (i = 1, 2, ..., n, where n is the total number of urban greenness patches). The aggregation index (AI) indicates the clustering pattern of urban greenness cover, which is a useful measurement for vegetation fragmentation. For the urban greenness landscape, AI is defined as the ratio of the edges shared by greenness pixels to the total number of edges, in the form of PLAND and AI can be processed globally and locally. When the entire region is selected, a global version of PLAND or AI will be defined. Conversely, a local version of the indices is processed for a predefined neighborhood. In the current study, a circular moving window with a radius of 200 pixels was defined to compute a local PLAND and AI, which is a local version capable of highlighting the spatially varied features of the indices.
The urban vegetation (through NDVI) and the urban greenness landscape (PLAND and AI) showed changes over the years. To explore the dynamics along the years, a linear regression analysis and Sen's slope were applied [56]. For a set of pairs (i, x i ) or (j, x j ), where x i (or x j ) is a time-series and i (or j) is the sequential position of the element in the dataset, Sen's slope β is calculated as A positive β suggests that the time-series show a rising trend, and vice versa. The significance of trend in the greenness landscape indices (PLAND and AI) and vegetation greenness level (NDVI) was evaluated using the Mann-Kendall test. The Kendall τ (−1 to +1) obtained from the Mann-Kendall Land 2020, 9, 500 7 of 21 statistics could be used to measure the significance of the trend [57]. Thus, the dynamics of NDVI, PLAND, and AI during 2000-2018 were tested using Sen's slope and Mann-Kendall trend analysis.
The spatial distribution patterns of the urban greenness landscape could be attributed to various factors. To examine the factors that may have impacted on the landscape patterns, three steps were taken to investigate the relationship between the temporally averaged results of the explanatory and target variables. First, Pearson correlation analysis, which takes two variables X and Y as input, computes correlation coefficient (P) or the ratio of the covariance between the variables (Cov(X,Y)) to their standard deviation products [58], using the formula The computed coefficient ranges between −1 to 1 and a value close to 0 indicates the least correlation between the variables [59]. By computing the correlation coefficients between the input factors and the target landscape indices, those explanatory factors that show significantly correlated to PLAND and AI could be retained for further analysis.
Ordinary least squares (OLS) regression is a global statistical method that estimates the relationship between independent variables and a dependent variable. With the input of the explanatory variables and the landscape indices, OLS can output the coefficients that indicate the direction and magnitude of the impact from the independent variables. Nevertheless, due to non-stationary or heterogeneous characteristics of most spatial patterns, a widely applied local regression model, namely, geographically weighted regression (GWR), has been applied in many applications to capture spatially varied relationships between multiple variables [60]. Compared to the global regression model such as OLS, GWR models the local relationships between predictors and a dependent variable, in the form of where y i is the dependent variable at location i, x ik (k = 1, 2, . . . , n, and n is the number of explanatory variables) is the independent variable, (u i , v i ) is the coordinate of the location i, β 0 (u i , v i ) is the intercept of a constant for spatial variability, β k (u i , v i ) is the estimated coefficient for the kth explanatory variable, and ε i is the random error. The kernel and bandwidth are key parameters in the GWR model that affect model accuracy. Gaussian kernel was decided in the current work. The optimal bandwidth was determined on estimated distance from spatial autocorrelation of the dependent variables (PLAND and AI). Explanatory variables input into OLS and the local model need to be free of multicollinearity [61].
To avoid multicollinearity among the variables, variance inflation factor (VIF) can be checked, where VIF values greater than 10 indicate multicollinearity [62]. The adjusted coefficient of determination (R Adj.

)
and corrected Akaike information criterion (AICc) serve as a measure of goodness of fit to determine whether a global or local model should be applied to have a better interpretation for the patterns in the analyzed data [60].

Temporal Dynamics of Urban Greenness
A time-series analysis using bivariate linear regression and Sen's slope against the years, as well as the Mann-Kendall trend test on the urban greenness landscape indices, i.e., PLAND and AI, revealed that not only the total area of urban greenness space showed a decrease, but also the urban greenness patches tended to be more fragmented, as both PLAND and AI presented significant decline over the years. Table 2 shows the results from the trend analysis by bivariate linear regression and Sen's slope. Both the linear slope and Sen's slope indicated a significantly negative value for PLAND and AI (p < 0.01, n = 19). The result revealed that more vegetated area was converted to non-vegetated land cover, most likely due to rapid urban expansion. At the same time, larger vegetated patches were found to be split into smaller ones, most likely due to intensified human activities. Nevertheless, the result demonstrated that the averaged NDVI still showed an overall increasing trend over the years, meaning the health of urban vegetation was improving at least for some areas that were not converted into built-up areas. Figure 3 shows two typical NDVI maps for the starting year (2000, Figure 3a) and the end (the year 2018, Figure 3b). The urbanized area expanded outside significantly from the urban center, explaining the reduction of the urban greenness space in the region. However, the overall greenness in terms of NDVI showed improvement in the year 2018 if compared to that in 2000. Though the detailed explanation remains undiscovered, the improvement of NDVI over vegetated areas might be attributed to recent reforestation programs in the city [63] and/or global warming [64].

Factor Analysis for the Dynamics of the Urban Greenness
The preliminary analysis selected a set of candidate factors, including the climatic variables (precipitation and surface temperature), terrain (DEM and slope), transportation impact (density of railways, roads, and waterways), and urbanization level (in the proxy of the normalized difference built-up index or NDBI). Those factors have proved to be related to the changes of the urban greenness landscape by many other studies [3,20]. The combination of those economic and environmental factors may provide a more comprehensive view for understanding the driving forces for the urban greenness space. First, Pearson correlation analysis was conducted to examine if there exists a strong bivariate correlation between any of the candidate variables and the urban greenness landscape indices. The correlation coefficients between the candidate variables and the indices (PLAND and AI) are shown in Table 3. Except for the temperature, all the other variables presented a significant correlation to the two urban greenness landscape indices. The significant and positive coefficients from the variables DEM, slope, and precipitation suggest strong positive effect to the improvement of the urban greenness landscape indices, while the negative values from the NDBI and transportation pathways (waterways, roads, and railways) indicated that they reduced the landscape indices.

Factor Analysis for the Dynamics of the Urban Greenness
The preliminary analysis selected a set of candidate factors, including the climatic variables (precipitation and surface temperature), terrain (DEM and slope), transportation impact (density of railways, roads, and waterways), and urbanization level (in the proxy of the normalized difference built-up index or NDBI). Those factors have proved to be related to the changes of the urban greenness landscape by many other studies [3,20]. The combination of those economic and environmental factors may provide a more comprehensive view for understanding the driving forces for the urban greenness space. First, Pearson correlation analysis was conducted to examine if there exists a strong bivariate correlation between any of the candidate variables and the urban greenness landscape indices. The correlation coefficients between the candidate variables and the indices (PLAND and AI) are shown in Table 3. Except for the temperature, all the other variables presented a significant correlation to the Land 2020, 9, 500 9 of 21 two urban greenness landscape indices. The significant and positive coefficients from the variables DEM, slope, and precipitation suggest strong positive effect to the improvement of the urban greenness landscape indices, while the negative values from the NDBI and transportation pathways (waterways, roads, and railways) indicated that they reduced the landscape indices. While the coefficients from the correlation analysis revealed the relationship between each individual variable and the indicators of the urban greenness landscape (i.e., PLAND and AI), multivariate analysis with input from the combination of the factors could be a better overview, as the selected individual variables may have multicollinearity in terms of the explanation effect to the dependent variables. Thus, the multicollinearity for the selected factors that are significantly correlated to the greenness landscape indices was further tested. Results indicated that the waterways and railway presented a serious multicollinearity (VIF > 8), as shown in Table 4. Transportation accessibility is usually regarded as a critical factor that may alter vegetation cover [23]. Out of the three selected transportation elements, the distance to waterways and the density of railways presented strong correlation with the road density. The multicollinearity suggests that the variable road alone could provide sufficient explanation power to the patterns of urban greenness space and that the waterways and railway could be replaced by the road network in the multivariate analysis. Therefore, the variables of waterway and railway were eliminated from further analysis.

Analysis of Driving Factors on Urban Greenness Landscape Based on GWR
The variables, namely, DEM, slope, NDVI, road, and precipitation, were selected based on the result from the correlation and multicollinearity analysis. Ordinary least squares (OLS) regression and GWR analysis were performed on the landscape indices PLAND and AI, respectively. Figure 4 shows that a moderate model performance was obtained from the OLS regression (adjusted R 2 = 0.604 and 0.593 for PLAND and AI, respectively), indicating a good fit of the OLS model [65], which further confirms the explanation power of the selected variables to the dependent variables. While OLS creates global estimates for the entire region, the GWR generates a set of spatially varying estimates defined by a local model. Thus, GWR differentiates location-dependent impact on the target indices from the variables. The result indicates that the GWR model produced a higher adjusted R 2 for both PLAND and AI indices, suggesting GWR is more reliable than OLS. Furthermore, the AICc from the GWR model was smaller than that from the OLS model for both PLAND and AI. Thus, applying the GWR model is more preferable in terms of the explanation power of the impact on the urban greenness landscape indices from the selected variables. and 0.593 for PLAND and AI, respectively), indicating a good fit of the OLS model [65], which further confirms the explanation power of the selected variables to the dependent variables. While OLS creates global estimates for the entire region, the GWR generates a set of spatially varying estimates defined by a local model. Thus, GWR differentiates location-dependent impact on the target indices from the variables. The result indicates that the GWR model produced a higher adjusted R 2 for both PLAND and AI indices, suggesting GWR is more reliable than OLS. Furthermore, the AICc from the GWR model was smaller than that from the OLS model for both PLAND and AI. Thus, applying the GWR model is more preferable in terms of the explanation power of the impact on the urban greenness landscape indices from the selected variables.  The elevation variable DEM presented either positive or negative impact on the landscape indices. The elevation and the regression coefficients with PLAND and AI are presented in Figure 5. The high elevated areas in the northern part did not show much impact on the two indices. The central part and southwestern areas imposed positive impact on both PLAND and AI, as revealed by the high positive coefficients, suggesting that areas at lower elevation tended to maintain good urban vegetation coverage (higher PLAND) and keep vegetation clustered (higher AI). However, the The elevation variable DEM presented either positive or negative impact on the landscape indices. The elevation and the regression coefficients with PLAND and AI are presented in Figure 5. The high elevated areas in the northern part did not show much impact on the two indices. The central part and southwestern areas imposed positive impact on both PLAND and AI, as revealed by the high positive coefficients, suggesting that areas at lower elevation tended to maintain good urban vegetation coverage (higher PLAND) and keep vegetation clustered (higher AI). However, the negative impact from elevation is more complicated. For example, the relatively high elevated area in the southern part showed a negative coefficient. It is noteworthy that the impact of elevation on PLAND and AI showed an identical pattern in terms of the spatial distribution of each coefficient, which suggests that elevation may impose similar impacts on both the abundance and fragmentation of urban greenness space in the current study area.
The slope variable, as well as its impact on the two landscape indices are presented in Figure 6. As the slope was derived from the elevation (DEM), the northern part, which was dominated by mountains, showed steep slopes. The coefficients of its impact on the landscape indices presented similar patterns to those from the DEM. A few clustered places in the central and southwestern parts of the region showed positive impact to PLAND and AI. The areas with steep slope (e.g., the northern part) did not present an obvious impact on the two indices. Those locations that showed significant negative coefficients should be given special attention as those areas tended to suffer substantial loss in urban greenness space or to fragment.   The urbanization process partially reflects the impact of human activities on urban ecosystems. The study area has experienced significant urban expansion, or increased built-up area, during the years. The urbanization of the city expanded from the center of the area and higher NDBI values were mainly distributed around the central part, as indicated in Figure 7a. The regression coefficients of the NDBI with PLAND and AI, as shown in Figure 7b,c, illustrated that the values were negative around the region of the central part and that urbanization in the greater part of the central area was likely to reduce urban greenness coverage and induce vegetation cover fragmentation. The areas that showed positive correlations between NDBI and the urban greenness indices mainly located close to water bodies. In fact, while considerable water bodies were found to be converted into urbanized land in Land 2020, 9, 500 13 of 21 the past two decades, it is not uncommon to see mosaicked vegetation in the newly built-up areas in the study area [66]. Thus, urban greenness space showed improvement around neighborhoods of water bodies.
Precipitation proves to be a significant factor to vegetation cover. The annually accumulated precipitation ranged between 1100 and 1300 mm, with an increasing trend from northeast to southwest, as shown in Figure 8a. Many studies have revealed a positive contribution of precipitation to vegetation growth in urban ecosystems [67]. We found that the impact of the precipitation on the urban greenness landscape was heterogeneously distributed with significant spatial variations. The negative and positive impact from precipitation on the landscape indices was found to present a clustering pattern, as shown in Figure 8b,c, suggesting the complicated impact on the urban greenness distribution from the precipitation variable. The variation in the impact from precipitation was likely related to the hundreds of water bodies (rivers and lakes) in the study area, making the effect of precipitation on PLAND and AI spatially varied. The distribution of water bodies may affect the environment of the nearby areas and thus the impact from precipitation did not present a uniform pattern over the whole study region.  The road network serves as another proxy for human activities which imposed spatially varied impact on urban greenness. Road network density was highest around the central region, as shown in Figure 9a, where the road density showed either negative or no obvious impact on the greenness indices, as shown in Figure 9b,c, meaning a highly densified transportation network tends to decrease urban vegetation space and fragment urban greenness patches. In contrast, the road network showed a strong positive effect on the landscape indices in some parts of the peri-urban areas, which means that the greenness space there could be improved with the building of a limited density transportation network. It is of policy implication to note that the road network might impose the opposite impact on the urban vegetation space, depending on the current road network density.

Discussions
Urban vegetation plays vital role for the sustainable development of urban ecosystems. With rapid urbanization and increasing intensified human-related activities in many cities, decreased ecological services have been observed due to significant loss of urban greenness [13,68]. Many studies have confirmed that urbanization imposed a direct impact on urban greenness space [10], and reduced exposure to greenness could cause environmental stresses to public health issues [1]. Because the world population continues to rise while cities provide more advantageous conditions, a significant amount of people in rural areas tend to live in urban areas to enjoy a convenient life [38]; thus, the increased population plays as a strong potential driver for urbanization, which may occupy urban greenness space.
The urban greenness landscape indices, i.e., PLAND and AI, prove to be useful in revealing the urban vegetation coverage and fragmentation/clustering patterns in distribution. Through spatial analysis, it is meaningful to know the changing trend in the landscape indices, as urban greenness space is closely related to public health of urban residents [7,58]. Moreover, investigating the underlying factors that are contributing to the changes would be even more important for urban land planning. Previous studies indicated that the influential factors of urban greenness were multiple and that misleading results might be drawn without a systematic analysis on its driving forces [17,42]. Those factors are composed of not only natural environments but also human-related activities [69]. Those influencing factors could impose a complex impact on the target urban greenness landscape [70].
The urban greenness changes and the possible driving factors for the changes can be explored through ordinary least square (OLS) regression and geographically weighted regression (GWR). Remotely sensed imagery, namely, the Landsat series, provides a valuable dataset to map the dynamics of the urban greenness space. There are a number of indices that could be applied to quantify urban greenness space. This study measures urban greenness patterns by two landscape indices, i.e., PLAND and AI. Since PLAND captures the proportion of urban greenness cover out of the total land area and AI measures the aggregation of greenness patches, the combination of the two could reflect the overall patterns of the urban greenness space. The time-series analysis of Sen's slope and Mann-Kendall trend testing showed usefulness in revealing the dynamics of the indices over the years. The results showed that the coverage of the urban greenness space had significantly decreased, which is in accordance with our previous study [19]. Moreover, the study also found that the urban greenness was increasingly fragmented over the years.
The identified driving forces of the landscape indices, as revealed from the OLS and GWR analysis, suggest that both the natural/environmental factors (including elevation, topographic slope, and precipitation) and human-related activities (e.g., urbanization intensity via NDBI and the road network) showed significant impact on urban greenness space in Wuhan city. Overall, the landscape indices were positively correlated with the natural factors but negatively correlated with the human-related activities. Because the natural and social factors demonstrated high variation across the region, GWR was applied to map the spatially varied impact, revealing the influential directions (either positive or negative) and magnitude from each of the factors. Those spatially varied patterns were illustrated in the coefficients of the factors against the target PLAND and AI. Urban planning policies targeting better urban ecosystem services should give more attention to urbanization control and location-dependent construction of the transportation networks and sensitive areas with urban greenness space that could be significantly affected by terrain characteristics and climatic factors. The identified driving factors and their spatially varied impact on the urban greenness indices can provide useful information to design policies targeting sustainable development of urban ecosystems.

Conclusions
This study has selected Wuhan, China as the case area to examine the temporal changes of urban greenness space and the driving factors explaining the patterns in the changes. The combination of the landscape indices, PLAND and AI, is useful to reveal both the abundance and clustering pattern of urban greenness space. It was found that the coverage of the urban greenness space had significantly decreased, and the urban greenness patches were more fragmented. Human and natural factors imposed considerable impact on urban greenness space. Road network density, normalized difference built-up index (NDBI), terrain elevation and slope, and precipitation were significantly correlated to the landscape indices. The study verified that GWR is more reliable than OLS in studying the relationship between the spatial variables. Applying the GWR model is more preferable in terms of the explanation power of the impact on the urban greenness landscape indices from the selected variables. GWR modeling also captures the spatially varied impact from the factors, which is more desirable for making location-specific urban planning.
The current work did not differentiate the vegetation cover types, for example, farmland crops or other plant types such as trees, shrubs, and grasses. The changes in those vegetation cover types could be induced from human activities (cultivation) or natural process (succession) and thus involves different mechanisms for the changes in greenness space. Future work may pay attention to the impact on different vegetation types (forestry, grassland, and croplands) from the driving factors.

Conflicts of Interest:
The authors declare no conflict of interest.