How Do Urban Parks Provide Bird Habitats and Birdwatching Service? Evidence from Beijing, China

: Parks are an important green infrastructure. Besides other beneﬁts for human and animals, parks provide important bird habitats and accommodate most human-bird interactions in cities. Understanding the complex dynamics among park characteristics, bird habitats and park attractiveness to birdwatchers will inform park designers and managers. However, previous studies often examined factors inﬂuencing bird habitats and birdwatching activities separately. To ﬁll this gap, we aim to study the whole picture of “parks, birds and birdwatchers” in Beijing, China for its spatial patterns and possible factors which inﬂuence bird habitat areas and birdwatching services. We conducted a three-month bird census in at 159 sites and mapped bird habitat areas in parks of Beijing through the maximum entropy method based on results of the bird survey as well as high-resolution remote sensing data. We derived the number of birdwatching records to describe birdwatching activities from the China Birdwatching Record Center website. We used correlation analysis, regression and analysis of variance to investigate factors that may inﬂuence areas of bird habitats and the number of birdwatching records for each park. Our results showed that among the 102 parks, 61 provide habitats to breeding birds with an average of 17 ha, and 26 parks generated a total of 330 birdwatching records. Park size, age, proportion of pavement, landscape connectedness, pavement largest patch index and woodland patch density explained 95% of the variation in habitat areas altogether. Bird habitat area alone explained 65% of the variation in the number of birdwatching records. Furthermore, parks with birdwatching records are signiﬁcantly larger, older, closer to the city center and more accessible than those have no reported birdwatching. These ﬁndings have important implications for park management. While park size or age cannot be easily changed, modifying landscape patterns can increase bird habitats in parks, and improving accessibility may attract more birdwatchers to parks that already have considerable bird habitats.


Introduction
Cities are densely populated and highly modified ecosystems [1], where plants and animals provide important ecosystem services to urban residents [1]. Birds, in particular, play an essential role in the wellbeing of urban dwellers. Researchers found that the presence of birds significantly contributes to people's mental health and self-rated wellbeing [2,3]. For example, birds' colorful look was found to deliver aesthetic enjoyment and make people feel pleasant [4,5]. Moreover, birdsong diversity increased people's self-reported wellbeing and appreciation of the landscape [6,7]. Beijing has a total of 21 orders, 74 families, 423 species and 173 breeding birds [67]. Among the 423 bird species in Beijing, 30 species were listed as threatened species by the International Union for Conservation of Nature [68]. The first organized birdwatching activity in China took place in Beijing in 1996 [69,70]. Nowadays, birdwatching is a popular activity. Every week, there is at least one public birdwatching event taking place in the parks of Beijing [71,72]. With a variety of parks, abundant bird species and vigorous birdwatching activities, Beijing provides a suitable study area to explore our research questions.

Bird Census
We conducted our bird census survey in seven parks. We selected the seven parks as samples out of the 102 parks, considering their size, location and age ( Figure 1, Table S1). We conducted a three-month bird survey during the breeding season, March-July 2017. For each park, we designed a survey route along the road, which covered most park areas. Along the route, we selected observation sites every 150 m for parks smaller than 150 ha and every 200 m for those larger than 200 ha ( Figure 2). As a result, we had 159 sample points in total. We used the sample point count method in the bird census [73]. At each sample point, we stayed for five minutes and counted the species and Beijing has a total of 21 orders, 74 families, 423 species and 173 breeding birds [67]. Among the 423 bird species in Beijing, 30 species were listed as threatened species by the International Union for Conservation of Nature [68]. The first organized birdwatching activity in China took place in Beijing in 1996 [69,70]. Nowadays, birdwatching is a popular activity. Every week, there is at least one public birdwatching event taking place in the parks of Beijing [71,72]. With a variety of parks, abundant bird species and vigorous birdwatching activities, Beijing provides a suitable study area to explore our research questions.

