Remotely Sensed Spectral Indices as Proxies of the Structure of Urban Bird Communities

Abstract

Abundant and diverse urban bird communities promote ecosystem and human health in cities. However, the estimation of bird community structure requires large amounts of resources. On the other hand, calculating remotely sensed spectral indices is cheap and easy. Such indices are directly related to vegetation cover, built-up cover, and temperature, factors that also affect the presence and abundance of bird species in urban areas. Therefore, spectral indices can be used as proxies of the structure of urban bird communities. We estimated the abundance, taxonomic, functional, and phylogenetic diversity of the bird community at each of 18 50 m radius survey stations in the urban core area of Kavala, Greece. We also calculated eight spectral indices (means and standard deviations, SDs) around survey stations at 50 m, 200 m, and 500 m spatial scales. The land surface temperature SD (LST) was the most important proxy, positively related to bird abundance at the 50 m and 200 m spatial scales. At the same time, the mean green normalized difference vegetation index (GNDVI) was the most important proxy, negatively related to abundance at the 500 m spatial scale. Means and SDs of vegetation indices, such as the normalized difference vegetation index (NDVI), enhanced vegetation index (EVI2), soil-adjusted vegetation index (SAVI), and atmospherically resistant vegetation index (ARVI), were the most important proxies, positively related to taxonomic and functional diversity at all the spatial scales. The mean and SDs of LST, normalized difference moisture index (NDMI), and normalized difference built-up index (NDBI) variously affected taxonomic and functional diversity. The mean and SDs of LST were the best proxies of phylogenetic diversity at the 50 m and 500 m spatial scales, while the SDs of NDBI and NDMI were the best proxies at the 200 m spatial scale. The results suggest that several spectral indices can be used as reliable proxies of various facets of urban bird diversity. Using such proxies is an easy and efficient way of informing successful urban planning and management.

1. Introduction

Urban areas have occupied formerly natural areas, resulting in the loss of local species populations [1]. On the other hand, urban areas host important diversity, especially of birds, including endemic and threatened species [1]. In urban areas, the replacement of native, urban-avoider bird species by sometimes non-native, dweller species; the exclusion of migrants; and the constant provision of resources, such as food and habitat, often lead to spatial and temporal homogenization in bird composition [2]. Also, high bird diversity and richness in urban green spaces have been connected with increased happiness and decreased stress and anxiety [3]. Therefore, bird diversity is important for both ecosystem stability and human health. The study of bird community structure is necessary to provide urban planners and wildlife managers with critical information for designing, creating, and maintaining bird-diverse urban green spaces [4,5,6].

Diversity has several facets, including taxonomic, functional, and phylogenetic diversity [7]. Taxonomic diversity refers to the abundance and richness of different taxonomic groups and the distribution and relationships among taxa, such as species, genera, families, and higher-level classifications within an area [8,9,10]. Higher taxonomic diversity implies richer, resilient ecosystems and has been widely used for protecting threatened species and areas with many species [4,5,6]. Functional diversity refers to the variety of biological traits and functions of species in a community [11,12,13]. It focuses on species’ roles in ecosystems and ecosystem processes and not on the number of species or their genetic relationships. Phylogenetic diversity measures how closely or distantly related species are on the tree of life, thus providing an understanding of the evolutionary relationships among these species [14,15,16,17]. This biotic homogenization affects all aspects of diversity. Urbanization has caused an overall decrease in taxonomic, functional, and phylogenetic diversity [1,18,19,20] in bird communities. Therefore, all facets of bird diversity should be examined to identify species and areas that face threats or hold important taxon richness, functional traits, and phylogenetic history for establishing conservation priorities.

The amount of vegetation cover and built-up areas are among the most important factors affecting bird abundance, richness, and diversity in urban green spaces [1,4,5,6,21]. However, the monitoring of bird communities and the measurement of environmental variables require time and trained personnel, hence large and constant funding, which is not always available. In contrast, acquiring satellite images is cheap, and calculating remotely sensed spectral indices is easy. Such indices have proven important for predicting various facets of the structure of urban bird communities [22,23,24]. The normalized difference vegetation index (NDVI) is the most used spectral index to evaluate vegetation attributes, such as canopy phenology, leaf area, and primary production [25,26]. Many studies have found that NDVI is positively associated with urban bird communities’ abundance, richness, and diversity [22,24,26,27,28,29,30]. The enhanced vegetation index (EVI) is another spectral index used to evaluate vegetation. NDVI and EVI are very similar, but the latter is more effective in low and high vegetation cover and less effective in hilly topography [31,32]. EVI has been used much less often than NDVI as proxies of bird community structure, both in natural areas [33,34,35,36] and urban areas [23]. Variants to NDVI and EVI include the green normalized difference vegetation index (GNDVI) [37], soil-adjusted vegetation index (SAVI) [38,39], atmospherically resistant vegetation index (ARVI) [40,41], and normalized difference moisture index (NDMI) [42]. Among these indices, only SAVI has been used for inferring bird richness and distribution during migration [43,44]. The normalized difference built-up index (NDBI) is an index used to identify urban and built-up areas from satellite imagery [45]. Gray infrastructure cover has been reported as a good predictor of bird diversity in urban areas [4,6,46,47,48]. However, NDBI has not been used as a proxy for bird community structure. Land surface temperature (LST) can change significantly inside a relatively small heterogeneous urban area and is measured by satellites that use thermal infrared sensors [49]. Differences in LST are closely connected to differences in vegetation density. They can be especially useful in the heterogeneous urban landscape with the great difference in temperature between green and gray infrastructure [50]. However, this index has not been used very often to infer bird diversity [33,43,51,52]. Means and standard deviations (SDs) of spectral indices are mainly used as proxies of bird community structure. SD values are surrogates of habitat heterogeneity, and previous studies have found a positive association between vegetation heterogeneity in a given space and bird diversity [23,53].

