Assessing Impacts of Urban Form on Landscape Structure of Urban Green Spaces in China Using Landsat Images Based on Google Earth Engine

The structure of urban green spaces (UGS) plays an important role in determining the ecosystem services that they support. Knowledge of factors shaping landscape structure of UGS is imperative for planning and management of UGS. In this study, we assessed the influence of urban form on the structure of UGS in 262 cities in China based on remote sensing data. We produced land cover maps for 262 cities in 2015 using 6673 scenes of Landsat ETM+/OLI images based on the Google Earth Engine platform. We analyzed the impact of urban form on landscape structure of UGS in these cities using boosted regression tree analysis with the landscape and urban form metrics derived from the land cover maps as response and prediction variables, respectively. The results showed that the three urban form metrics—perimeter area ratio, road density, and compound terrain complexity index—were all significantly correlated with selected landscape metrics of UGS. Cities with high road density had less UGS area and the UGS in those cities was more fragmented. Cities with complex built-up boundaries tended to have more fragmented UGS. Cities with high terrain complexity had more UGS but the UGS were more fragmented. Our results for the first time revealed the importance of urban form on shaping landscape structure of UGS in 262 cities at a national scale.


Introduction
Urban green spaces (UGS) contribute significantly to quality of life [1,2] and sustainable development in cities [3] because they can deliver multiple ecosystem services to urban inhabitants [4,5].UGS can regulate climate [6][7][8], remove air pollution [9], reduce stormwater run-off [10], provide recreational opportunities [5], and improve public health [11].For example, the urban forest in Beijing could remove 772 tons of particulate matters annually [9].In UK, urban inhabitants were found to have lower mental distresses and higher wellbeing in places with more UGS coverage [11].Extensive studies have shown that landscape structure of UGS is closely related to ecosystem services that they can deliver [12][13][14][15].For example, clustered UGS were found to enhance local cooling while dispersed UGS can achieve larger overall regional cooling [13].Carbon densities were found to be strongly and positively related to the total tree cover and the proportion of tree cover in the surrounding areas [14].Patch sizes and corridors of UGS were found to have strong positive effects on biodiversity [16].Loss and fragmentation of UGS were found to cause declines in bird species richness [15].Therefore, understanding the pattern of landscape structure of UGS and its shaping forces is important for urban ecosystem service management.
Remote sensing techniques have been widely used to investigate landscape structure of UGS.For example, at the single city scale, changes of landscape structure of UGS have been quantified for Jinan, China [17], Shenzhen, China [18], Hong Kong, China [19], and Hanoi, Vietnam [20] using medium (e.g., Landsat) or high resolution (e.g., Quickbird) satellite images.Spatiotemporal patterns of UGS of nine Chinese cities were investigated using ALOS, SPOT, and Landsat images [21].The results showed that UGS in these cities were highly fragmented and dynamic.At the regional scale, landscape structure of UGS of 111 Southeast Asian cities was investigated using Landsat images and considerable variation of UGS among these cities was identified [22].Landscape structure of urban vegetation of 100 cities around the world was quantified using Landsat images and varied percentage of urban vegetation, and different degrees of fragmentation were found for these cities [23].So far, the majority of studies focused at single city scale and studies of large spatial scales are still limited, especially for countries and regions that are experiencing rapid urbanization such as China [24].This might be due to the fact that obtaining landscape structure of UGS at large spatial scales using the traditional remote sensing data processing tools is labor-intensive and time-consuming.Recently, the emergence of high-performance cloud platforms like Google Earth Engine (GEE) [25] makes it possible to process massive remote sensing data more efficiently [26][27][28], and provides opportunities to study the general pattern and driving mechanism of landscape structure of UGS at large spatial scales.
Except for limited studies on landscape structure of UGS at large spatial scales, our knowledge on shaping factors of landscape structure of UGS are also limited.Existing studies have shown that cities with higher population densities tend to have smaller [22,23] and more fragmented UGS [23].In Southeast Asia, rich cities have more UGS than poor cities [22].Within cities, unequal economic development leads to more fragmented urban vegetation [23].Megacities with lower annual mean temperature and higher annual precipitation were found to have more UGS relative to other large cities [28].Apart from these climatic and socioeconomic factors, there are indications that urban form (i.e., spatial configuration of fixed elements within a city [29]) may have important influences on landscape structure of UGS.For example, a negative relationship between the extent of UGS and total length of the road network has been identified in Sheffield, UK [30].Terrain roughness was found to be positively correlated with levels of tree canopy in Cincinnati, USA [31].However, very few studies addressed the impact of urban form on landscape structure of UGS, especially at a large spatial scale.How urban form can affect landscape structure of UGS is still an open question.
In this study, we assessed the impact of urban form on landscape structure of UGS through a case study in China.We selected 262 Chinese cities (population >500,000) as our study sites.With information on landscape structure of UGS and urban form extracted out from Landsat images using the GEE platform, we intend to answer the following questions.What are features of landscape structure of UGS in these cities? How does the urban form of these cities affect landscape structure of UGS in these cities? Specifically, how does landscape structure of UGS vary with topography, road density, and boundary complexity?