Bird Census
We conducted our bird census survey in seven parks. We selected the seven parks as samples out of the 102 parks, considering their size, location and age ( Figure 1, Table S1). We conducted a three-month bird survey during the breeding season, March-July 2017. For each park, we designed a Remote Sens. 2020, 12, 3166 4 of 15 survey route along the road, which covered most park areas. Along the route, we selected observation sites every 150 m for parks smaller than 150 ha and every 200 m for those larger than 200 ha ( Figure 2). As a result, we had 159 sample points in total. We used the sample point count method in the bird census [73]. At each sample point, we stayed for five minutes and counted the species and the number of individual birds sited within a 25 m radius. We conducted each survey within the first four hours after sunrise.

Possible Influencing Factors: Park Characteristics and Land Cover
We considered factors which may influence bird habitats and birdwatching from several aspects: accessibility and urbanization degree of a park and its surrounding area, park size and age and the landscape composition and configuration in a park (Table S2).
We used high spatial resolution land cover data extracted from the Pleiades satellite imagery to quantify the composition and configuration of park landscapes. These landscape composition and configuration variables were further used to investigate their potential association with bird habitat area and number of birdwatching records in parks. The Pleiades image was acquired on 19 September 2015. It has one 0.5 m spatial resolution panchromatic band (470-830 nm) and four multispectral bands with 2 m spatial resolution (blue: 430-550 nm; green: 500-620 nm; red: 590-710 nm; nearinfrared: 740-940 nm). The image was preprocessed using ERDAS TM . The land cover classification included seven types: woodland (i.e., woody vegetation), grassland (i.e., herbaceous vegetation), building, road, pavement, waterbody and bare land [61]. An object-based approach was used for the land cover classification. We first segmented the image into objects, using the fractal net evolution approach that is embedded in the software eCognition TM . The segmentation algorithm uses a bottomup region merging method, starting with each pixel in the image as an individual segment which is then merged into larger ones. After segmentation, the image objects were then classified using a rulebased classifier built on a class hierarchy and associated features and rules for classifications [61]. Extensive manual editing was conducted to further refine the classification. The overall accuracy of the classification was 84.56%, with a Kappa value of 0.82. The producer's and user's accuracies ranged from 71.43% to 100% and from 81.08% to 88.24%, respectively.
In addition, we also used medium resolution (30m) Landsat ETM+ images to calculate the normalized difference vegetation index (NDVI) and a Digital Elevation Model (DEM 30) obtained from the U.S. Geological Survey (USGS) to calculate slope through 3D analysis (Table 1). Land cover types, NDVI and slope were used as input of MaxEnt model, which will be explained in Section 2.3. The Landsat ETM+ image was retrieved in summer 2016. The DEM 30 was retrieved in 2011. Park boundaries and entrances were manually mapped according to an online mapping platform (Amap.com) and validated through a ground-truth survey from October 14th to November 20th 2015

Possible Influencing Factors: Park Characteristics and Land Cover
We considered factors which may influence bird habitats and birdwatching from several aspects: accessibility and urbanization degree of a park and its surrounding area, park size and age and the landscape composition and configuration in a park (Table S2).
We used high spatial resolution land cover data extracted from the Pleiades satellite imagery to quantify the composition and configuration of park landscapes. These landscape composition and configuration variables were further used to investigate their potential association with bird habitat area and number of birdwatching records in parks. The Pleiades image was acquired on 19 September 2015. It has one 0.5 m spatial resolution panchromatic band (470-830 nm) and four multispectral bands with 2 m spatial resolution (blue: 430-550 nm; green: 500-620 nm; red: 590-710 nm; near-infrared: 740-940 nm). The image was preprocessed using ERDAS TM . The land cover classification included seven types: woodland (i.e., woody vegetation), grassland (i.e., herbaceous vegetation), building, road, pavement, waterbody and bare land [61]. An object-based approach was used for the land cover classification. We first segmented the image into objects, using the fractal net evolution approach that is embedded in the software eCognition TM . The segmentation algorithm uses a bottom-up region merging method, starting with each pixel in the image as an individual segment which is then merged into larger ones. After segmentation, the image objects were then classified using a rule-based classifier built on a class hierarchy and associated features and rules for classifications [61]. Extensive manual editing was conducted to further refine the classification. The overall accuracy of the classification was 84.56%, with a Kappa value of 0.82. The producer's and user's accuracies ranged from 71.43% to 100% and from 81.08% to 88.24%, respectively.
In addition, we also used medium resolution (30 m) Landsat ETM+ images to calculate the normalized difference vegetation index (NDVI) and a Digital Elevation Model (DEM 30) obtained from the U.S. Geological Survey (USGS) to calculate slope through 3D analysis (Table 1). Land cover types, NDVI and slope were used as input of MaxEnt model, which will be explained in Section 2.3. The Landsat ETM+ image was retrieved in summer 2016. The DEM 30 was retrieved in 2011. Park boundaries and entrances were manually mapped according to an online mapping platform (Amap.com) and validated through a ground-truth survey from 14 October to 20 November 2015 [74]. Since land cover and vegetation structure in parks did not experience significant change over recent years, we believe that the time gap between different data acquisitions had little significant impact on the results of the study.