The aim of this study was to examine if several spectral indices measuring vegetation cover, built-up cover, and temperature can be used as proxies of the taxonomic, functional, and phylogenetic diversity of the urban bird community in Kavala, Greece. Inferring various facets of bird community structure in urban areas through cheap and easily retrieved and measured methods is critical for their effective conservation management.

2. Materials and Methods

2.1. Study Area

The study was carried out in Kavala, Northern Greece (Figure 1). The urban area of the municipality of Kavala occupies about 8.0 km² and has about 56,300 inhabitants [54]. It is delineated by a Turkish pine (Pinus brutia) forest to the north and the Aegean Sea to the south. Most green spaces of Kavala are smaller than 3 ha, except for the 17 ha Panagiouda pine woodland, including square gardens, playgrounds, hedgerows, and median road strips [6].

2.2. Fieldwork

Birds were counted in 18 green spaces, including woodlands (5; survey stations 2, 5, 7, 8, and 16), square gardens (5; 1, 9, 11, 12, and 15), playgrounds (3; 3, 6, and 10), and median road strips (5; 4, 13, 14, 17, and 18) [6]. One survey station was located at the center of each green space and at least 250 m away from other stations to avoid counting a bird twice. Birds were counted within a 50 m radius from the survey station center [55]. We visited the survey stations in 2016—once in April, once in May, and once in June to cover the whole breeding season and avoid missing early or late nesters. Counts were carried out within four hours from dawn when birds were more active. When the observer (E.V.) arrived at the center of the station, he remained silent for 5 min and then counted birds for another 5 min [55]. Only one researcher carried out counts to avoid observer effects. Every bird seen within the 50 m radius was recorded. Birds flying overhead were counted only if connected to the green space (e.g., feeding, courting).

2.3. Bird Community Indices

Bird abundance, taxonomic, functional, and phylogenetic diversity indices were used to assess whether spectral indices can be used as proxies in the urban green spaces of Kavala. Abundance was estimated as the maximum number of individuals counted in each survey station during the three visits (Table S1).

Rarefied indices were used to measure taxonomic diversity because they account for the non-detectability of rare species: the Chao1 richness estimator, the probability of interspecific encounter index (PIE), and Shannon entropy (Shannon–Wiener index) (see

Table 1. Taxonomic, functional, and phylogenetic diversity indices used in this study.

Index	Code	Definition
Taxonomic diversity
Species richness	Chao 1 richness	Chao1 estimator, the lower bound of undetected species richness in terms of the numbers of singletons and doubletons [8]
Probability of Interspecific Encounter	PIE	The probability that two randomly sampled individuals from the assemblage represent two different species; a diversity metric not sensitive to rare species; high PIE values indicate high species evenness [9]
Shannon entropy (Shannon–Wiener index)	Chao Shannon	A diversity index that weighs species exactly by their frequencies, without favoring rare or common species [10]
Functional diversity
Functional richness	FRic	It is the amount of functional space occupied by a community; high values suggest that most of the existing niches are used by the species in that community [11]
Functional evenness	FEve	It represents the evenness of distribution of abundance in a functional trait space [11,12]
Functional divergence	FDiv	It is the distance of species abundances to the center of the functional space in a community; it represents the overall functional diversity in that community [11,12]
Functional dispersion	FDis	The mean distance in multidimensional trait space of individual species to the centroid of all species; can account for species abundances by shifting the position of the centroid toward the more abundant species and weighting distances of individual species by their relative abundances [13]
Rao’s quadratic entropy	Rao’s Q	Includes both the relative abundances of species and a measure of the pairwise functional differences between species in a community; conceptually similar to FDis [15]
Phylogenetic diversity
Community evolutionary distinctiveness	CED	The mean evolutionary distinctiveness of the species in a community; a species’ evolutionary distinctiveness score is a measure of the species’ uniqueness [16,17]
Maximum community evolutionary distinctiveness	maxCED	The maximum value of a species’ evolutionary distinctiveness score recorded in a community [16,17]

for definitions and Table S1 for data). The Chao1 estimator was calculated with the function ChaoSpecies and Shannon entropy with the function Diversity of the SpadeR R package [56], while the PIE index was calculated with the function calc_PIE of the mobr package in R [57].

Functional diversity was measured with the functional richness (FRic), functional evenness (FEve), functional divergence (FDiv), functional dispersion (FDis), and Rao’s quadratic entropy (Rao’s Q) indices (see Table 1 for definitions and Table S1 for data). Body mass, food type (seeds, fruits, vegetation, invertebrates, fish, mammals, birds, carrion), feeding substrate (water surface, water, ground, canopy, shrub, vegetation, air) (16 traits taken from the avian niche database of Pearman et al. [58]), and resident or migrant status (2 traits; see Table S2), totaling 18 traits, were used to measure functional diversity indices. All functional diversity indices were calculated with the function dbFD of the FD package in R 4.2.2 [13,59].

The community evolutionary distinctiveness (CED) and the maximum community evolutionary distinctiveness (maxCED) were used as measures of phylogenetic diversity (see Table 1 for definitions and Table S1 for data). The median ED score for each species was taken from the database in Jetz et al. [14], and the CED of each survey station was calculated as the mean ED for all species observed in the survey station [18]. The maxCED was the maximum value of ED in each survey station [20].

2.4. Remotely Sensed Spectral Indices