Study Area
We selected 262 cities in China that had an urban population of more than 500,000 by 2014 (Figure 1, Table S1, Supplementary file 1).In total, these 262 cities accounted for 59% of China's urban population [32].Most of these cities are located in east China while a few cities are located in west China.

Landsat Image Processing on GEE Platform
We used 3921 scenes of standard terrain corrected (L1T) Landsat 7 ETM+ and 8 OLI/TIRS imagery with a cloud cover <70% acquired in 2015.Winter images were not used as we mainly focused on the growing season.We also included 2752 images acquired in 2014 and 2016 for 63 cities where cloud contamination was severe.The Landsat data format is the atmospherically corrected surface reflectance Landsat Collection 1 Level-1 data, which is produced by the United States Geological Survey (USGS/EROS).All images were hosted on GEE.We used the pixel_qa band of the Landsat images to mask the cloud, cloud shadow, and snow pixels [33][34][35].
We selected the blue, green, red, near infrared, and two shortwave infrared bands of Landsat images for further processing [28].We calculated the normalized difference vegetation index (NDVI) [36], normalized difference built-up index (NDBI) [37], and modified normalized difference water index (MNDWI) [38] to enhance the information on vegetation, impervious surface, and water [28].These calculated bands were combined with the original six bands.
We used greenest pixel compositing method [27] and percentile-based image compositing method [39] to create composited images.The greenest pixel compositing method considers associated NDVI band of each Landsat scene and selects band values from Landsat scenes that have the highest NDVI values [27].This procedure generated a nine-band composite image for each city.For each pixel, values of NDVI and other eight bands were extracted from the same image.The percentile-based image compositing method is also a pixel-based compositing method [28].This method first ranks all available pixel values from minimum to maximum for each specific pixel location.Then values for specific percentile can be selected.For example, value of percentile 0 represents the minimum value for a specific pixel from all available pixel values.In this study, we chose percentile values at 0%, 10%, 25%, 50%, 75%, 90%, and 100% [39] for each band.This step generated a 63-band (nine bands multiplied by seven percentile values) composite image for each city.
Images generated using the greenest pixel compositing method reflect the status of the peak growing season and were used to separate vegetation from nonvegetation.The percentile-based compositing method can capture phenology information without any explicit assumption and prior knowledge of the timing of phenology [40].It was used to differentiate different vegetation types such as crops and grasses.