Birdwatching Record
We obtained birdwatching records from the website of the China Birdwatching Record Center [75]. The China Birdwatching Record Center is a citizen science website where registered users can upload birdwatching records. A typical record contains the following information: location, date, recorder ID, weather, watching equipment, watching route, bird species and the number of bird individuals. The most active time for the China Birdwatching Record Center website was 2008 to 2014. Then, it experienced a decline with the growth of smartphone-based social media. We derived all the birdwatching records submitted within the 102 parks between 1 January 2009 and 31 December 2013 and used the total number of birdwatching records in each park to represent its service for birdwatching.

Modeling Bird Habitats in Parks
We used the maximum entropy method (MaxEnt) to estimate bird habitats in the 102 parks based on results from the bird census survey and land cover dataset. In this study, we focused on habitats for breeding birds to represent overall bird habitats because breeding birds have a relatively stable nest location, territory and well-documented breeding traits (nest location and nest type, body length and body weight) [26,67,76,77]. The point of bird occurrence is likely to be within their breeding territories, making a strong link with the spatial properties of the habitat.
First, we identified the breeding birds according to the bird checklist in Beijing [67]. Then, we grouped them into different guilds by traits including body length, body weight, nest location and nest type. Next, we ran the MaxEnt model to estimate distribution for each guild. Finally, we merged the distribution maps for all guilds to generate a bird habitat map.
In order to group breeding birds into different guilds, we input body length, body weight, nest location and nest type data to draw a dendrogram in the R TM program [76,77]. We used the Functional Diversity (FD) package to calculate the Gower distance matrix of these traits, which evaluated the overall dissimilarity of species or taxa [78]. We chose the average method of hierarchical clustering to generate a trait dendrogram. The branch lengths in the resulting dendrogram between clusters were calculated by pairwise distances. We decided to use the distance of 0.18 as a threshold to balance the number of guilds and their dissimilarity. As a result, we divided the breeding birds into thirteen functional guilds (Figure 3) and generated a presence-absence dataset for each guild. We have included the command syntax for R in the Supplementary Materials. Detailed steps on how to generate a dendrogram can be found in this book [79]. Finally, we input the bird presence-absence dataset along with land cover type, normalized difference vegetation index, slope, distance to city center and distance to the nearest water body into the MaxEnt model (Table 1). It worth noting that the MaxEnt model requires a minimum sample size of 14-that is, birds of a guild observed at 14 sites or more. Therefore, we excluded the water surface open nest guild (observed at five sites) and the reed open nest guild (observed at one site) from further analysis.
Remote Sens. 2020, 11, x FOR PEER REVIEW 6 of 16  As a result, the MaxEnt model generated a probability map of distribution probability for each guild at the resolution of 30 m × 30 m. We used 25% presence data for model testing. For each guild, we ran the Digital Elevation Model14 times to generate the distribution probability map. We assessed the model accuracy by the threshold-independent evaluation statistic method. The threshold is the As a result, the MaxEnt model generated a probability map of distribution probability for each guild at the resolution of 30 m × 30 m. We used 25% presence data for model testing. For each guild, we ran the Digital Elevation Model14 times to generate the distribution probability map. We assessed the model accuracy by the threshold-independent evaluation statistic method. The threshold is the area under the curve (AUC) of the receiver's operating characteristic [80]. After this, we defined habitat as a distribution probability value higher than 0.60. Then, we merged the suitable habitats for all the guilds (Table S3) to represent bird habitat and calculated the total area.