Eight remotely sensed spectral indices (mean values and standard deviations, SDs) were assessed as proxies of the structure of the urban bird community of Kavala: NDVI, GNDVI, EVI2, SAVI, ARVI, NDMI, NDBI, and LST (see

Table 2. Remotely sensed spectral indices (means and standard deviations, SDs) used to assess their relationship with the abundance and diversity of the bird community of Kavala.

Index	Code	Definition
Normalized difference vegetation index	NDVI Mean	It uses an algorithm that extracts information from two channels of a satellite image, red and near-infrared (NIR); its values range from −1 to +1 depending on the relative reflectance of geographic features in the two spectral bands; vegetated areas tend to give high NDVI values due to the high reflectance of green vegetation in NIR and low reflectance in the red band; rocks and impervious surfaces have similar reflectance in both bands and give values close to zero; open water gives negative values [25].
	NDVI SD
Green normalized difference vegetation index	GNDVI Mean	An index of the plant’s “greenness” or photosynthetic movement; uses the green instead of the red band used in NDVI; ranges between −1 and +1 [37].
	GNDVI SD
Enhanced vegetation index 2	EVI2 Mean	Quantifies vegetation greenness; similar to NDVI; calculated using the near-infrared (NIR) and red spectral bands; corrects for some atmospheric conditions and canopy background noise; more sensitive in areas with dense vegetation; reduces the limitations imposed by the blue band in the EVI; ranges between −2.5 and +2.5 [31,32].
	EVI2 SD
Soil-adjusted vegetation index	SAVI Mean	Developed to reduce the influence of soil on canopy spectra by incorporating a soil adjustment factor (L) into the denominator of the NDVI equation; ranges between −1.5 and +1.5 [38,39].
	SAVI SD
Atmospherically resistant vegetation index	ARVI Mean	An enhanced vegetation index calculated using the near-infrared (NIR), red, and blue spectral bands; resistant to atmospheric effects due to a self-correction process for the atmospheric effect on the red band; ranges between −1 and +1 [40,41].
	ARVI SD
Normalized difference moisture index	NDMI Mean	Assesses vegetation water content in plants; calculated using the near-infrared (NIR) and short-wave infrared (SWIR) spectral bands; ranges between −1 and +1; the lowest values indicate low vegetation water content and the highest values correspond to high water content [42].
	NDMI SD
Normalized difference built-up index	NDBI Mean	Used to identify urban and built-up areas from satellite imagery; leverages the near-infrared (NIR) and short-wave infrared (SWIR) spectral bands; ranges between −1 and +1 [45].
	NDBI SD
Land surface temperature	LST Mean	The temperature of the land surface in degrees Celsius (°C); can change significantly inside a relatively small heterogeneous urban area; satellites use thermal infrared (TIR) sensors to measure the heat emitted from the Earth’s surface [49].
	LST SD

for definitions and Table S3 for data).

The spectral indices, except LST, were calculated on a Sentinel-2 image sensed on 7 October 2016. We did not use multiple images because of the very small changes in the urban landscape during the short period of the study. The Sentinel-2 mission consists of two polar-orbiting satellites operated by the European Space Agency [60]. They carry a multispectral instrument sensor that can produce data at spatial resolutions from 10 m to 60 m. The image was retrieved from the Copernicus Scientific Hub and processed at Level 2. The spectral bands in the visual, near-infrared (NIR), and short-wave infrared (SWIR) parts of the electromagnetic spectrum, at spatial resolutions of 10 and 20 m, were employed to calculate the spectral indices. The land surface temperature was retrieved using Band 10 of a Landsat 8 image Level 1C image, sensed on 25 June 2016. The original B10 Digital Numbers were converted to Top-of-Atmosphere (ToA) radiance using the rescaling coefficients provided in the Landsat 8 metadata file. The ToA radiance was then converted to Brightness Temperature (BT) in degrees Celsius using the band-specific thermal conversion constants, provided in the Landsat 8 metadata file. The BT was then corrected using the surface emissivity estimated based on the NDVI value, calculated using Bands 5 and 4 of the Landsat 8 image [61].

The mean and SD values of the 8 spectral indices were calculated in three spatial scales at each survey station. These buffer zones were circles with radius 50 m, 200 m, and 500 m, centered at each station. The size of the buffer zones was selected based on the bird species’ territory and colony, and the size and configuration of Kavala’s gray and green infrastructure [22]. When a buffer extended into the sea, its surface was truncated for measurement accuracy.

2.5. Data Analysis

Generalized Linear Models (GLMs) were used to examine the relationships between bird abundance; taxonomic, functional, and phylogenetic diversity (response variables); and the 8 spectral indices (predictors). Predictors were modeled separately, and there was a model for each predictor (mean and SD) and each response variable, resulting in 16 models for each spatial scale. This analysis design was selected because we aimed to assess the single capacity of each spectral index as a proxy of each diversity index.

Abundance was a count variable modeled with Poisson distribution and log link function. All other response variables were continuous and modeled with Gaussian distributions and identity link functions. Models were fitted using the function glm of the stats package in R 4.2.2 [62]. Rate ratios, equivalent to odds ratios in logistic regression, were calculated to assess the rate of change in the response variable for a one-unit change in the predictor variable [63,64]. Model selection was performed using the Akaike information criterion (AIC_c) with the function model.sel of the MuMIn package in R [63]. Models differing less than 2 AIC_c units from the best model (ΔAIC_c ≤ 2) were selected as the most parsimonious, while the model with the lowest AIC_c was considered the best model [64]. Plots of the best models were produced with the visreg package in R [65].

3. Results

3.1. Spectral Indices as Proxies of Bird Community Structure at the 50 m Spatial Scale

The most parsimonious models for the 50 m buffer zone are shown in

Table 3. The most parsimonious GLM models (ΔAIC_c ≤ 2) selected to assess the relationship of remotely sensed spectral indices, at a 50 m buffer around survey stations, on the abundance and taxonomic, functional, and phylogenetic diversity of the bird community in Kavala.

	LogLik	AICc	ΔAICc	wi	Cumulative wi	Rate Ratio (95% CI)
