Spatial Characteristics of Urban Green Spaces and Human Health: An Exploratory Analysis of Canonical Correlation

In highly urbanized areas, urban green spaces (UGSs) are important natural and cultural entities. Previous studies have shown some evidence of positive relationships between UGSs and human health. Most of these studies relied on self-reported health data and often used institutional quantitative measures of UGSs instead of the spatial characteristics of UGSs. This study analyzed the relationships between the spatial characteristics of UGSs and morbidity of diseases, which were considered variables of human health in the Taipei Metro. The Longitudinal Health Insurance Database of Taiwan was applied as the source of morbidity of diseases. A canonical correlation analysis was performed by using the six variables of spatial characteristics of UGSs as predictors and three variables of morbidity as criterion variables to evaluate the multivariate shared relationships between the two variable sets. The results found a strong canonical correlation between the spatial characteristics of UGSs and human health. Furthermore, the results revealed that living in districts with a higher area percentage of green spaces and denser vegetation cover, as well as exposure to more aggregative and irregular-shape green spaces, can reduce the morbidity of diseases.


Introduction
Green space can be defined as lands that consist predominantly of unsealed, pervious, soft surfaces, such as soil, grass, shrubs, and trees, with few structures present [1]. In highly urbanized areas sealed by artificial constructions, green spaces are important natural and cultural entities for urban residents. Green spaces in urban areas provide important ecosystem services, such as air quality improvement, atmospheric carbon dioxide reduction, and recreational and cultural values [2][3][4]. However, urbanization, environmental degradation, and lifestyle changes are diminishing opportunities for human contact with nature [5].
Amenities in green spaces encourage urban residents to increase physical activities, to pursue more leisure and recreation, and to have contact with the natural environment. Increasing evidence indicates that the presence of green spaces can enhance the quality of life [6,7]. For example, easy access to green spaces increases physical activity [5,8,9] and promotes physical and mental health of urban residents [6,[10][11][12][13]. More greenness in urban residential areas is associated with more physical activity and lower risk of disease mortality [14]. Moreover, urban green spaces (UGSs) can facilitate social interactions, social networking, social support, social inclusion, and improve the mental health of residents [15][16][17].

Study Area
In this study, the largest metropolitan area in Taiwan, Taipei Metro Area, was selected for the analysis of the relationships between spatial patterns of green spaces and human health. Figure 1 shows that Taipei Metro Area has three cities, i.e., Taipei City, New Taipei City, and Keelung City, and 48 districts and has a total area of 2457 km 2 , with a population of 7.04 million at the end of 2017, representing 30% of the total population of Taiwan [39]. The mean annual temperature of Taipei Metro Area is 23 • C, with the variation between the lowest monthly average temperature of 16 • C in January and the highest monthly average temperature of 29.6 • C in July. The annual rainfall varies by terrain. The rainfall approximates 2400 mm in plain areas and is over 4500 mm in mountain areas (The information of rainfall and temperature data about the Taipei Metro Area was collected from [40][41][42]) This study analyzed UGSs in 37 districts that had a population larger than 25,000 and population density above 500 pop/km 2 . The remote mountain districts and rural districts with lower population or population density were excluded. section concludes with the significant correlations between UGSs and human health and the advantages of using the data derived from the national health insurance system.

Study Area
In this study, the largest metropolitan area in Taiwan, Taipei Metro Area, was selected for the analysis of the relationships between spatial patterns of green spaces and human health. Figure 1 shows that Taipei Metro Area has three cities, i.e., Taipei City, New Taipei City, and Keelung City, and 48 districts and has a total area of 2457 km 2 , with a population of 7.04 million at the end of 2017, representing 30% of the total population of Taiwan [39]. The mean annual temperature of Taipei Metro Area is 23 °C, with the variation between the lowest monthly average temperature of 16 °C in January and the highest monthly average temperature of 29.6 °C in July. The annual rainfall varies by terrain. The rainfall approximates 2400 mm in plain areas and is over 4500 mm in mountain areas (The information of rainfall and temperature data about the Taipei Metro Area was collected from [40][41][42]) This study analyzed UGSs in 37 districts that had a population larger than 25,000 and population density above 500 pop/km 2 . The remote mountain districts and rural districts with lower population or population density were excluded.

Spatial Characteristics of UGSs
The digital map of UGSs of the study area was derived from the Taiwan land use investigation maps for 2015 produced by the National Land Surveying and Mapping Center, Ministry of the Interior. UGSs were selected by the attributes of land use/land cover and reclassified into five types: farmland, forest, wetland, grassland, and park and open spaces. This study also selected residential land use as the focus area to assess the spatial accessibility to UGSs. The ArcGIS tool was applied to convert the map from vector format to raster format with a pixel size of 30 m × 30 m for metric analysis using the tool for measuring landscape metrics, FRAGSTATS 4.2. Satellite data of SPOT 5 on 8 May 2014 and 5 February 2012 were also used for analyzing the normalized difference vegetation index (NDVI) of UGSs by using ERDAS IMAGINE.
A large set of metrics for analyzing landscape composition and configuration has been developed, modified, and tested in the past three decades [42,43]. In this study, we targeted the landscape metrics of UGSs regarding human health based on the literature review. The proportion

Spatial Characteristics of UGSs
The digital map of UGSs of the study area was derived from the Taiwan land use investigation maps for 2015 produced by the National Land Surveying and Mapping Center, Ministry of the Interior. UGSs were selected by the attributes of land use/land cover and reclassified into five types: farmland, forest, wetland, grassland, and park and open spaces. This study also selected residential land use as the focus area to assess the spatial accessibility to UGSs. The ArcGIS tool was applied to convert the map from vector format to raster format with a pixel size of 30 m × 30 m for metric analysis using the tool for measuring landscape metrics, FRAGSTATS 4.2. Satellite data of SPOT 5 on 8 May 2014 and 5 February 2012 were also used for analyzing the normalized difference vegetation index (NDVI) of UGSs by using ERDAS IMAGINE.
A large set of metrics for analyzing landscape composition and configuration has been developed, modified, and tested in the past three decades [42,43]. In this study, we targeted the landscape metrics of UGSs regarding human health based on the literature review. The proportion of green spaces (PLAND) was calculated to assess the availability of UGSs in the Taipei Metro area. This study also selected the shape of UGS (SHAPE_MN), mean distance of UGS patches to the nearest residential patch (ENN_MN), and continuity of UGS patches (COHESION) as the indicators of the spatial pattern of UGSs. All these landscape metrics were implemented by using the statistical package FRAGSTATS 4.2 [44] (See Table 1). The quality of UGSs was measured by using NDVI to assess the greenness of vegetation in the Taipei Metro Area. The index of NDVI correlates with biophysical properties of the vegetation canopy, such as leaf area index, fractional vegetation cover, vegetation condition, and biomass [46]. NDVI is often used as a metric to assess exposure to surrounding greenness. Previous studies also showed that a higher NDVI value is related to less depression [47] and better mental health [48]. Additionally, the analysis of Network Distance of ArcGIS was used to calculate the shortest network distance from the residential area to UGS as the index of accessibility.
The landscape spatial characteristics of UGSs were analyzed for 37 districts of Taipei Metro Area. The results of the spatial analysis were used to examine the correlation between spatial characteristics of UGSs and morbidities of diseases in the study area.

National Health Insurance Research Database
The database of Taiwan's National Health Insurance program contains registration files and original claim data for reimbursement. The National Health Insurance Administration, Ministry of Health and Welfare, Taiwan, derives large computerized databases from this program. The databases are maintained by the National Health Research Institutes, Taiwan, and are provided to scientists for research purposes. The data files from the National Health Insurance program are de-identified by scrambling the identification codes of both patients and medical facilities and sent to the National Health Research Institutes to form the original files of the National Health Insurance Research Database (NHIRD). This study applied the Longitudinal Health Insurance Database (LHID), which is a data subset of the NHIRD and contains all original claim data for 1,000,000 beneficiaries randomly sampled from the Registry for Beneficiaries (ID) of the NHIRD. All the registration files and original claim data of these 1,000,000 individuals were collected and constituted the LHID [49].
In Taiwan, according to the Personal Information Protection Act, residential locations of patients are confidential. Dataset of outpatients reveals only the spatial locations of the hospital. To consider the influences of cross-district outpatients, we combined the outpatient ID of the common cold (ICD9 460) and influenza (ICD9 487) with the outpatient ID of cardiovascular disease, mental disorders, and chronic respiratory disease. The locations of medical institutions for common cold and influenza are determined as the living places of outpatients based on the combined data.
This study applied the ambulatory care expenditures by visits (CD) of original claim data as the source of morbidity of diseases. The CD data contains the date of ambulatory care, ICD-9-CM codes, and ID of medical facilities (HOSP_ID). By combining ID of medical facilities (HOSP_ID) and registry for contracted medical facilities (HOSB) which includes area code of medical facilities (AREA_NO_H), the number of outpatients can be summed up by district level. We included all registered outpatient database between 1 January 2010 and 31 December 2012 and selected the database of outpatients for the three types of disease related to exposure of green space, i.e., cardiovascular diseases (ICD9 390-459), mental disorders (ICD9 290-319), and chronic respiratory diseases (ICD9 470-478, 490-519).

Canonical Correlation Analysis
In 1936, Hotelling developed canonical correlation analysis (CCA) to examine the correlation between a synthetic criterion variable and a synthetic predictor variable that are weighted based on the relationships between the variables within the sets [50]. Instead of analyzing a single dataset, however, the goal of CCA is to analyze the relationship between a pair of datasets [51]. Each dataset can contain several variables, and CCA calculates a linear combination for each set, called a canonical variable, such that any correlation between two canonical variables is maximized [52]. These correlations are called canonical correlations, and the linear combinations are canonical variates. The fundamental principle behind CCA is to create a number of canonical variates, each consisting of a linear combination of one set of variables (X i ) in the following equation: Additionally, a linear combination of the other set of variables (Y i ) has the following equation: The goal of CCA is to determine the coefficients, or canonical weights (a ij and b ij ), that maximize the correlation between canonical variates U i and V i . The first canonical correlation, Corr. (U 1 , V 1 ), is the strongest possible correlation between a linear combination of variables in the exposure set and a linear combination of variables in the outcome set. Other pairs of maximally correlated linear combinations are chosen in turn and are orthogonal to those previously identified. The process of constructing canonical variates continues until the number of pairs of canonical variates equals the number of variables in the smaller of the two sets.
Because CCA assesses the existence of correlations between two sets of variables, this exploratory statistical method is useful for analyzing relationships between the combinations of urban form and different variables of human health. Although CCA is widely used in social science research [53], the method is rarely used in studies on UGSs [54]. This study applied CCA to examine multivariate relationships between spatial characteristics of UGSs and human health. The proposed canonical correlation model is shown in Figure 2. This study attempted to quantify the independent relationships of spatial characteristics of UGSs and levels of morbidity of three diseases. Table 2 shows the descriptive statistics of the studied variables. It provides information about the number of samples, minimum, maximum, mean, and standard deviation of the studied variables.
A canonical correlation analysis was performed with IBM SPSS Statistics for Windows, Version 22.0 (IBM Corp. released 2013, Armonk, NY, USA). Six variables of spatial characteristics of UGSs were used as predictors of three observed variables of morbidity to evaluate the multivariate shared relationship between these two variable sets. Before the analysis results show the fitness of the model, the Pearson correlation matrix describes the associations between the variables of the spatial characteristics of UGSs. It also detects the possibility of multicollinearity of studied variables. Table 3 shows that all the correlations (r) < 0.765 (in absolute value). The general rule of thumb is that if the correlation (r) > 0.8, then severe multicollinearity may be present. Therefore, it is not necessary to remove the redundant variables.
This study uses a multivariate test of significance of CCA (i.e., Pillai's, Hotelling's, Wilk's, and Roy's multivariate criteria) to show the general fit of the model. Furthermore, the canonical correlation coefficients and the eigenvalues are used to explain the fitness of the canonical roots of the model. Based on the tests, the significance of each of the roots in the dimension reduction analysis of SPSS, the possible significant roots are found. Since our model contains three disease morbidities and six indicators of spatial characteristics of UGSs, SPSS extracts three canonical roots or dimensions. In the dimension reduction analysis, the first test of significance tests all three canonical roots of significance; the second test excludes the first root and tests roots two to three; the last test tests root three by itself. . The canonical correlation model that was used in our study to assess the relationship of the predictor variables (spatial characteristics of urban green spaces (UGSs)) with criterion variables (morbidities of CVD (cardiovascular diseases), mental disorders, and respiratory diseases). A canonical correlation analysis was performed with IBM SPSS Statistics for Windows, Version 22.0 (IBM Corp. released 2013, Armonk, NY, USA). Six variables of spatial characteristics of UGSs were used as predictors of three observed variables of morbidity to evaluate the multivariate shared relationship between these two variable sets. Before the analysis results show the fitness of the model, the Pearson correlation matrix describes the associations between the variables of the spatial characteristics of UGSs. It also detects the possibility of multicollinearity of studied variables. Table  3 shows that all the correlations (r) < 0.765 (in absolute value). The general rule of thumb is that if the correlation (r) > 0.8, then severe multicollinearity may be present. Therefore, it is not necessary to remove the redundant variables. Figure 2. The canonical correlation model that was used in our study to assess the relationship of the predictor variables (spatial characteristics of urban green spaces (UGSs)) with criterion variables (morbidities of CVD (cardiovascular diseases), mental disorders, and respiratory diseases).

Spatial Pattern of the Characteristics of UGSs and the Morbidity of Diseases
We analyzed the spatial characteristics of UGSs of the districts in the study area based on the land use investigation map and the FRAGSTATS 4.2 software. Figure 3a shows that the districts with a low percentage of green spaces roughly clustered in the center of the Taipei Metro Area. The level of NDVI of the districts (Figure 3b) shows that UGS quality was lower in the urban central districts than that in peripheral districts. The level SHAPE_MN of the districts (Figure 3c) similarly shows that the shape of UGSs in the urban central districts was more regular compared to peripheral districts. In highly urbanized districts, the distances between green spaces (ENN_MN) were longer than those in the lightly urbanized out-ring districts (Figure 3d). Conversely, the continuity of green spaces (COHESION) was higher in peripheral districts compared to central districts (Figure 3e). The Network Distances from residential areas to green spaces were longer in the central districts (Figure 3f). Generally, availability, quality, and accessibility of UGS were higher in southern and northern west parts of the Taipei Metro Area. Inversely, availability, quality, and accessibility of UGSs of the districts in the Taipei basin area and Keelung city were poorer. The connectivity of UGS was lower in highly urbanized districts.
We analyzed the spatial characteristics of UGSs of the districts in the study area based on the land use investigation map and the FRAGSTATS 4.2 software. Figure 3a shows that the districts with a low percentage of green spaces roughly clustered in the center of the Taipei Metro Area. The level of NDVI of the districts (Figure 3b) shows that UGS quality was lower in the urban central districts than that in peripheral districts. The level SHAPE_MN of the districts (Figure 3c) similarly shows that the shape of UGSs in the urban central districts was more regular compared to peripheral districts. In highly urbanized districts, the distances between green spaces (ENN_MN) were longer than those in the lightly urbanized out-ring districts (Figure 3d). Conversely, the continuity of green spaces (COHESION) was higher in peripheral districts compared to central districts (Figure 3e). The Network Distances from residential areas to green spaces were longer in the central districts ( Figure  3f). Generally, availability, quality, and accessibility of UGS were higher in southern and northern west parts of the Taipei Metro Area. Inversely, availability, quality, and accessibility of UGSs of the districts in the Taipei basin area and Keelung city were poorer. The connectivity of UGS was lower in highly urbanized districts. Morbidities of cardiovascular disease, mental disorder, and chronic respiratory disease were high, and the quality of green spaces was low in the central districts of the Taipei Basin. Conversely, western and southern districts had low morbidity and high quality and connectivity of green spaces. The spatial patterns of morbidity of cardiovascular disease, mental disorder, and chronic respiratory Morbidities of cardiovascular disease, mental disorder, and chronic respiratory disease were high, and the quality of green spaces was low in the central districts of the Taipei Basin. Conversely, western and southern districts had low morbidity and high quality and connectivity of green spaces. The spatial patterns of morbidity of cardiovascular disease, mental disorder, and chronic respiratory disease ( Figure 4) were roughly opposite to the quality, greenness, and connectivity of green spaces in our study area.
(e) (f) Morbidities of cardiovascular disease, mental disorder, and chronic respiratory disease were high, and the quality of green spaces was low in the central districts of the Taipei Basin. Conversely, western and southern districts had low morbidity and high quality and connectivity of green spaces. The spatial patterns of morbidity of cardiovascular disease, mental disorder, and chronic respiratory disease ( Figure 4) were roughly opposite to the quality, greenness, and connectivity of green spaces in our study area.

Relationships Between Spatial Characteristics of UGSs and Morbidity
The output of CCA shows the general fit of the model reporting Pillai's, Helling's, Wilk's and Roy's multivariate criteria. The commonly used test is Wilk's, but we found that all of these tests were significant, with p < 0.001 (see Table 4).

Relationships Between Spatial Characteristics of UGSs and Morbidity
The output of CCA shows the general fit of the model reporting Pillai's, Helling's, Wilk's and Roy's multivariate criteria. The commonly used test is Wilk's, but we found that all of these tests were significant, with p < 0.001 (see Table 4).  Table 5 reports the canonical correlation coefficients and the eigenvalues of the canonical roots. The CCA analysis yielded three canonical functions with squared canonical correlation (Sq Cor.) of 0.62941, 0.40733, and 0.09710 for each root. The first pair of variates, a linear combination of spatial characteristics of UGSs and a linear combination of the variables of morbidities of diseases, had a canonical correlation coefficient of 0.79336 and an eigenvalue of 1.698. The analytical results indicate that the spatial characteristics of UGSs and human health are correlated. The first root was created to maximize the canonical correlation (Pearson r) between the two synthetic variables. The second pair of canonical variates was constructed out of the residuals of the first pair to maximize the correlation between them. The first and second roots had explained variances of 68.12% and 27.57%. The third root had a low explained variance of 4.31% (see Table 5). SPSS extracted three canonical roots or dimensions for our model. The dimension reduction analysis tested the significance of each of the roots (see Table 6). The first test of significance showed the full model across root 1 to 3 was significant (Wilk's λ = 0.19831, F = 3.41680, p < 0.01). The second test excluded the first root and tests roots 2 to 3 and showed it was also significant (Wilk's λ = 0.53512, F = 2.21872, p = 0.036 < 0.05). The root 3 itself was not significant (Wilk's λ = 0.90290, F = 0.80659, p = 0.531) (see Table 6). Next, we turned to the research question of which variables of the spatial characteristic of UGSs contribute to this relationship between the variable sets across the two functions. Table 7 shows the statistics of the canonical solution of Function 1 and Function 2 for the relationships between spatial characteristics of UGSs and morbidities of diseases. The table presents the structure coefficients (r s ), the squared structure coefficients (r s 2 ), and cross-structure coefficients (r c ) between a set of variables and the synthetic variable created by another set of variables. In Function 1, the first pair of canonical variates groups the variables in the way that the correlation between them is maximized. The structure coefficients (r s ) of PLAND, NDVI, SHAPE_MN, ENN_MN, and Network Distance indicated that they were the primary contributors (with r s = −0.53151, −0.60826, −0.48613, 0.53458, and 0.88403) to the synthetic variate of characteristics of green spaces. The crossstructure coefficients (r c ) of PLAND, NDVI, and SHAPE_MN showed that they were negatively related to the synthetic variate of morbidity of diseases. Inversely, the cross-structure coefficients (r c ) for ENN_MN and Network Distance show that they were positively related to the synthetic variate of morbidity of diseases. The other side of Function 1, the structure coefficients (r s ) of criterion variables indicate that cardiovascular diseases, mental disorders, and respiratory diseases were positively correlated with the synthetic variate of morbidity (with very high r s = 0.98556, 0.98352, and 0.99449). The cross-structure coefficients (r c ) of the three criterion variables, all of which were higher than 0.7, revealed that morbidity of cardiovascular diseases, mental disorders, and respiratory diseases all had positive correlations with the synthetic variate of spatial characteristics of UGSs.
In the Function 2, the values of structure coefficients (r s ) indicated that the variables of PLAND, NDVI, and COHESION were the primary contributors to the synthetic variate of predictors. The crossstructure coefficients (r c ) of PLAND, NDVI, and COHESION showed that they were positively related to the synthetic variate of disease morbidities with moderate r c = 0.34186, 0.33197, and 0.47943. On the other side of Function 2, the structure coefficients (r c ) of criterion variables indicated that cardiovascular diseases, mental disorders, and respiratory diseases had very low correlation coefficients with the synthetic variate of morbidity of diseases (with very low r s < 0.2). The cross-structure coefficients (r c ) of all the three criterion variables were lower than 0.12, revealing that morbidity of cardiovascular diseases, mental disorders, and respiratory diseases all had very low correlations with the synthetic variate of spatial characteristics of UGSs.