Land Cover Classification
To collect training samples, we first divided the 262 cities into three city groups according to the 'urban ecoregions' proposed by Schneider et al. (2010) [41].The three city groups were Arid, semi-arid steppe in Central Asia, temperate forest in East Asia, and tropical, subtropical forest in Asia.Then we collected representative training samples of six types of land covers (i.e., water, impervious surface, bare soil, tree/shrub, grass, and crop) by referencing to high-resolution images on Google Earth (GE) [28].We collected 500 to 5000 polygons for each class in each city group.In total, we collected 41,281 polygons as training samples.Specifically, 4900 polygons of bare soil, 6401 polygons of crop, 4319 polygons of grass, 13,073 polygons of impervious surface, 9838 polygons of tree/shrub, and 2750 polygons of water samples were collected.
We performed land cover classification on the GEE platform using the composited images, training samples, and a Random Forest (RF) classifier [42].We first classified the percentile-based composited images into water and nonwater type.Then we classified the nonwater type into vegetation and nonvegetation type using the greenest pixel composited images.We further classified the vegetation class into tree/shrub, grass, and crop, and classified the nonvegetation type into impervious surface, and bare soil using the percentile-based composited images.The number of trees for RF was set to 100.
We conducted accuracy assessment by extracting 100 pixels for each city from the classified land cover maps as validation samples by using a stratified random sampling method [43].The land cover classes were used as strata and the number of pixels sampled for each stratum was in proportion to the percentage of each land cover class [28].We manually interpreted these validation samples by referencing to high-resolution images on GE and estimated classification accuracies based on the interpretation results.

UGS Extraction
Based on the land cover maps, we generated borders of urban areas for each city.We first created 1 km by 1 km impervious surface intensity (0-1) maps using the Block Statistics tool in ArcMap 10.1 (Esri Inc., Redlands, CA, USA).Next we divided the intensity maps into built/nonbuilt maps using a threshold value of 0.2 [44].Then we converted the built/nonbuilt maps to polygons and chose the largest polygon as the core urban area of each city.We filled holes within the largest polygon using Eliminate tool in ArcMap 10.1 to form the continuous boundary of urban area.We extracted out tree/shrub and grass cover within the boundary of urban area and merged them as UGS.

Calculation of Landscape Structure Metrics
To evaluate landscape structure of UGS for each city, we calculated five landscape metrics commonly used in describing landscape structure of UGS [17,19,20,45,46].The five landscape metrics were percentage of landscape (PLAND), patch density (PD), largest patch index (LPI), mean Euclidian nearest neighbor distance (ENN_MN), and mean patch shape index (SHAPE_MN) (Table 1).PLAND was used to quantify abundance of UGS.Cities with larger PLAND values have more UGS.LPI was used to measure dominance of UGS in the landscape.Large LPI values mean that the largest urban green space patch is more dominant in cities. PD was used to quantify the fragmentation of UGS.Cities with larger PD values have more fragmented UGS.ENN_MN was used to quantify isolation of urban green space patches.Cities with larger ENN_MN values have more dispersed UGS.SHAPE_MN was used to measure complexity of patch shapes of UGS.Cities with higher SHAPE_MN values have more complex UGS patch shapes.We calculated the landscape metrics using Fragstats 4.2 and the extracted maps of UGS [47].We examined the correlations between these metrics using the Pearson correlation test and found they were not highly correlated with each other (R < 0.6) (Table S2).

Calculation of Urban form Metrics
We chose three commonly used urban form metrics: compound terrain complexity index (CTCI) [48,49], road density (RD) [50], and perimeter-area ratio (PARA) [45].CTCI indicates terrain complexity of each city.It was calculated using the 30-m Shuttle Radar Topography Mission (SRTM) data [51].CTCI was calculated as follows.
where N SE , N RE , N TC and N RU represent normalized standard deviation of elevation, normalized range of elevation, normalized terrain curvature, and normalized rugosity.The normalized index (N SE , N RE , N TC , N RU and N RU ) was calculated as follows, RD represents density dimension of urban form for each city, which was calculated using the road vector data downloaded from the OpenStreetMap (OSM) website (www.openstreetmap.org).The OSM is developed as a volunteered geographic information project thus the data set may include uncertainties.However, it is the most comprehensive data set on urban roads around the world that is freely available to us.RD was calculated as follows.
where L represents total length of roads, A urban represents the size of urban area.PARA represents shape complexity of each urban area, which was calculated using the vector boundary of continuous urban area.