Abundance
LST SD	−90.52	185.85	0.00	1.00	1.00	1.80 (1.51–2.14)
Richness (Chao1)
NDVI Mean	−41.49	90.69	0.00	0.29	0.29	31.96 (16.67–47.25)
SAVI Mean	−41.49	90.69	0.00	0.29	0.58	21.31 (11.12–31.50)
EVI2 Mean	−41.74	91.19	0.50	0.22	0.80	16.86 (8.57–25.14)
Evenness (PIE)
NDVI Mean	23.02	−38.33	0.00	0.22	0.22	1.69 (1.10–2.58)
SAVI Mean	23.02	−38.33	0.00	0.22	0.44	1.42 (1.07–1.88)
EVI2 Mean	22.90	−38.09	0.24	0.20	0.64	1.32 (1.05–1.65)
GNDVI Mean	22.82	−37.93	0.40	0.17	0.81	1.74 (1.09–2.79)
Diversity (Shannon entropy)
NDVI Mean	−3.76	15.24	0.00	0.29	0.29	36.74 (5.61–67.79)
SAVI Mean	−3.76	15.24	0.00	0.29	0.58	11.06 (3.16–38.72)
EVI2 Mean	−4.01	15.73	0.49	0.23	0.81	6.66 (2.40–18.45)
ARVI Mean	−4.38	16.48	1.23	0.17	0.98	37.09 (4.85–83.64)
Functional richness (FRic)
NDBI SD	−87.00	181.71	0.00	0.17	0.17	2.78 (1.39–4.45)
ARVI Mean	−87.58	182.88	1.17	0.09	0.26	1.83 (1.47–2.71)
ARVI SD	−87.63	182.98	1.27	0.09	0.35	1.52 (1.25–1.99)
EVI2 SD	−87.77	183.26	1.54	0.08	0.43	1.33 (1.04–2.47)
Functional evenness (FEve)
GNDVI Mean	20.15	−32.59	0.00	0.19	0.19	1.73 (1.03–2.14)
LST SD	19.72	−31.73	0.86	0.12	0.31	0.89 (0.78–0.99)
SAVI Mean	19.66	−31.61	0.99	0.11	0.42	1.34 (1.01–1.84)
NDVI Mean	19.66	−31.61	0.99	0.11	0.53	1.55 (1.02–2.58)
EVI2 Mean	19.59	−31.46	1.13	0.11	0.64	1.26 (0.95–1.65)
Functional diversity (FDiv)
LST SD	26.13	−44.55	0.00	0.18	0.18	0.91 (0.85–0.98)
GNDVI Mean	25.98	−44.25	0.30	0.16	0.34	1.43 (1.06–2.02)
LST Mean	25.62	−43.53	1.02	0.11	0.44	1.01 (1.00–1.03)
Functional dispersion (FDis)
ARVI SD	11.23	−14.75	0.00	0.28	0.28	46.12 (7.5–84.23)
EVI2 SD	11.02	−14.33	0.42	0.23	0.51	5.17 (1.27–21.09)
NDVI SD	10.67	−13.63	1.12	0.16	0.67	26.92 (1.53–72.9)
SAVI SD	10.33	−12.95	1.80	0.11	0.78	8.98 (1.33–40.72)
Rao’s quadratic entropy (Q)
EVI2 SD	−37.29	82.29	0.00	0.20	0.20	16.84 (3.68–29.99)
ARVI SD	−37.44	82.59	0.30	0.17	0.37	35.23 (6.88–63.58)
NDVI SD	−37.46	82.63	0.34	0.18	0.55	33.36 (6.40–60.31)
SAVI SD	−37.46	82.63	0.34	0.18	0.73	22.24 (4.27–40.21)
Community evolutionary distinctiveness (CED)
NDVI Mean	−11.01	29.73	0.00	0.22	0.22	62.05 (3.74–131.99)
SAVI Mean	−11.01	29.73	0.00	0.22	0.44	15.68 (2.41–52.18)
GNDVI Mean	−11.09	29.90	0.17	0.20	0.64	87.62 (3.99–192.68)
EVI2 Mean	−11.32	30.36	0.64	0.16	0.80	8.40 (1.82–38.81)
Maximum community evolutionary distinctiveness (maxCED)
LST SD	−48.56	104.83	0.00	0.22	0.22	8.40 (2.59–14.42)

Each predictor was modeled separately for each response variable. Abundance models were built under Poisson distributions; diversity models were built under Gaussian distributions. LogLik, log-likelihood; AICc, second-order Akaike information criterion; ΔAICc, AICc difference between the model under concern and the best model; wi, Akaike weight; rate ratio, equivalent to odds ratio in logistic regression. See Table 1 for diversity index codes, Table 2 for spectral index codes, and Table S4 for the full list of models.

, while the full list of models is shown in Table S4. Rate ratios ranged from 0.89 to 87.62, showing a considerable effect of spectral indices on bird abundance and diversity (Table 3). LST SD was the best predictor of abundance, FDiv, and maxCED (Table 3). Mean GNDVI and mean LST were also parsimonious predictors of FDiv. Abundance and maxCED increased, while FDiv decreased, with increasing LST SD (Figure 2).

Mean NDVI was the best predictor of Chao 1 richness, PIE evenness, Shannon entropy, and CED (Table 3). Mean SAVI and EVI2 were also parsimonious predictors of these diversity indices, while mean GNDVI was a parsimonious predictor of PIE evenness and CED, and ARVI was a parsimonious predictor of Shannon entropy. Chao 1 richness, PIE evenness, Shannon entropy, and CED increased with increasing mean NDVI (Figure 2).

Mean NDBI was the best predictor of FRic (Table 3). Mean ARVI, ARVI SD, and EVI2 SD were also parsimonious predictors of FRic. FRic increased with increasing mean NDBI values (Figure 2).