Discussion
In this study, we correlated the district-level spatial characteristics of UGSs with the morbidities of specific diseases. The results of CCA analysis showed a significant correlation between spatial characteristics of UGSs and human health in the Taipei Metro Area. The cross-structure coefficients of predictor variables revealed that the morbidity of diseases was decreased by residence in the districts with higher area percentage, denser vegetable cover, more aggregative urban greens, and irregular-shape green spaces. In comparison with the variables of spatial characteristics of UGSs, accessibility (measured by Network Distance) was the major contributor to the relationship between spatial characteristics and human health. The secondary contributors were the variables of quality of UGSs (measured by NDVI) and spatial pattern of UGSs (measured by SHAPE_MN, ENN_MN, and COHESION). The availability of UGSs (measured by PLAND) was less influential in the relationship between spatial characteristics of UGSs and human health.
The results of correlation analyses validated the evidences of health benefits of exposure to green spaces that have been found by previous studies (e.g., [21,35,[55][56][57]). This study also revealed the key variables of spatial characteristics of UGSs in the relationship between UGSs and human health. The analytical results showed that the variables of quality (measured by NDVI) and accessibility played key roles in influencing the morbidities of diseases in the Taipei Metro Area. These findings confirm the evidence found by [9,12,34] in their study areas. In this study, the strategy of using a database of national health insurance and conducting district-level analyses is similar to the newly published study of [29]. However, the previous study usually concerned only with a single disease, i.e., cardiovascular disease. Our study reveals the correlations between the spatial characteristics of UGSs and the morbidity of three diseases.
A major limitation of our study is that the data of the National Health Insurance Research did not reveal the real addresses of the patients but only the districts where the patients were insured and outpatient clinics located. However, the locations of medical institutions they visited sometimes may be different from that they lived. To reduce the effects of cross-district outpatients, we cross-matched the address of the medical institution which the patients visited for the common cold and flu. This match was based on the assumption that urban residents tend to visit medical institutions near to their dwelling places to get timely treatment for a light ailment, such as cold and flu.
The database of residential health used in this study is different from self-reported data used by the previous studies that investigate the association between exposures to UGSs and health status of interviewees (e.g., [25,27,32,33,58]) We do not know their real contact time to green spaces and individual behavior of urban residents based on the database of the National Health Insurance. However, this approach is better than that using self-reported data in the case of dealing with the relationship between UGSs and health status for whole residents in a big metropolitan area.

Conclusions
Many previous studies focusing on the health benefits of urban green spaces relied on self-reported health and socio-economic data from a relatively small number of respondents. Using the database of the National Health Insurance Research of Taiwan and applying the statistical method of canonical correlation analysis distinguish this study from previous ones. This study shows the detailed relationships between spatial characteristics of UGSs and three diseases that humans suffer from based on the data of morbidities of three diseases. The analytical results revealed a significant positive correlation between the spatial characteristics of UGSs and human health. The results also showed that the districts with lower availability, quality, accessibility to green spaces had a high level of morbidity of cardiovascular diseases, mental disorders, and respiratory diseases. Conversely, districts with high aggregation, availability, quality, and accessibility to green spaces apparently had a low level of morbidity.
Since previous studies of UGSs and human health have focused on the risks of a specific disease, this study investigated the relationships between the spatial characteristics of UGSs and levels of general human health. Instead of relying on the morbidity of a specific disease, this study simultaneously investigated the relationships between the morbidities of three diseases and the spatial characteristics of UGSs. The analytical results reveal that human health generally has a strong correlation with the spatial characteristics of urban green spaces. The results also reveal that better biophysical property of vegetation canopy in urban areas is very important to improve human health. We emphasize not only the importance of the availability of UGSs for human health but the contributions of quality, accessibility, and spatial patterns of UGSs to human health. This study revealed that Network Distance and NDVI of UGSs are the main contributors to human health. Respiratory disease is shown by our results as the most influenced disease by spatial characteristics of urban green spaces. In addition to increasing the number of UGSs, improving the quality, accessibility, and spatial patterns of UGSs can reduce the morbidities of diseases and promote better urban human health. The findings of this study can be useful for maintaining precious green spaces in cities and for urban planning to improve green infrastructure in urban areas.