Analyses
We first conducted a correlation analysis between the number of birdwatching records and habitat areas generated by MaxEnt modeling. We excluded parks with no habitat or birdwatching record in this analysis.
Then, we used stepwise regression analysis to investigate how possible influencing factors may explain variations in habitat area and birdwatching records. We conducted two stepwise regression analyses with bird habitat area and the number of birdwatching records as dependent variables, respectively. Independent variables for both regressions included three aspects: park size and age, accessibility and urbanization degree of a park and its surrounding area and the landscape composition and configuration in a park (Table S2). In addition, we added bird habitat area as an independent variable in the regression which had the number of birdwatching records as a dependent variable. We tested collinearity between the variables by correlation analysis and variance inflation coefficient (VIF). The results were included in Tables S4 and S5, respectively. Stepwise regression was able to identify independent variables that were highly correlated and the final model only contained independent variables with a relatively low degree of autocorrelation.
Finally, we conducted a one-way analysis of variance (ANOVA) to investigate the differences between parks with and without birdwatching records (Table S2). We separated parks with and without birdwatching records into two groups and compared them to determine whether any of the variables in terms of park size and age, accessibility, urbanization degree and landscape composition and configuration were significantly different across the two groups. All the analyses were conducted in SPSS TM . Additionally, we used ArcGIS TM 10.0 to visualize the spatial distributions of bird habitat area and the number of birdwatching records.

Bird Census and Modeling
We found 79 species from the bird survey, among which 53 were breeding species. We identified 13 functional guilds for the 53 breeding species based on their breeding traits (Table S3). As the MaxEnt model requires a minimum sample size of 14, 11 out of the 13 functional guilds were entered into the model. The model passed the AUC accuracy test with the values ranging from 0.71 to 0.85 (Table S6). More details on the accuracy test are included in the Supplementary Materials.

Spatial Patterns and Relationships of Park's Bird Habitats and Birdwatching Records
Our results showed that 61 parks had habitats at least for one breeding guild. The area of habitat in each park ranged from 0.09 ha to 169 ha, with an average of 17 ha. The habitat proportion ranged from 0.45% to 71%, with an average of 19%. Most parks with large bird habitat areas were located in the central and northern regions of the city (Figure 4a)

Park Characteristics Explaining Variations in Bird Habitat Area and Birdwatching Activities
The stepwise regression results showed that park size, age, proportion of pavement, landscape connectedness, pavements largest patch index and woodland patch density ( Table 2) explained 95% of the variation in bird habitat area. Parks that are larger and older and have less pavement, better connected landscapes and less woodland patches tend to have larger areas of suitable bird habitats.
In terms of the number of birdwatching records, our results showed that park size and the woodland landscape shape index together explained 80% of its variation ( Table 2). Parks that are larger and have woodland in relatively regular shapes tend to have more birdwatching records. However, after we added bird habitat area as an independent variable, park size and woodland landscape shape index were no longer significant. Results indicated that bird habitat area explained 65% of the variation in birdwatching records. Table 2. Results of stepwise regression models of factors that may impact bird habitat area and number of birdwatching records in parks.

Park Characteristics Explaining Variations in Bird Habitat Area and Birdwatching Activities
The stepwise regression results showed that park size, age, proportion of pavement, landscape connectedness, pavements largest patch index and woodland patch density ( Table 2) explained 95% of the variation in bird habitat area. Parks that are larger and older and have less pavement, better connected landscapes and less woodland patches tend to have larger areas of suitable bird habitats. In terms of the number of birdwatching records, our results showed that park size and the woodland landscape shape index together explained 80% of its variation ( Table 2). Parks that are larger and have woodland in relatively regular shapes tend to have more birdwatching records. However, after we added bird habitat area as an independent variable, park size and woodland landscape shape index were no longer significant. Results indicated that bird habitat area explained 65% of the variation in birdwatching records.
Results from the one-way ANOVA analysis revealed significant differences in park size between parks with and without birdwatching records (Table 3, Figure 5). Firstly, parks with birdwatching records were significantly larger than their counterparts. After the effect of park size was considered, results showed that parks with birdwatching records were older, closer to city center and more accessible than those that had no reported birdwatching.