Mean GNDVI was the best predictor of FEve (Table 3). LST SD and mean SAVI, NDVI, and EVI2 were also parsimonious predictors of FEve. FEve increased with increasing mean GNBVI values (Figure 2).

ARVI SD was the best predictor of FDis (Table 3). EVI2 SD, NDVI SD, and SAVI SD were also parsimonious predictors of FDis. FDis increased with increasing ARVI SD (Figure 2).

EVI2 SD was the best predictor of Rao’s Q (Table 3). ARVI SD, NDVI SD, and SAVI SD were also parsimonious predictors of Rao’s Q. Rao’s Q increased with increasing EVI2 SD (Figure 2).

3.2. Spectral Indices as Proxies of Bird Community Structure at the 200 m Spatial Scale

The most parsimonious models for the 200 m buffer zone are shown in

Table 4. The most parsimonious GLM models (ΔAIC_c ≤ 2) selected to assess the relationship of remotely sensed spectral indices, at a 200 m buffer around survey stations, on the abundance and taxonomic, functional, and phylogenetic diversity of the bird community in Kavala.

	LogLik	AICc	ΔAICc	wi	Cumulative wi	Rate Ratio (95% CI)
Abundance
LST SD	−91.77	188.33	0.00	1.00	1.00	1.51 (1.32–1.72)
Richness (Chao1)
ARVI SD	−39.26	86.24	0.00	0.38	0.38	129.45 (79.47–179.23)
EVI2 SD	−39.89	87.50	1.26	0.20	0.58	59.63 (35.33–93.83)
Evenness (PIE)
ARVI SD	23.83	−39.95	0.00	0.28	0.28	8.55 (1.92–28.14)
EVI2 SD	22.99	−38.26	1.69	0.12	0.40	2.47 (1.18–5.18)
SAVI SD	22.89	−38.07	1.89	0.11	0.50	4.01 (1.26–12.74)
NDVI SD	22.89	−38.07	1.89	0.11	0.61	8.02 (1.41–25.44)
Diversity (Shannon entropy)
ARVI SD	−3.67	15.04	0.00	0.33	0.33	13.36 (6.47–20.25)
EVI2 SD	−4.52	16.75	1.71	0.14	0.47	5.94 (2.54–9.34)
Functional richness (FRic)
GNDVI SD	−86.95	181.62	0.00	0.20	0.20	79.91 (13.87–141.98)
SAVI SD	−87.31	182.33	0.71	0.14	0.34	65.32 (9.74–111.22)
NDVI SD	−87.31	182.33	0.71	0.14	0.48	67.24 (12.45–103.76)
ARVI SD	−87.57	182.86	1.24	0.11	0.58	66.64 (7.88–105.43)
Functional evenness (FEve)
ARVI SD	20.10	−32.49	0.00	0.30	0.30	6.16 (1.18–12.74)
EVI2 SD	19.54	−31.37	1.12	0.18	0.48	2.08 (0.85–5.08)
GNDVI Mean	19.51	−31.31	1.18	0.17	0.65	1.63 (0.89–2.97)
Functional diversity (FDiv)
NDMI SD	27.21	−46.71	0.00	0.28	0.28	0.05 (0.01–0.54)
LST Mean	27.06	−46.41	0.30	0.24	0.52	1.04 (1.01–1.06)
LST SD	26.38	−45.05	1.66	0.12	0.65	0.94 (0.88–0.99)
Functional dispersion (FDis)
LST Mean	18.34	−28.97	0.00	0.27	0.27	1.09 (1.02–1.22)
GNDVI Mean	17.97	−28.23	0.74	0.19	0.46	0.41 (0.08–2.19)
EVI2 Mean	17.73	−27.75	1.22	0.15	0.61	0.69 (0.29–1.66)
NDVI Mean	17.65	−27.59	1.38	0.14	0.75	0.53 (0.11–1.66)
SAVI Mean	17.43	−27.15	1.82	0.11	0.86	0.65 (0.22–1.89)
Rao’s quadratic entropy (Q)
GNDVI SD	−39.49	86.69	0.00	0.27	0.27	1.88 (1.31–2.67)
GNDVI Mean	−39.68	87.08	0.39	0.22	0.49	1.61 (1.19–2.83)
LST SD	−39.82	87.35	0.66	0.19	0.68	1.14 (0.09–3.41)
Community evolutionary distinctiveness (CED)
NDBI SD	−8.39	24.49	0.00	0.86	0.86	21.80 (11.27–32.33)
Maximum community evolutionary distinctiveness (maxCED)
NDBI SD	−48.03	103.78	0.00	0.36	0.36	11.8 (5.73–21.34)
NDMI SD	−48.62	104.95	1.17	0.20	0.56	8.15 (2.33–22.82)

Each predictor was modeled separately for each response variable. Abundance models were built under Poisson distributions; diversity models were built under Gaussian distributions. LogLik, log-likelihood; AICc, second-order Akaike information criterion; ΔAICc, AICc difference between the model under concern and the best model; wi, Akaike weight; rate ratio, equivalent to odds ratio in logistic regression. See Table 1 for diversity index codes, Table 2 for spectral index codes, and Table S5 for the full list of models.

, while the full list of models is shown in Table S5. Rate ratios ranged from 0.41 to 129.45, showing a considerable effect of spectral indices on bird abundance and diversity (Table 4). LST SD was the best predictor of abundance, while mean LST was the best predictor of FDis (Table 4). Mean GNDVI, EVI2, NDVI, and SAVI were also parsimonious predictors of FDis. Abundance and FDis increased with increasing LST SD and mean LST, respectively (Figure 3).

ARVI SD was the best predictor of Chao 1 richness, PIE evenness, Shannon entropy, and FEve (Table 4). EVI2 SD was a parsimonious predictor of these diversity indices, while SAVI SD and NDVI SD were parsimonious predictors of PIE evenness, and mean GNDVI was a parsimonious predictor of FEve. Chao 1 richness, PIE evenness, Shannon entropy, and FEve increased with increasing ARVI SD (Figure 3).