PARA =
P urban A urban (4) where P urban represents perimeter of an urban area, A urban represents total area of an urban area.

Statistical Analysis
Associations between landscape metrics and urban form metrics were explored using the maximal information-based nonparametric exploration (MINE) method [52].We calculated the maximum information coefficients (MIC), Pearson correlation coefficients (ρ), and the nonlinearity coefficients (MIC-ρ 2 ) between each landscape metric and the urban form metrics using the Java version of the MINE application [52].The MIC values are between 0 and 1, with values near 0 representing statistically independent relationship and values near 1 representing noiseless functional relationships [52].The MIC-ρ 2 values range from 0 to 1, with values near 0 representing linear relationships and large (usually >0.2) values representing nonlinear relationships [52].
We used the boosted regression tree (BRT) model to analyze the relative importance of urban form metrics on each landscape metric.We chose the BRT model because it can handle nonlinearities and different type of variables.The BRT model was fitted using the gbm.step function of dismo package with R (version 3.3.1;R Development Core Team, Vienna, Austria).Model parameters of the learning rate, tree complexity, and bag fraction were set to 0.001, 6, and 0.75, respectively.

Landscape Structure of UGS in 262 Chinese Cities
The mean overall classification accuracy of land cover was 88.76 ± 3.17%.The mean kappa coefficient was 0.84 ± 0.05.Confusion matrix of land cover classification for each city was listed in Supplementary file 2.
Landscape metrics of UGS for the 262 Chinese cities were shown in Figure 2. Values of the landscape metrics varied among these cities.The ENN_MN values varied from 70.73 m to 130.53 m, averaging 82.74 m (SD = 6.54 m), which showed varied dispersion of UGS patches in these cities.The mean LPI was 4.34% (SD = 3.62%), ranging from 0.22% to 21.18%.The values suggested that the dominance of the largest UGS patches is weak in most cities and variation is high among these cities.The mean SHAPE_MN was 1.32 (SD = 0.04), ranging from 1.20 to 1.46.The shape complexity of UGS did not vary much across cities.The mean PD was 11.49 count/km 2 (SD = 2.48 count/km 2 ), ranging from 5.21 to 19.06 count/km 2 .This wide range of PD values indicated that levels of fragmentation of UGS varied greatly in studied cities.The mean PLAND was 30.46% (SD = 8.11%), ranging from 6.17% to 57.05% (Figure 2).Clearly, abundance of UGS varied substantially among these cities.
Remote Sens. 2018, 10, x FOR PEER REVIEW 6 of 14 We used the boosted regression tree (BRT) model to analyze the relative importance of urban form metrics on each landscape metric.We chose the BRT model because it can handle nonlinearities and different type of variables.The BRT model was fitted using the gbm.step function of dismo package with R (version 3.3.1;R Development Core Team, Vienna, Austria).Model parameters of the learning rate, tree complexity, and bag fraction were set to 0.001, 6, and 0.75, respectively.

Landscape Structure of UGS in 262 Chinese Cities
The mean overall classification accuracy of land cover was 88.76 ± 3.17%.The mean kappa coefficient was 0.84 ± 0.05.Confusion matrix of land cover classification for each city was listed in Supplementary file 2.
Landscape metrics of UGS for the 262 Chinese cities were shown in Figure 2. Values of the landscape metrics varied among these cities.The ENN_MN values varied from 70.73 m to 130.53 m, averaging 82.74 m (SD = 6.54 m), which showed varied dispersion of UGS patches in these cities.The mean LPI was 4.34% (SD = 3.62%), ranging from 0.22% to 21.18%.The values suggested that the dominance of the largest UGS patches is weak in most cities and variation is high among these cities.The mean SHAPE_MN was 1.32 (SD = 0.04), ranging from 1.20 to 1.46.The shape complexity of UGS did not vary much across cities.The mean PD was 11.49 count/km 2 (SD = 2.48 count/km 2 ), ranging from 5.21 to 19.06 count/km 2 .This wide range of PD values indicated that levels of fragmentation of UGS varied greatly in studied cities.The mean PLAND was 30.46% (SD = 8.11%), ranging from 6.17% to 57.05% (Figure 2).Clearly, abundance of UGS varied substantially among these cities.   Figure 3 shows the spatial distribution of PLAND values for the 262 Chinese cities. Cities in Southwest and Central South of China had high PLAND values, which indicated that cities in these regions had high UGS coverage.However, cities in Northwestern China and North China had lower PLAND values than many other cities, which indicated that these cities had low UGS coverage (Figure 3).Spatial distributions of other four landscape metrics can be found in the supplementary file (Figures S1-S4).