Spatial Patterns of Parks with Bird Habitats and Birdwatching Record
Our results indicated that parks with relatively more habitats and/or birdwatching activities clustered along the "central axis", which started from the Forbidden City (the center of Beijing in both ancient and modern times), stretched all the way to the north and ended in the Olympic Forestry

Spatial Patterns of Parks with Bird Habitats and Birdwatching Record
Our results indicated that parks with relatively more habitats and/or birdwatching activities clustered along the "central axis", which started from the Forbidden City (the center of Beijing in both ancient and modern times), stretched all the way to the north and ended in the Olympic Forestry Park (designed with a dragon-shaped waterbody to echo the Forbidden City). Parks along the central axis often have historical and cultural legacies that can be traced back to dynasties. For example, the Beihai Park was built in 1163 and has served as an imperial home and garden since 1179. The Beijing Zoo was built in the Qing Dynasty upon royal gardens. Researchers found that high density of old trees contributes to bird richness [24], which could explain the connections between ancient parks and suitable bird habitats. In contrast, parks between the third and the fifth ring roads were built more recently, provided less bird habitats and attracted less birdwatchers.
Results showed that numbers of birdwatching records in parks located in the northwest area were high. The northwest area in Beijing hosts the IT companies (known as the Chinese "Silicon Valley"), quite a few well-known universities and residents with good education (Figure 4). The spatial overlap between highly educated residents and parks with active birdwatching activities may not be a coincidence. Birdwatchers usually have a relatively high education level [70]. Moreover, living close to parks tended to increase birdwatching and interacting with birds [15,17,32]. The spatial overlaps observed in this study confirmed these previous findings regarding the education level of birdwatchers and "nearby birdwatching" in urban areas.

Relationship between Parks with Bird Habitats and Birdwatching Record
We found the area of bird habitat was the most important independent variable explaining the variation in the number of birdwatching records. Habitat area is usually considered to be associated with bird species richness [30], so this finding is consistent with previous studies which found that bird species richness is a major indicator of locations where birdwatching activities may take place [17,35,36]. With that being said, bird habitat was clearly not a sufficient condition to attract birdwatchers. Among the 66 parks with bird habitats, only 22 of them had birdwatching records. Accessibility played an important role in attracting birdwatchers. Our analysis showed that parks with birdwatching records were closer to the downtown area with more roads, bus stops and metro stations nearby, compared with parks with no birdwatching records.
Interestingly, bird habitat was not a necessary condition for birdwatching either. Among the 26 parks with birdwatching records, four of them were estimated to have few bird habitats according to our results. We did not conduct the one-way ANOVA analysis to examine differences in park characteristics' due to the small sample size. These four parks were all relatively small (0.8 to 8 ha) compared to the average park size of 41 ha. We might underestimate bird habitat in the four small parks due to the limitation of MaxEnt modeling, which will be discussed later in the limitation section. In addition, there are other parks within a 1-mile radius of these four parks, which may provide habitats for birds. Lastly, birdwatching records indicated that birds which appeared in the four parks were Eurasian Tree Sparrow (Passer montanus), Common Magpie (Pica pica), Azure-winged Magpie (Cyanopica cyanus), Large-billed Crow (Corvus macrorhynchos) and Carrion Crow (Corvus corone). These species are resident birds well adapted to urban environments, which contain a large proportion of impervious surfaces (such as buildings, pavements and roads). Future on-site study is needed to investigate parks with few habitats hosting birdwatchers by using high-resolution spatial data.
It is also worth noting that environmental non-government organizations (NGOs) played an essential role in introducing and promoting birdwatching activities [18,70]. The three parks with the most birdwatching records all host birdwatching activities organized by NGOs on a regular basis.