GNDVI SD was the best predictor of FRic and Rao’s Q (Table 4). SAVI SD, NDVI SD, and ARVI SD were also parsimonious predictors of FRic, while mean GNDVI and LST SD were parsimonious predictors of Rao’s Q. FRic and Rao’s Q increased with increasing GNDVI SD (Figure 3).

NDBI SD was the best predictor of CED and maxCED (Table 4). NDMI SD was also a parsimonious predictor of maxCED. CED and maxCED increased with increasing NDBI SD (Figure 3).

NDMI SD was the best predictor of FDiv (Table 4). Mean LST and LST SD were also parsimonious predictors of FDiv. FDiv decreased with increasing NDMI SD (Figure 3).

3.3. Spectral Indices as Proxies of Bird Community Structure at the 500 m Spatial Scale

The most parsimonious models for the 500 m buffer zone are shown in

Table 5. The most parsimonious GLM models (ΔAIC_c ≤ 2) selected to assess the relationship of remotely sensed spectral indices, at a 500 m buffer around survey stations, on the abundance and taxonomic, functional, and phylogenetic diversity of the bird community in Kavala.

	LogLik	AICc	ΔAICc	wi	Cumulative wi	Rate Ratio (95% CI)
Abundance
GNDVI Mean	−94.67	194.14	0.00	0.64	0.64	0.05 (0.02–0.14)
Richness (Chao1)
NDMI Mean	−44.38	96.47	0.00	0.24	0.24	74.93 (22.30–127.56)
NDBI Mean	−45.11	97.93	1.46	0.11	0.35	0.23 (0.04–0.88)
ARVI Mean	−45.12	97.95	1.49	0.11	0.46	27.74 (5.36–50.13)
NDVI Mean	−45.23	98.17	1.70	0.10	0.56	26.50 (4.62–48.37)
SAVI Mean	−45.23	98.17	1.70	0.10	0.66	17.67 (3.08–32.25)
Evenness (PIE)
NDMI Mean	21.36	−35.01	0.00	0.18	0.18	2.79 (1.44–7.49)
LST Mean	21.35	−34.99	0.02	0.18	0.36	0.95 (0.89–1.02)
NDBI Mean	20.89	−34.06	0.95	0.11	0.47	0.63 (0.28–1.42)
SAVI Mean	20.85	−33.99	1.02	0.11	0.58	1.23 (0.85–1.78)
NDVI Mean	20.85	−33.99	1.02	0.11	0.69	1.36 (0.78–2.37)
Diversity (Shannon entropy)
NDMI Mean	−7.02	21.74	0.00	0.19	0.19	7.52 (1.28–14.13)
ARVI Mean	−7.69	23.09	1.35	0.10	0.29	14.29 (0.87–74.67)
SAVI Mean	−7.72	23.15	1.41	0.10	0.39	5.52 (0.89–33.93)
NDVI Mean	−7.72	23.15	1.41	0.10	0.49	12.98 (0.85–197.53)
EVI2 Mean	−7.8	23.31	1.56	0.09	0.58	3.94 (0.88–17.54)
Functional richness (FRic)
NDMI Mean	−88.23	184.18	0.00	0.17	0.17	348.44 (33.32–750.22)
LST Mean	−88.47	184.66	0.48	0.14	0.31	0.02 (0.00–35.42)
NDBI Mean	−88.56	184.84	0.65	0.12	0.43	0.01 (0.00–57.84)
ARVI Mean	−88.59	184.89	0.71	0.12	0.55	1.67 (0.65–3.85)
SAVI Mean	−88.63	184.98	0.80	0.12	0.67	1.45 (0.32–5.11)
NDVI Mean	−88.63	184.98	0.80	0.12	0.79	1.12 (0.45–4.38)
Functional evenness (FEve)
GNDVI Mean	18.70	−29.69	0.00	0.18	0.18	2.66 (1.37–3.26)
GNDVI SD	18.60	−29.49	0.20	0.17	0.35	2.43 (1.11–2.58)
EVI2 Mean	18.19	−28.67	1.02	0.11	0.46	1.16 (0.82–1.64)
SAVI Mean	18.19	−28.67	1.02	0.11	0.57	1.19 (0.69–1.82)
NDVI Mean	18.19	−28.67	1.02	0.11	0.68	1.31 (0.69–2.45)
ARVI Mean	17.88	−28.05	1.64	0.08	0.76	1.26 (0.66–2.43)
Functional diversity (FDiv)
NDVI SD	27.49	−47.26	0.00	0.22	0.22	3.45 (1.35–8.81)
SAVI SD	27.49	−47.26	0.00	0.22	0.44	2.28 (1.22–4.27)
EVI2 SD	27.43	−47.14	0.11	0.21	0.65	1.73 (1.14–2.62)
ARVI SD	27.20	−46.69	0.56	0.16	0.81	3.53 (1.51–9.67)
Functional dispersion (FDis)
NDMI Mean	−11.32	30.35	0.00	0.27	0.27	2.48 (0.07–204.35)
NDMI SD	−11.44	30.59	0.24	0.24	0.51	0.00 (0.00–100.71)
LST SD	−11.89	31.49	1.14	0.15	0.66	0.87 (0.61–1.24)
NDBI Mean	−12.15	32.01	1.66	0.12	0.78	0.41 (0.03–5.16)
Rao’s quadratic entropy (Q)
NDMI Mean	−39.76	87.23	0.00	0.23	0.23	15.44 (1.09–67.73)
NDBI Mean	−39.96	87.63	0.40	0.19	0.42	0.01 (0.00–37.32)
ARVI Mean	−40.1	87.92	0.69	0.17	0.59	4.78 (0.11–54.88)
SAVI Mean	−40.59	88.89	1.66	0.10	0.69	2.32 (0.08–57.41)
NDVI Mean	−40.59	88.89	1.66	0.10	0.79	3.20 (0.00–254.65)
Community evolutionary distinctiveness (CED)
LST Mean	−13.9	35.51	0.00	0.25	0.25	0.87 (0.65–0.97)
GNDVI Mean	−14.26	36.24	0.72	0.17	0.42	9.41 (0.09–542.00)
SAVI Mean	−14.33	36.38	0.86	0.16	0.58	3.26 (0.23–44.88)
NDVI Mean	−14.33	36.38	0.86	0.16	0.74	5.89 (0.11–300.45)
Maximum community evolutionary distinctiveness (maxCED)
LST SD	−49.89	107.49	0.00	0.32	0.32	59.08 (2.45–420.77)
LST Mean	−50.07	107.85	0.35	0.27	0.32	1.12 (0.22–3.85)