Landscape Structure of UGS in China
Our results revealed features of landscape structure of UGS in Chinese cities at the national scale.Dobbs et al. (2017) [23] investigated landscape structures of UGS for 100 cities worldwide and found that density of vegetation patch, nearest neighbor, and green cover of these cities were 0.33 count/ha (SD = 0.13 count/ha), 253.5 m (SD = 61.9m), and 32.6% (SD = 12.1%), respectively.The percentage of UGS (i.e., PLAND) of Chinese cities was similar to world cities reported by Dobbs et al. (2017) [23].However, UGS in Chinese cities were more fragmented (higher PD) and less isolated (lower ENN_MN).This agrees with a case study of nine major Chinese cities in which they found UGS of the nine cities were highly fragmented and heterogeneous [21].This might be due to the fact that many Chinese cities have high population density (e.g., Hong Kong, Guangzhou, and Shanghai [53]), thus there are fewer spaces for large and homogeneous UGS to exist.

Impacts of Urban forms on Landscape Structure of UGS
We found that all three urban form metrics had significant effects on landscape structure of UGS.Cities with higher RD tended to have less UGS (lower PLAND value).This might because high RD is a feature of compact cities, which may have limited spaces for UGS [54].Our results are in agreement with the study in Sheffield, UK [30] in which they found a negative relationship between extent of green space and total length of the roads network.We found that cities with higher RD tended to have more fragmented UGS (higher PD value).Urban landscape is often highly dissected by roads [55] and has fragmented UGS as a result.
Terrain complexity is a measure of terrain heterogeneity [49], which may positively affect the spatial distribution of UGS [30].We found that CTCI affected landscape structure of UGS significantly.Cities with high terrain complexity tended to have more UGS (higher PLAND value), but UGS in these cities tended to be more fragmented (higher PD value).This observation was in agreement with the situation in Cincinnati, USA, in which the authors found that terrain roughness was positively correlated with the levels of tree canopy [31].A study in 60 urban areas of Central Indiana, USA also found a positive relationship between urban tree canopy cover and slope of urban areas [56].
PARA indicates shape complexity of urban areas.We found that cities with complex shapes (higher PARA value) also had more fragmented UGS (higher PD value).In the USA, megapolitan areas with more complex shapes were found to have more complex and fragmented urban land covers than metropolitan areas with less complex shapes [57].The same pattern was observed in Chinese cities through our study.