Impacts of Park Characteristics on Bird Habitats and Birdwatching
Results showed that park size was an important factor. Larger parks tend to have more bird habitats and birdwatching records. Park size is commonly used as an indicator to describe the quality of a park in studies evaluating their attractiveness to people as well as their ecosystem services in general [81,82]. Our findings provide evidence that parks with larger sizes contribute to both supporting bird species as well as hosting birdwatching activities. Despite these benefits of large parks, it is not always feasible to enlarge a park. Our results indicate that certain landscape patterns, such as less pavements, may increase bird habitats and birdwatching records, which could be achieved through proper design. Previous studies also found that larger parks usually have high bird richness [9,16,21] and also found that less pavements can promote bird richness [29].
Furthermore, our results also showed that accessibility determined whether a park would be attractive to birdwatchers, which is consistent with previous studies conducted in protected areas and national parks [15,17]. We found that small parks with limited bird habitats but good accessibility still hosted some birdwatching activities. For example, Liuyin Park is 15 ha in size, with 3 ha bird habitat (compared with the mean values of 41 and 10, respectively). With a score of 51 in accessibility (number of metro and bus stops, compared with the mean value of 31), it contributed six out of 330 birdwatching records. So, increasing the number of bus stops and subway stations can help to attract more birdwatchers.

Limitations
This study has several limitations. Firstly, we only considered birdwatchers for the benefits brought by birds. Park visitors who did not intend to watch birds or submit birdwatching records may also enjoy the presence of birds in parks, but they were not considered in this approach. Furthermore, using the number of birdwatching records may oversample areas with rare or exotic birds. For example, the Bluethroat (Luscinia svecica) appeared in the Olympic Forest Park and attracted three records within four days between May 16 and 19, 2010. Future studies should consider the possibilities of eliminating such impact, such as surveying visitors to estimate how they enjoy birds in parks.
Secondly, the MaxEnt model worked on a resolution of 30m * 30m and was not able to capture habitats smaller than this scale. This may result in an underestimation of the total bird habitat area in small parks. While it did generate habitat estimates successfully in the 11 out of 33 parks smaller than 8 ha (Table S7), four parks smaller than 8 ha that reported birdwatching records were estimated to have few habitats by the model. Future studies should investigate what causes the inconsistency between suitable bird habitats and birdwatching records in the parks as well as apply high-resolution data and combine them with small-scale surveys to make the habitat estimation of small-scale parks more accurate.
Lastly, we employed land cover based on high-resolution remote sensing data as input for the MaxEnt model to estimate bird habitats in this study. Therefore, we carry the limitations of remote sensing images, which hardly detect the vertical vegetation structure, such as conditions of shrubs under arbor cover and length or quality of grass. Such 3D vegetation structures that are captured by remote sensing data may impact bird habitats [9,16,21,83].

Conclusions
We investigated the factors that impact bird habitats and birdwatching activities in urban parks of Beijing. Our results showed that parks with large areas and a long history tend to provide more bird habitats and birdwatching services. Besides size and age, park landscape pattern also has an impact on bird habitat area. Parks with less pavement cover, less clustered pavements and better-connected landscapes tend to have more bird suitable habitats. Birdwatching activities were highly associated with suitable bird habitats, as the area of habitat can explain 65% of the variation in birdwatching records. In addition, park location and accessibility also contributed to the popularity of a park as a birdwatching destination. Providing habitats in cities and connecting urban residents with nature are the fundamental goals of urban parks. Here, we focus on birds and birdwatching to examine what factors may help or impede parks in fulfilling these objectives. Our findings will inform park planners and managers about strategies and measures that can help them better achieve one or both goals.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2072-4292/12/19/3166/s1, The command syntax of drawing a trait dendrogram in R TM . Table S1. Location, size, age and waterbody of sample parks, Table S2. Independent variables for statistical analysis, Table S3. Name list of breeding guilds, Table S4. Correlation matrix of landscape index, Table S5. The variance inflation coefficient (VIF) of index in regression, Table S6. The area under the curve (AUC) of each guild habitat assessment in MaxEnt model, Table S7. Parks with bird habitats of less than 1.6 ha.