Each predictor was modeled separately for each response variable. Abundance models were built under Poisson distributions; diversity models were built under Gaussian distributions. LogLik, log-likelihood; AICc, second-order Akaike information criterion; ΔAICc, AICc difference between the model under concern and the best model; wi, Akaike weight; rate ratio, equivalent to odds ratio in logistic regression. See Table 1 for diversity index codes, Table 2 for spectral index codes, and Table S6 for the full list of models.

, while the full list of models is shown in Table S6. Rate ratios ranged from 0.00 to 384.44, showing a considerable effect of spectral indices on bird abundance and diversity (Table 5). Mean GNDVI was the best predictor of abundance and FEve (Table 5). GNDVI SD and mean EVI2, SAVI, NDVI, and ARVI were also parsimonious predictors of FEve. Abundance decreased and FEve increased with increasing mean GNDVI values (Figure 4).

Mean NDMI was the best predictor of Chao 1 richness, PIE evenness, Shannon entropy, FRic, FDis, and Rao’s Q (Table 5). Mean LST was a parsimonious predictor of PIE evenness and FRic, while LST SD was a parsimonious predictor of FDis. Mean NDVI and SAVI were parsimonious predictors of Chao 1 richness, PIE evenness, Shannon entropy, FRic, and Rao’s Q. Mean NDBI was a parsimonious predictor of Chao 1 richness, PIE evenness, FRic, FDis, and Rao’s Q. Mean NDBI was a parsimonious predictor of Chao 1 richness, Shannon entropy, FRic, and Rao’s Q. Mean EVI2 was a parsimonious predictor of Shannon entropy, while NDMI SD was a parsimonious predictor of FDis. Chao 1 richness, PIE evenness, Shannon entropy, FRic, FDis, and Rao’s Q increased with increasing mean NDMI values (Figure 4).

NDVI SD was the best predictor of FDiv (Table 5). SAVI SD, EVI2 SD, and ARVI SD were also parsimonious predictors of NDVI SD. FDiv increased with increasing mean NDMI values (Figure 4).

Mean LST was the best predictor of CED, while LST SD was the best predictor of maxCED. Mean GNDVI, SAVI, and NDVI were also parsimonious predictors of CED. Mean LST was also a parsimonious predictor of maxCED (Table 5). SAVI SD, EVI2 SD, and ARVI SD were also parsimonious predictors of NDVI SD. CED decreased with increasing mean LST values, while maxCED increased with increasing LST SD (Figure 4).

4. Discussion

4.1. Spectral Indices and Bird Abundance

The LST SD was the most important spectral index, positively associated with abundance at the 50 m and 200 m spatial scales. On the other hand, the GNDVI mean was the most important at the 500 m scale, negatively affecting abundance. The urban heat island refers to the higher air and surface temperature in cities compared to adjacent rural areas [66]. This thermal heterogeneity is also observed within cities and between green spaces and built-up areas [67]. Surface temperature differences are higher in built-up areas than in green spaces due to higher heat gains during the day and higher heat losses during the night in the latter. This heterogeneity has been captured by the variation in standard deviations of LST at the 50 m and 200 m spatial scales. Bird abundance increased with increasing temperature heterogeneity, suggesting a higher bird abundance in built-up areas. Previous studies have associated bird abundance with gray infrastructure. Urban dwellers, such as the feral pigeon (Columba livia), collared dove (Streptopelia decaocto), house sparrow (Passer domesticus), and northern house martin (Delichon urbicum) abound in heavily built-up city centers, while forest specialist species, such as the common chaffinch (Fringilla coelebs), European greenfinch (Chloris chloris), Eurasian jay (Garrulus glandarius), great tit (Parus major), and Eurasian blue tit (Cyanistes caeruleus) are found in lower numbers or become practically extinct as impervious cover increases [68,69]. The numerous urban dwellers have been responsible for the increasing abundance with increasing impervious cover in Kavala [4,5,6] and elsewhere [1,21,69]. A vegetation index, the mean GNDVI, was the most important negative predictor of bird abundance at the 500 m scale. Vegetated areas favor forest specialists, which occur in low numbers compared to urban dwellers, resulting in lower abundance with increasing green vegetation cover in Kavala’s green spaces [4,5,6].

4.2. Spectral Indices and Bird Diversity