Implications for UGS Planning
Our study can provide some useful information for planning and management of UGS in China and other regions and countries facing rapid urbanization.Our results showed that urban form has a significant influence on landscape structure of UGS and must be considered in UGS planning.First, cities with high RD tended to have less and fragmented UGS.China is experiencing a rapid urbanization process and new roads are being built up consistently.When planning for new roads in Chinese cities, especially in cities with high RD, planners should pay attention to the existing UGS patches to avoid further fragmentation of UGS.Furthermore, to construct roadside greenway networks whenever possible may be a way to alleviate the negative effects associated with fragmentation caused by roads on UGS in cities with high RD.Roadside greenway networks have been found to act as important corridors for urban residents and wildlife [58].In addition, innovative designs such as vegetated wildlife overpasses [59] and elevated linear parks [60] can be used to alleviate the fragmentation effects caused by roads.Besides, a vegetated drainage system can also be designed as a supplement to reverse habitat fragmentation and act as wildlife corridors and buffer zones [61].
Second, we found that cities with complex shapes often have more fragmented UGS.Urban planners should pay attention to the restriction caused by urban form when planning and designing the UGS system.Finally, our results showed that cities with high terrain complexities have more UGS but that these UGS are also more fragmented.There are both advantages and disadvantages when planning for UGS in cities with complex terrains.On the one hand, complex terrain may provide an opportunity to conserve green spaces in cities as it is difficult to develop on sites with complex terrain.On the other hand, planners must realize that corridors among green spaces are important for sustaining urban biodiversity [16] and should be preserved as much as possible.This task may be more difficult in cities with complex terrains because some sites are more isolated by the terrain.We call for actions to preserve natural pathways and corridors among green spaces in these cities during the development.The increased connectivity of green spaces can improve their ecological functions, such as their use as habitats by wildlife [62].
While our findings provide new knowledge on the influences of urban form on landscape structure of UGS, they should be interpreted with the limitations of the current study in mind.We only used a set of basic urban form and landscape metrics due to data availability.To use more metrics in future studies may reveal more details on the relationship between urban form and landscape structure of UGS.In addition, we only focused on cities with medium size and above where UGS are important due to their population sizes.The findings from our study may not be readily applicable to small cities.Despite these limitations, we showed that urban form has a significant influence on landscape structure of UGS using data retrieved from satellite images.The finding can contribute to better planning and management of UGS.

Conclusions
Landscape structure of UGS is closely related to urban ecosystem services [1,2].Quantifying landscape structure of UGS and analyzing its influencing factors are important for managing ecosystem services for urban inhabitants' health and wellbeing.We assessed the influences of urban form on landscape structure of UGS of 262 Chinese cities using Landsat images and the GEE platform.The GEE platform allowed us to obtain, preprocess, and classify 6673 scenes of Landsat images for 262 cities in a relatively short time.With the help of this geospatial analysis platform, we revealed the unique features of landscape structure of UGS in Chinese cities.We also showed that urban form had significant influences on landscape structure of UGS.Based on our results, we suggested that attention should be paid to urban form when planning and designing UGS.More urban form and landscape metrics can be considered in future studies to reveal more details of this relationship.

Figure 1 .
Figure 1.Locations of the 262 studied Chinese cities.

Figure 2 .
Figure 2. Summary statistics of landscape metrics of urban green spaces, (a) percentage of landscape, (b) patch density, (c) largest patch index, (d) mean patch shape index, and (e) mean Euclidian nearest neighbor distance.The box shows median, 25th and 75th percentile, and the whiskers illustrate the 5th and 95th percentile values.

Figure 3
Figure3shows the spatial distribution of PLAND values for the 262 Chinese cities. Cities in Southwest and Central South of China had high PLAND values, which indicated that cities in these regions had high UGS coverage.However, cities in Northwestern China and North China had lower PLAND values than many other cities, which indicated that these cities had low UGS

Figure 2 .
Figure 2. Summary statistics of landscape metrics of urban green spaces, (a) percentage of landscape, (b) patch density, (c) largest patch index, (d) mean patch shape index, and (e) mean Euclidian nearest neighbor distance.The box shows median, 25th and 75th percentile, and the whiskers illustrate the 5th and 95th percentile values.

Figure 5 .
Figure 5. Partial dependence plots for boosted regression tree analysis showing the effects of urban form predictors on: (a) ENN_MN (mean Euclidian nearest neighbor distance), (b) LPI (largest patch index), (c) SHAPE_MN (mean patch shape index), (d) PD (patch density), and (e) PLAND (percentage of landscape).Numbers enclosed inside parenthesis indicated the relative importance of predictors.CTCI, compound terrain complexity index; PARA, perimeter-area ratio; RD, road density.

Table 2 .
Associations between landscape metrics and the urban form metrics.