All spectral vegetation indices were important predictors of taxonomic and functional diversity with their importance varying between indices and spatial scales. Diversity indices increased, with increasing mean vegetation cover and habitat heterogeneity (SDs) at all the spatial scales. NDBI and LST were less important as predictors of taxonomic and functional diversity. On the other hand, NDBI and LST were the most important predictors of phylogenetic diversity, although several vegetation indices were also equally important. Heterogeneity in built-up cover and surface temperature positively influenced phylogenetic diversity, while mean surface temperature had a negative effect. Bird taxonomic, functional, and phylogenetic diversity indices have been commonly associated with high vegetation cover in urban areas, either tree, shrub, or grass [1,4,5,6,21,69,70,71,72,73]. Among spectral indices used to quantify vegetation, NDVI has been the most commonly used proxy of the taxonomic, functional, and phylogenetic diversity of bird communities in urban areas [22,23,24,26,27,28,29,30]. These studies reported a positive relationship between NDVI and all the facets of bird diversity. Other spectral indices, such as EVI and SAVI, have also been used, although in few studies, as proxies of bird diversity in urban areas. [23,74]. These indices were positive predictors of bird diversity. When mean and SD NDVI and EVI values calculated in a 50 m buffer zone around the survey stations were compared, the mean NDVI was the most important predictor of the taxonomic, functional, and phylogenetic diversity of the bird communities in 15 European cities [23].

Our study revealed that different vegetation and urbanization indices were important for predicting bird abundance and diversity in Kavala’s green spaces. Vegetation indices measure primary productivity and are closely related, but they use different spectral bands, algorithms, and corrections. NDVI yields information about vegetation density and health [25]. GNDVI is an index of the plant’s “greenness” or photosynthetic movement, uses the green instead of the red band used in NDVI, and is related to plant vigor [37]. EVI2 is similar to NDVI but corrects for the influence of atmospheric conditions and is more sensitive in densely vegetated areas [31,32]. SAVI incorporates a soil adjustment factor to minimize soil brightness and is particularly useful in areas with sparse vegetation [38,39]. ARVI is an enhanced vegetation index, resistant to atmospheric effects, and useful in areas with high atmospheric variability [40,41]. NDMI measures vegetation water content in plants and is useful in assessing drought conditions and plant stress [42]. NDBI is used to identify urban and built-up areas from satellite imagery [45]. LST measures surface temperature and is useful in climate studies and urban heat island effect analysis and monitoring [49]. Although similar, different spectral indices have differing sensitivity in measuring the different aspects, quality, and configuration of the vegetated and non-vegetated areas in cities. Urban areas are very heterogeneous both at the local and landscape scales, and this heterogeneity can be differently captured by different spectral indices [1,2].

Kavala’s green spaces consist of woodlands, square gardens, playgrounds, and median road strips [4,5,6]. Besides one 18 ha woodland (survey station 7), all other green spaces were smaller than 3 ha, irregular in size, and irregularly dispersed among built-up areas [4]. The woodlands have a Turkish pine (Pinus brutia)-dominated tree storey and a shrub understorey mainly consisting of lentisk (Pistacia lentiscus), kermes oak (Quercus coccifera), the green olive tree (Phillyrea latifolia), and flowering ash (Fraxinus ornus). The other green spaces consisted of both native (e.g., silver lime (Tilia tomentosa), bay laurel (Laurus nobilis), box (Buxus sempervirens), and tamarisk (Tamarix sp.)) and non-native plants (e.g., Chinaberry (Melia azedarach), oriental arborvitae (Platycladus orientalis), paper mulberry (Brussonetia papyrifera), and Canary Island date palm (Phoenix canariensis)). Most green spaces include paved walks, cafés, and restaurants. Also, all the buffer zones used included buildings surrounding the green spaces.

The degree of habitat fragmentation at various landscape levels and the size of habitat patches are important determinants of bird species abundance, richness, and diversity [75]. The green spaces of Kavala were small fragments of vegetation, irregularly spaced among built-up areas [4,5,6]. This variable mix of green and gray infrastructure, in combination with the differences among spectral indices, might explain the varying importance of these indices in predicting bird abundance, richness, and diversity, both within and between different spatial scales. EVI2 would perform better in areas with denser vegetation, while SAVI would more reliably measure areas with sparser vegetation. ARVI would perform better in areas with high atmospheric variability, NDBI would better separate green from built-up areas, and LST would better capture big differences in temperature among areas, especially in city centers.

Mean values of spectral indices and habitat heterogeneity, implied by SDs, were important predictors of diversity indices in Kavala. Areas with higher plant biomass and habitat heterogeneity supported more species were more diverse, both taxonomically and functionally, and hosted more evolutionarily distinct species. Benedetti et al. [23] reported that means were more important predictors of bird diversity than SDs of NDVI and EVI in 15 European cities. SD values can capture habitat heterogeneity [76]. SDs of spectral indices were important predictors of most facets of bird diversity in Kavala because heterogeneous habitats contain a higher variety of microhabitats and resources potentially suitable for many bird species [77,78]. In agreement with Benedetti et al. [44], NDVI was a better predictor of bird diversity in Kavala than EVI2 in most cases. However, we also used other spectral indices that were more important than NDVI and EVI2 in predicting several facets of bird diversity. Also, we used EVI2, which is an improvement on EVI used in Benedetti et al. [23], as it better assesses dense vegetation [31,32].

5. Conclusions

We investigated the ability of eight vegetation and urbanization spectral indices to predict the abundance, species richness, and taxonomic, functional, and phylogenetic diversity of the bird community of Kavala. Our results showed that most indices could be used as proxies of the bird community structure. However, different indices performed better for different facets of bird diversity and at different spatial scales. This trend could be attributed to the differences among spectral indices, each of them being sensitive to different spectra and vegetation structures, the different vegetation types, health, and density among green spaces, and the different mixes of green and gray infrastructure among spatial scales. The performance of these spectral indices should also be investigated in other urban areas of various sizes and spatial configurations of green and built-up areas. The main outcome of our study was that researchers should use several easy- and cheap-to-measure spectral indices to determine the most suitable for predicting urban bird community structure. The use of easy- and cheap-to-measure spectral indices will provide a valuable tool to conservationists and urban planners. They will be able to readily infer the status of urban bird communities and design and implement successful management plans to maintain and enhance urban biodiversity, a key element for ecosystem stability and human health [11,18,19,20].